lemon-project-template-glpk

comparison deps/glpk/src/glpnpp06.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:146ef40247a9
1 /* glpnpp06.c (translate feasibility problem to CNF-SAT) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
24
25 #include "glpnpp.h"
26
27 /***********************************************************************
28 * npp_sat_free_row - process free (unbounded) row
29 *
30 * This routine processes row p, which is free (i.e. has no finite
31 * bounds):
32 *
33 * -inf < sum a[p,j] x[j] < +inf. (1)
34 *
35 * The constraint (1) cannot be active and therefore it is redundant,
36 * so the routine simply removes it from the original problem. */
37
38 void npp_sat_free_row(NPP *npp, NPPROW *p)
39 { /* the row should be free */
40 xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX);
41 /* remove the row from the problem */
42 npp_del_row(npp, p);
43 return;
44 }
45
46 /***********************************************************************
47 * npp_sat_fixed_col - process fixed column
48 *
49 * This routine processes column q, which is fixed:
50 *
51 * x[q] = s[q], (1)
52 *
53 * where s[q] is a fixed column value.
54 *
55 * The routine substitutes fixed value s[q] into constraint rows and
56 * then removes column x[q] from the original problem.
57 *
58 * Substitution of x[q] = s[q] into row i gives:
59 *
60 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
61 * j
62 *
63 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
64 * j!=q
65 *
66 * L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==>
67 * j!=q
68 *
69 * L~[i] <= sum a[i,j] x[j] <= U~[i],
70 * j!=q
71 *
72 * where
73 *
74 * L~[i] = L[i] - a[i,q] s[q], (2)
75 *
76 * U~[i] = U[i] - a[i,q] s[q] (3)
77 *
78 * are, respectively, lower and upper bound of row i in the transformed
79 * problem.
80 *
81 * On recovering solution x[q] is assigned the value of s[q]. */
82
83 struct sat_fixed_col
84 { /* fixed column */
85 int q;
86 /* column reference number for variable x[q] */
87 int s;
88 /* value, at which x[q] is fixed */
89 };
90
91 static int rcv_sat_fixed_col(NPP *, void *);
92
93 int npp_sat_fixed_col(NPP *npp, NPPCOL *q)
94 { struct sat_fixed_col *info;
95 NPPROW *i;
96 NPPAIJ *aij;
97 int temp;
98 /* the column should be fixed */
99 xassert(q->lb == q->ub);
100 /* create transformation stack entry */
101 info = npp_push_tse(npp,
102 rcv_sat_fixed_col, sizeof(struct sat_fixed_col));
103 info->q = q->j;
104 info->s = (int)q->lb;
105 xassert((double)info->s == q->lb);
106 /* substitute x[q] = s[q] into constraint rows */
107 if (info->s == 0)
108 goto skip;
109 for (aij = q->ptr; aij != NULL; aij = aij->c_next)
110 { i = aij->row;
111 if (i->lb != -DBL_MAX)
112 { i->lb -= aij->val * (double)info->s;
113 temp = (int)i->lb;
114 if ((double)temp != i->lb)
115 return 1; /* integer arithmetic error */
116 }
117 if (i->ub != +DBL_MAX)
118 { i->ub -= aij->val * (double)info->s;
119 temp = (int)i->ub;
120 if ((double)temp != i->ub)
121 return 2; /* integer arithmetic error */
122 }
123 }
124 skip: /* remove the column from the problem */
125 npp_del_col(npp, q);
126 return 0;
127 }
128
129 static int rcv_sat_fixed_col(NPP *npp, void *info_)
130 { struct sat_fixed_col *info = info_;
131 npp->c_value[info->q] = (double)info->s;
132 return 0;
133 }
134
135 /***********************************************************************
136 * npp_sat_is_bin_comb - test if row is binary combination
137 *
138 * This routine tests if the specified row is a binary combination,
139 * i.e. all its constraint coefficients are +1 and -1 and all variables
140 * are binary. If the test was passed, the routine returns non-zero,
141 * otherwise zero. */
142
143 int npp_sat_is_bin_comb(NPP *npp, NPPROW *row)
144 { NPPCOL *col;
145 NPPAIJ *aij;
146 xassert(npp == npp);
147 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
148 { if (!(aij->val == +1.0 || aij->val == -1.0))
149 return 0; /* non-unity coefficient */
150 col = aij->col;
151 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
152 return 0; /* non-binary column */
153 }
154 return 1; /* test was passed */
155 }
156
157 /***********************************************************************
158 * npp_sat_num_pos_coef - determine number of positive coefficients
159 *
160 * This routine returns the number of positive coefficients in the
161 * specified row. */
162
163 int npp_sat_num_pos_coef(NPP *npp, NPPROW *row)
164 { NPPAIJ *aij;
165 int num = 0;
166 xassert(npp == npp);
167 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
168 { if (aij->val > 0.0)
169 num++;
170 }
171 return num;
172 }
173
174 /***********************************************************************
175 * npp_sat_num_neg_coef - determine number of negative coefficients
176 *
177 * This routine returns the number of negative coefficients in the
178 * specified row. */
179
180 int npp_sat_num_neg_coef(NPP *npp, NPPROW *row)
181 { NPPAIJ *aij;
182 int num = 0;
183 xassert(npp == npp);
184 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
185 { if (aij->val < 0.0)
186 num++;
187 }
188 return num;
189 }
190
191 /***********************************************************************
192 * npp_sat_is_cover_ineq - test if row is covering inequality
193 *
194 * The canonical form of a covering inequality is the following:
195 *
196 * sum x[j] >= 1, (1)
197 * j in J
198 *
199 * where all x[j] are binary variables.
200 *
201 * In general case a covering inequality may have one of the following
202 * two forms:
203 *
204 * sum x[j] - sum x[j] >= 1 - |J-|, (2)
205 * j in J+ j in J-
206 *
207 *
208 * sum x[j] - sum x[j] <= |J+| - 1. (3)
209 * j in J+ j in J-
210 *
211 * Obviously, the inequality (2) can be transformed to the form (1) by
212 * substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
213 * negation of variable x[j]. And the inequality (3) can be transformed
214 * to (2) by multiplying both left- and right-hand sides by -1.
215 *
216 * This routine returns one of the following codes:
217 *
218 * 0, if the specified row is not a covering inequality;
219 *
220 * 1, if the specified row has the form (2);
221 *
222 * 2, if the specified row has the form (3). */
223
224 int npp_sat_is_cover_ineq(NPP *npp, NPPROW *row)
225 { xassert(npp == npp);
226 if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
227 { /* row is inequality of '>=' type */
228 if (npp_sat_is_bin_comb(npp, row))
229 { /* row is a binary combination */
230 if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
231 { /* row has the form (2) */
232 return 1;
233 }
234 }
235 }
236 else if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
237 { /* row is inequality of '<=' type */
238 if (npp_sat_is_bin_comb(npp, row))
239 { /* row is a binary combination */
240 if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
241 { /* row has the form (3) */
242 return 2;
243 }
244 }
245 }
246 /* row is not a covering inequality */
247 return 0;
248 }
249
250 /***********************************************************************
251 * npp_sat_is_pack_ineq - test if row is packing inequality
252 *
253 * The canonical form of a packing inequality is the following:
254 *
255 * sum x[j] <= 1, (1)
256 * j in J
257 *
258 * where all x[j] are binary variables.
259 *
260 * In general case a packing inequality may have one of the following
261 * two forms:
262 *
263 * sum x[j] - sum x[j] <= 1 - |J-|, (2)
264 * j in J+ j in J-
265 *
266 *
267 * sum x[j] - sum x[j] >= |J+| - 1. (3)
268 * j in J+ j in J-
269 *
270 * Obviously, the inequality (2) can be transformed to the form (1) by
271 * substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
272 * negation of variable x[j]. And the inequality (3) can be transformed
273 * to (2) by multiplying both left- and right-hand sides by -1.
274 *
275 * This routine returns one of the following codes:
276 *
277 * 0, if the specified row is not a packing inequality;
278 *
279 * 1, if the specified row has the form (2);
280 *
281 * 2, if the specified row has the form (3). */
282
283 int npp_sat_is_pack_ineq(NPP *npp, NPPROW *row)
284 { xassert(npp == npp);
285 if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
286 { /* row is inequality of '<=' type */
287 if (npp_sat_is_bin_comb(npp, row))
288 { /* row is a binary combination */
289 if (row->ub == 1.0 - npp_sat_num_neg_coef(npp, row))
290 { /* row has the form (2) */
291 return 1;
292 }
293 }
294 }
295 else if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
296 { /* row is inequality of '>=' type */
297 if (npp_sat_is_bin_comb(npp, row))
298 { /* row is a binary combination */
299 if (row->lb == npp_sat_num_pos_coef(npp, row) - 1.0)
300 { /* row has the form (3) */
301 return 2;
302 }
303 }
304 }
305 /* row is not a packing inequality */
306 return 0;
307 }
308
309 /***********************************************************************
310 * npp_sat_is_partn_eq - test if row is partitioning equality
311 *
312 * The canonical form of a partitioning equality is the following:
313 *
314 * sum x[j] = 1, (1)
315 * j in J
316 *
317 * where all x[j] are binary variables.
318 *
319 * In general case a partitioning equality may have one of the following
320 * two forms:
321 *
322 * sum x[j] - sum x[j] = 1 - |J-|, (2)
323 * j in J+ j in J-
324 *
325 *
326 * sum x[j] - sum x[j] = |J+| - 1. (3)
327 * j in J+ j in J-
328 *
329 * Obviously, the equality (2) can be transformed to the form (1) by
330 * substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
331 * negation of variable x[j]. And the equality (3) can be transformed
332 * to (2) by multiplying both left- and right-hand sides by -1.
333 *
334 * This routine returns one of the following codes:
335 *
336 * 0, if the specified row is not a partitioning equality;
337 *
338 * 1, if the specified row has the form (2);
339 *
340 * 2, if the specified row has the form (3). */
341
342 int npp_sat_is_partn_eq(NPP *npp, NPPROW *row)
343 { xassert(npp == npp);
344 if (row->lb == row->ub)
345 { /* row is equality constraint */
346 if (npp_sat_is_bin_comb(npp, row))
347 { /* row is a binary combination */
348 if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
349 { /* row has the form (2) */
350 return 1;
351 }
352 if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
353 { /* row has the form (3) */
354 return 2;
355 }
356 }
357 }
358 /* row is not a partitioning equality */
359 return 0;
360 }
361
362 /***********************************************************************
363 * npp_sat_reverse_row - multiply both sides of row by -1
364 *
365 * This routines multiplies by -1 both left- and right-hand sides of
366 * the specified row:
367 *
368 * L <= sum x[j] <= U,
369 *
370 * that results in the following row:
371 *
372 * -U <= sum (-x[j]) <= -L.
373 *
374 * If no integer overflow occured, the routine returns zero, otherwise
375 * non-zero. */
376
377 int npp_sat_reverse_row(NPP *npp, NPPROW *row)
378 { NPPAIJ *aij;
379 int temp, ret = 0;
380 double old_lb, old_ub;
381 xassert(npp == npp);
382 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
383 { aij->val = -aij->val;
384 temp = (int)aij->val;
385 if ((double)temp != aij->val)
386 ret = 1;
387 }
388 old_lb = row->lb, old_ub = row->ub;
389 if (old_ub == +DBL_MAX)
390 row->lb = -DBL_MAX;
391 else
392 { row->lb = -old_ub;
393 temp = (int)row->lb;
394 if ((double)temp != row->lb)
395 ret = 2;
396 }
397 if (old_lb == -DBL_MAX)
398 row->ub = +DBL_MAX;
399 else
400 { row->ub = -old_lb;
401 temp = (int)row->ub;
402 if ((double)temp != row->ub)
403 ret = 3;
404 }
405 return ret;
406 }
407
408 /***********************************************************************
409 * npp_sat_split_pack - split packing inequality
410 *
411 * Let there be given a packing inequality in canonical form:
412 *
413 * sum t[j] <= 1, (1)
414 * j in J
415 *
416 * where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable.
417 * And let J = J1 U J2 is a partition of the set of literals. Then the
418 * inequality (1) is obviously equivalent to the following two packing
419 * inequalities:
420 *
421 * sum t[j] <= y <--> sum t[j] + (1 - y) <= 1, (2)
422 * j in J1 j in J1
423 *
424 * sum t[j] <= 1 - y <--> sum t[j] + y <= 1, (3)
425 * j in J2 j in J2
426 *
427 * where y is a new binary variable added to the transformed problem.
428 *
429 * Assuming that the specified row is a packing inequality (1), this
430 * routine constructs the set J1 by including there first nlit literals
431 * (terms) from the specified row, and the set J2 = J \ J1. Then the
432 * routine creates a new row, which corresponds to inequality (2), and
433 * replaces the specified row with inequality (3). */
434
435 NPPROW *npp_sat_split_pack(NPP *npp, NPPROW *row, int nlit)
436 { NPPROW *rrr;
437 NPPCOL *col;
438 NPPAIJ *aij;
439 int k;
440 /* original row should be packing inequality (1) */
441 xassert(npp_sat_is_pack_ineq(npp, row) == 1);
442 /* and nlit should be less than the number of literals (terms)
443 in the original row */
444 xassert(0 < nlit && nlit < npp_row_nnz(npp, row));
445 /* create new row corresponding to inequality (2) */
446 rrr = npp_add_row(npp);
447 rrr->lb = -DBL_MAX, rrr->ub = 1.0;
448 /* move first nlit literals (terms) from the original row to the
449 new row; the original row becomes inequality (3) */
450 for (k = 1; k <= nlit; k++)
451 { aij = row->ptr;
452 xassert(aij != NULL);
453 /* add literal to the new row */
454 npp_add_aij(npp, rrr, aij->col, aij->val);
455 /* correct rhs */
456 if (aij->val < 0.0)
457 rrr->ub -= 1.0, row->ub += 1.0;
458 /* remove literal from the original row */
459 npp_del_aij(npp, aij);
460 }
461 /* create new binary variable y */
462 col = npp_add_col(npp);
463 col->is_int = 1, col->lb = 0.0, col->ub = 1.0;
464 /* include literal (1 - y) in the new row */
465 npp_add_aij(npp, rrr, col, -1.0);
466 rrr->ub -= 1.0;
467 /* include literal y in the original row */
468 npp_add_aij(npp, row, col, +1.0);
469 return rrr;
470 }
471
472 /***********************************************************************
473 * npp_sat_encode_pack - encode packing inequality
474 *
475 * Given a packing inequality in canonical form:
476 *
477 * sum t[j] <= 1, (1)
478 * j in J
479 *
480 * where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable,
481 * this routine translates it to CNF by replacing it with the following
482 * equivalent set of edge packing inequalities:
483 *
484 * t[j] + t[k] <= 1 for all j, k in J, j != k. (2)
485 *
486 * Then the routine transforms each edge packing inequality (2) to
487 * corresponding covering inequality (that encodes two-literal clause)
488 * by multiplying both its part by -1:
489 *
490 * - t[j] - t[k] >= -1 <--> (1 - t[j]) + (1 - t[k]) >= 1. (3)
491 *
492 * On exit the routine removes the original row from the problem. */
493
494 void npp_sat_encode_pack(NPP *npp, NPPROW *row)
495 { NPPROW *rrr;
496 NPPAIJ *aij, *aik;
497 /* original row should be packing inequality (1) */
498 xassert(npp_sat_is_pack_ineq(npp, row) == 1);
499 /* create equivalent system of covering inequalities (3) */
500 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
501 { /* due to symmetry only one of inequalities t[j] + t[k] <= 1
502 and t[k] <= t[j] <= 1 can be considered */
503 for (aik = aij->r_next; aik != NULL; aik = aik->r_next)
504 { /* create edge packing inequality (2) */
505 rrr = npp_add_row(npp);
506 rrr->lb = -DBL_MAX, rrr->ub = 1.0;
507 npp_add_aij(npp, rrr, aij->col, aij->val);
508 if (aij->val < 0.0)
509 rrr->ub -= 1.0;
510 npp_add_aij(npp, rrr, aik->col, aik->val);
511 if (aik->val < 0.0)
512 rrr->ub -= 1.0;
513 /* and transform it to covering inequality (3) */
514 npp_sat_reverse_row(npp, rrr);
515 xassert(npp_sat_is_cover_ineq(npp, rrr) == 1);
516 }
517 }
518 /* remove the original row from the problem */
519 npp_del_row(npp, row);
520 return;
521 }
522
523 /***********************************************************************
524 * npp_sat_encode_sum2 - encode 2-bit summation
525 *
526 * Given a set containing two literals x and y this routine encodes
527 * the equality
528 *
529 * x + y = s + 2 * c, (1)
530 *
531 * where
532 *
533 * s = (x + y) % 2 (2)
534 *
535 * is a binary variable modeling the low sum bit, and
536 *
537 * c = (x + y) / 2 (3)
538 *
539 * is a binary variable modeling the high (carry) sum bit. */
540
541 void npp_sat_encode_sum2(NPP *npp, NPPLSE *set, NPPSED *sed)
542 { NPPROW *row;
543 int x, y, s, c;
544 /* the set should contain exactly two literals */
545 xassert(set != NULL);
546 xassert(set->next != NULL);
547 xassert(set->next->next == NULL);
548 sed->x = set->lit;
549 xassert(sed->x.neg == 0 || sed->x.neg == 1);
550 sed->y = set->next->lit;
551 xassert(sed->y.neg == 0 || sed->y.neg == 1);
552 sed->z.col = NULL, sed->z.neg = 0;
553 /* perform encoding s = (x + y) % 2 */
554 sed->s = npp_add_col(npp);
555 sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
556 for (x = 0; x <= 1; x++)
557 { for (y = 0; y <= 1; y++)
558 { for (s = 0; s <= 1; s++)
559 { if ((x + y) % 2 != s)
560 { /* generate CNF clause to disable infeasible
561 combination */
562 row = npp_add_row(npp);
563 row->lb = 1.0, row->ub = +DBL_MAX;
564 if (x == sed->x.neg)
565 npp_add_aij(npp, row, sed->x.col, +1.0);
566 else
567 { npp_add_aij(npp, row, sed->x.col, -1.0);
568 row->lb -= 1.0;
569 }
570 if (y == sed->y.neg)
571 npp_add_aij(npp, row, sed->y.col, +1.0);
572 else
573 { npp_add_aij(npp, row, sed->y.col, -1.0);
574 row->lb -= 1.0;
575 }
576 if (s == 0)
577 npp_add_aij(npp, row, sed->s, +1.0);
578 else
579 { npp_add_aij(npp, row, sed->s, -1.0);
580 row->lb -= 1.0;
581 }
582 }
583 }
584 }
585 }
586 /* perform encoding c = (x + y) / 2 */
587 sed->c = npp_add_col(npp);
588 sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
589 for (x = 0; x <= 1; x++)
590 { for (y = 0; y <= 1; y++)
591 { for (c = 0; c <= 1; c++)
592 { if ((x + y) / 2 != c)
593 { /* generate CNF clause to disable infeasible
594 combination */
595 row = npp_add_row(npp);
596 row->lb = 1.0, row->ub = +DBL_MAX;
597 if (x == sed->x.neg)
598 npp_add_aij(npp, row, sed->x.col, +1.0);
599 else
600 { npp_add_aij(npp, row, sed->x.col, -1.0);
601 row->lb -= 1.0;
602 }
603 if (y == sed->y.neg)
604 npp_add_aij(npp, row, sed->y.col, +1.0);
605 else
606 { npp_add_aij(npp, row, sed->y.col, -1.0);
607 row->lb -= 1.0;
608 }
609 if (c == 0)
610 npp_add_aij(npp, row, sed->c, +1.0);
611 else
612 { npp_add_aij(npp, row, sed->c, -1.0);
613 row->lb -= 1.0;
614 }
615 }
616 }
617 }
618 }
619 return;
620 }
621
622 /***********************************************************************
623 * npp_sat_encode_sum3 - encode 3-bit summation
624 *
625 * Given a set containing at least three literals this routine chooses
626 * some literals x, y, z from that set and encodes the equality
627 *
628 * x + y + z = s + 2 * c, (1)
629 *
630 * where
631 *
632 * s = (x + y + z) % 2 (2)
633 *
634 * is a binary variable modeling the low sum bit, and
635 *
636 * c = (x + y + z) / 2 (3)
637 *
638 * is a binary variable modeling the high (carry) sum bit. */
639
640 void npp_sat_encode_sum3(NPP *npp, NPPLSE *set, NPPSED *sed)
641 { NPPROW *row;
642 int x, y, z, s, c;
643 /* the set should contain at least three literals */
644 xassert(set != NULL);
645 xassert(set->next != NULL);
646 xassert(set->next->next != NULL);
647 sed->x = set->lit;
648 xassert(sed->x.neg == 0 || sed->x.neg == 1);
649 sed->y = set->next->lit;
650 xassert(sed->y.neg == 0 || sed->y.neg == 1);
651 sed->z = set->next->next->lit;
652 xassert(sed->z.neg == 0 || sed->z.neg == 1);
653 /* perform encoding s = (x + y + z) % 2 */
654 sed->s = npp_add_col(npp);
655 sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
656 for (x = 0; x <= 1; x++)
657 { for (y = 0; y <= 1; y++)
658 { for (z = 0; z <= 1; z++)
659 { for (s = 0; s <= 1; s++)
660 { if ((x + y + z) % 2 != s)
661 { /* generate CNF clause to disable infeasible
662 combination */
663 row = npp_add_row(npp);
664 row->lb = 1.0, row->ub = +DBL_MAX;
665 if (x == sed->x.neg)
666 npp_add_aij(npp, row, sed->x.col, +1.0);
667 else
668 { npp_add_aij(npp, row, sed->x.col, -1.0);
669 row->lb -= 1.0;
670 }
671 if (y == sed->y.neg)
672 npp_add_aij(npp, row, sed->y.col, +1.0);
673 else
674 { npp_add_aij(npp, row, sed->y.col, -1.0);
675 row->lb -= 1.0;
676 }
677 if (z == sed->z.neg)
678 npp_add_aij(npp, row, sed->z.col, +1.0);
679 else
680 { npp_add_aij(npp, row, sed->z.col, -1.0);
681 row->lb -= 1.0;
682 }
683 if (s == 0)
684 npp_add_aij(npp, row, sed->s, +1.0);
685 else
686 { npp_add_aij(npp, row, sed->s, -1.0);
687 row->lb -= 1.0;
688 }
689 }
690 }
691 }
692 }
693 }
694 /* perform encoding c = (x + y + z) / 2 */
695 sed->c = npp_add_col(npp);
696 sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
697 for (x = 0; x <= 1; x++)
698 { for (y = 0; y <= 1; y++)
699 { for (z = 0; z <= 1; z++)
700 { for (c = 0; c <= 1; c++)
701 { if ((x + y + z) / 2 != c)
702 { /* generate CNF clause to disable infeasible
703 combination */
704 row = npp_add_row(npp);
705 row->lb = 1.0, row->ub = +DBL_MAX;
706 if (x == sed->x.neg)
707 npp_add_aij(npp, row, sed->x.col, +1.0);
708 else
709 { npp_add_aij(npp, row, sed->x.col, -1.0);
710 row->lb -= 1.0;
711 }
712 if (y == sed->y.neg)
713 npp_add_aij(npp, row, sed->y.col, +1.0);
714 else
715 { npp_add_aij(npp, row, sed->y.col, -1.0);
716 row->lb -= 1.0;
717 }
718 if (z == sed->z.neg)
719 npp_add_aij(npp, row, sed->z.col, +1.0);
720 else
721 { npp_add_aij(npp, row, sed->z.col, -1.0);
722 row->lb -= 1.0;
723 }
724 if (c == 0)
725 npp_add_aij(npp, row, sed->c, +1.0);
726 else
727 { npp_add_aij(npp, row, sed->c, -1.0);
728 row->lb -= 1.0;
729 }
730 }
731 }
732 }
733 }
734 }
735 return;
736 }
737
738 /***********************************************************************
739 * npp_sat_encode_sum_ax - encode linear combination of 0-1 variables
740 *
741 * PURPOSE
742 *
743 * Given a linear combination of binary variables:
744 *
745 * sum a[j] x[j], (1)
746 * j
747 *
748 * which is the linear form of the specified row, this routine encodes
749 * (i.e. translates to CNF) the following equality:
750 *
751 * n
752 * sum |a[j]| t[j] = sum 2**(k-1) * y[k], (2)
753 * j k=1
754 *
755 * where t[j] = x[j] (if a[j] > 0) or t[j] = 1 - x[j] (if a[j] < 0),
756 * and y[k] is either t[j] or a new literal created by the routine or
757 * a constant zero. Note that the sum in the right-hand side of (2) can
758 * be thought as a n-bit representation of the sum in the left-hand
759 * side, which is a non-negative integer number.
760 *
761 * ALGORITHM
762 *
763 * First, the number of bits, n, sufficient to represent any value in
764 * the left-hand side of (2) is determined. Obviously, n is the number
765 * of bits sufficient to represent the sum (sum |a[j]|).
766 *
767 * Let
768 *
769 * n
770 * |a[j]| = sum 2**(k-1) b[j,k], (3)
771 * k=1
772 *
773 * where b[j,k] is k-th bit in a n-bit representation of |a[j]|. Then
774 *
775 * m n
776 * sum |a[j]| * t[j] = sum 2**(k-1) sum b[j,k] * t[j]. (4)
777 * j k=1 j=1
778 *
779 * Introducing the set
780 *
781 * J[k] = { j : b[j,k] = 1 } (5)
782 *
783 * allows rewriting (4) as follows:
784 *
785 * n
786 * sum |a[j]| * t[j] = sum 2**(k-1) sum t[j]. (6)
787 * j k=1 j in J[k]
788 *
789 * Thus, our goal is to provide |J[k]| <= 1 for all k, in which case
790 * we will have the representation (1).
791 *
792 * Let |J[k]| = 2, i.e. J[k] has exactly two literals u and v. In this
793 * case we can apply the following transformation:
794 *
795 * u + v = s + 2 * c, (7)
796 *
797 * where s and c are, respectively, low (sum) and high (carry) bits of
798 * the sum of two bits. This allows to replace two literals u and v in
799 * J[k] by one literal s, and carry out literal c to J[k+1].
800 *
801 * If |J[k]| >= 3, i.e. J[k] has at least three literals u, v, and w,
802 * we can apply the following transformation:
803 *
804 * u + v + w = s + 2 * c. (8)
805 *
806 * Again, literal s replaces literals u, v, and w in J[k], and literal
807 * c goes into J[k+1].
808 *
809 * On exit the routine stores each literal from J[k] in element y[k],
810 * 1 <= k <= n. If J[k] is empty, y[k] is set to constant false.
811 *
812 * RETURNS
813 *
814 * The routine returns n, the number of literals in the right-hand side
815 * of (2), 0 <= n <= NBIT_MAX. If the sum (sum |a[j]|) is too large, so
816 * more than NBIT_MAX (= 31) literals are needed to encode the original
817 * linear combination, the routine returns a negative value. */
818
819 #define NBIT_MAX 31
820 /* maximal number of literals in the right hand-side of (2) */
821
822 static NPPLSE *remove_lse(NPP *npp, NPPLSE *set, NPPCOL *col)
823 { /* remove specified literal from specified literal set */
824 NPPLSE *lse, *prev = NULL;
825 for (lse = set; lse != NULL; prev = lse, lse = lse->next)
826 if (lse->lit.col == col) break;
827 xassert(lse != NULL);
828 if (prev == NULL)
829 set = lse->next;
830 else
831 prev->next = lse->next;
832 dmp_free_atom(npp->pool, lse, sizeof(NPPLSE));
833 return set;
834 }
835
836 int npp_sat_encode_sum_ax(NPP *npp, NPPROW *row, NPPLIT y[])
837 { NPPAIJ *aij;
838 NPPLSE *set[1+NBIT_MAX], *lse;
839 NPPSED sed;
840 int k, n, temp;
841 double sum;
842 /* compute the sum (sum |a[j]|) */
843 sum = 0.0;
844 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
845 sum += fabs(aij->val);
846 /* determine n, the number of bits in the sum */
847 temp = (int)sum;
848 if ((double)temp != sum)
849 return -1; /* integer arithmetic error */
850 for (n = 0; temp > 0; n++, temp >>= 1);
851 xassert(0 <= n && n <= NBIT_MAX);
852 /* build initial sets J[k], 1 <= k <= n; see (5) */
853 /* set[k] is a pointer to the list of literals in J[k] */
854 for (k = 1; k <= n; k++)
855 set[k] = NULL;
856 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
857 { temp = (int)fabs(aij->val);
858 xassert((int)temp == fabs(aij->val));
859 for (k = 1; temp > 0; k++, temp >>= 1)
860 { if (temp & 1)
861 { xassert(k <= n);
862 lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
863 lse->lit.col = aij->col;
864 lse->lit.neg = (aij->val > 0.0 ? 0 : 1);
865 lse->next = set[k];
866 set[k] = lse;
867 }
868 }
869 }
870 /* main transformation loop */
871 for (k = 1; k <= n; k++)
872 { /* reduce J[k] and set y[k] */
873 for (;;)
874 { if (set[k] == NULL)
875 { /* J[k] is empty */
876 /* set y[k] to constant false */
877 y[k].col = NULL, y[k].neg = 0;
878 break;
879 }
880 if (set[k]->next == NULL)
881 { /* J[k] contains one literal */
882 /* set y[k] to that literal */
883 y[k] = set[k]->lit;
884 dmp_free_atom(npp->pool, set[k], sizeof(NPPLSE));
885 break;
886 }
887 if (set[k]->next->next == NULL)
888 { /* J[k] contains two literals */
889 /* apply transformation (7) */
890 npp_sat_encode_sum2(npp, set[k], &sed);
891 }
892 else
893 { /* J[k] contains at least three literals */
894 /* apply transformation (8) */
895 npp_sat_encode_sum3(npp, set[k], &sed);
896 /* remove third literal from set[k] */
897 set[k] = remove_lse(npp, set[k], sed.z.col);
898 }
899 /* remove second literal from set[k] */
900 set[k] = remove_lse(npp, set[k], sed.y.col);
901 /* remove first literal from set[k] */
902 set[k] = remove_lse(npp, set[k], sed.x.col);
903 /* include new literal s to set[k] */
904 lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
905 lse->lit.col = sed.s, lse->lit.neg = 0;
906 lse->next = set[k];
907 set[k] = lse;
908 /* include new literal c to set[k+1] */
909 xassert(k < n); /* FIXME: can "overflow" happen? */
910 lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
911 lse->lit.col = sed.c, lse->lit.neg = 0;
912 lse->next = set[k+1];
913 set[k+1] = lse;
914 }
915 }
916 return n;
917 }
918
919 /***********************************************************************
920 * npp_sat_normalize_clause - normalize clause
921 *
922 * This routine normalizes the specified clause, which is a disjunction
923 * of literals, by replacing multiple literals, which refer to the same
924 * binary variable, with a single literal.
925 *
926 * On exit the routine returns the number of literals in the resulting
927 * clause. However, if the specified clause includes both a literal and
928 * its negation, the routine returns a negative value meaning that the
929 * clause is equivalent to the value true. */
930
931 int npp_sat_normalize_clause(NPP *npp, int size, NPPLIT lit[])
932 { int j, k, new_size;
933 xassert(npp == npp);
934 xassert(size >= 0);
935 new_size = 0;
936 for (k = 1; k <= size; k++)
937 { for (j = 1; j <= new_size; j++)
938 { if (lit[k].col == lit[j].col)
939 { /* lit[k] refers to the same variable as lit[j], which
940 is already included in the resulting clause */
941 if (lit[k].neg == lit[j].neg)
942 { /* ignore lit[k] due to the idempotent law */
943 goto skip;
944 }
945 else
946 { /* lit[k] is NOT lit[j]; the clause is equivalent to
947 the value true */
948 return -1;
949 }
950 }
951 }
952 /* include lit[k] in the resulting clause */
953 lit[++new_size] = lit[k];
954 skip: ;
955 }
956 return new_size;
957 }
958
959 /***********************************************************************
960 * npp_sat_encode_clause - translate clause to cover inequality
961 *
962 * Given a clause
963 *
964 * OR t[j], (1)
965 * j in J
966 *
967 * where t[j] is a literal, i.e. t[j] = x[j] or t[j] = NOT x[j], this
968 * routine translates it to the following equivalent cover inequality,
969 * which is added to the transformed problem:
970 *
971 * sum t[j] >= 1, (2)
972 * j in J
973 *
974 * where t[j] = x[j] or t[j] = 1 - x[j].
975 *
976 * If necessary, the clause should be normalized before a call to this
977 * routine. */
978
979 NPPROW *npp_sat_encode_clause(NPP *npp, int size, NPPLIT lit[])
980 { NPPROW *row;
981 int k;
982 xassert(size >= 1);
983 row = npp_add_row(npp);
984 row->lb = 1.0, row->ub = +DBL_MAX;
985 for (k = 1; k <= size; k++)
986 { xassert(lit[k].col != NULL);
987 if (lit[k].neg == 0)
988 npp_add_aij(npp, row, lit[k].col, +1.0);
989 else if (lit[k].neg == 1)
990 { npp_add_aij(npp, row, lit[k].col, -1.0);
991 row->lb -= 1.0;
992 }
993 else
994 xassert(lit != lit);
995 }
996 return row;
997 }
998
999 /***********************************************************************
1000 * npp_sat_encode_geq - encode "not less than" constraint
1001 *
1002 * PURPOSE
1003 *
1004 * This routine translates to CNF the following constraint:
1005 *
1006 * n
1007 * sum 2**(k-1) * y[k] >= b, (1)
1008 * k=1
1009 *
1010 * where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
1011 * or constant false (zero), b is a given lower bound.
1012 *
1013 * ALGORITHM
1014 *
1015 * If b < 0, the constraint is redundant, so assume that b >= 0. Let
1016 *
1017 * n
1018 * b = sum 2**(k-1) b[k], (2)
1019 * k=1
1020 *
1021 * where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
1022 * therefore cannot be represented in the form (2), the constraint (1)
1023 * is infeasible.) In this case the condition (1) is equivalent to the
1024 * following condition:
1025 *
1026 * y[n] y[n-1] ... y[2] y[1] >= b[n] b[n-1] ... b[2] b[1], (3)
1027 *
1028 * where ">=" is understood lexicographically.
1029 *
1030 * Algorithmically the condition (3) can be tested as follows:
1031 *
1032 * for (k = n; k >= 1; k--)
1033 * { if (y[k] < b[k])
1034 * y is less than b;
1035 * if (y[k] > b[k])
1036 * y is greater than b;
1037 * }
1038 * y is equal to b;
1039 *
1040 * Thus, y is less than b iff there exists k, 1 <= k <= n, for which
1041 * the following condition is satisfied:
1042 *
1043 * y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] < b[k]. (4)
1044 *
1045 * Negating the condition (4) we have that y is not less than b iff for
1046 * all k, 1 <= k <= n, the following condition is satisfied:
1047 *
1048 * y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] >= b[k]. (5)
1049 *
1050 * Note that if b[k] = 0, the literal y[k] >= b[k] is always true, in
1051 * which case the entire clause (5) is true and can be omitted.
1052 *
1053 * RETURNS
1054 *
1055 * Normally the routine returns zero. However, if the constraint (1) is
1056 * infeasible, the routine returns non-zero. */
1057
1058 int npp_sat_encode_geq(NPP *npp, int n, NPPLIT y[], int rhs)
1059 { NPPLIT lit[1+NBIT_MAX];
1060 int j, k, size, temp, b[1+NBIT_MAX];
1061 xassert(0 <= n && n <= NBIT_MAX);
1062 /* if the constraint (1) is redundant, do nothing */
1063 if (rhs < 0)
1064 return 0;
1065 /* determine binary digits of b according to (2) */
1066 for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
1067 b[k] = temp & 1;
1068 if (temp != 0)
1069 { /* b >= 2**n; the constraint (1) is infeasible */
1070 return 1;
1071 }
1072 /* main transformation loop */
1073 for (k = 1; k <= n; k++)
1074 { /* build the clause (5) for current k */
1075 size = 0; /* clause size = number of literals */
1076 /* add literal y[k] >= b[k] */
1077 if (b[k] == 0)
1078 { /* b[k] = 0 -> the literal is true */
1079 goto skip;
1080 }
1081 else if (y[k].col == NULL)
1082 { /* y[k] = 0, b[k] = 1 -> the literal is false */
1083 xassert(y[k].neg == 0);
1084 }
1085 else
1086 { /* add literal y[k] = 1 */
1087 lit[++size] = y[k];
1088 }
1089 for (j = k+1; j <= n; j++)
1090 { /* add literal y[j] != b[j] */
1091 if (y[j].col == NULL)
1092 { xassert(y[j].neg == 0);
1093 if (b[j] == 0)
1094 { /* y[j] = 0, b[j] = 0 -> the literal is false */
1095 continue;
1096 }
1097 else
1098 { /* y[j] = 0, b[j] = 1 -> the literal is true */
1099 goto skip;
1100 }
1101 }
1102 else
1103 { lit[++size] = y[j];
1104 if (b[j] != 0)
1105 lit[size].neg = 1 - lit[size].neg;
1106 }
1107 }
1108 /* normalize the clause */
1109 size = npp_sat_normalize_clause(npp, size, lit);
1110 if (size < 0)
1111 { /* the clause is equivalent to the value true */
1112 goto skip;
1113 }
1114 if (size == 0)
1115 { /* the clause is equivalent to the value false; this means
1116 that the constraint (1) is infeasible */
1117 return 2;
1118 }
1119 /* translate the clause to corresponding cover inequality */
1120 npp_sat_encode_clause(npp, size, lit);
1121 skip: ;
1122 }
1123 return 0;
1124 }
1125
1126 /***********************************************************************
1127 * npp_sat_encode_leq - encode "not greater than" constraint
1128 *
1129 * PURPOSE
1130 *
1131 * This routine translates to CNF the following constraint:
1132 *
1133 * n
1134 * sum 2**(k-1) * y[k] <= b, (1)
1135 * k=1
1136 *
1137 * where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
1138 * or constant false (zero), b is a given upper bound.
1139 *
1140 * ALGORITHM
1141 *
1142 * If b < 0, the constraint is infeasible, so assume that b >= 0. Let
1143 *
1144 * n
1145 * b = sum 2**(k-1) b[k], (2)
1146 * k=1
1147 *
1148 * where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
1149 * therefore cannot be represented in the form (2), the constraint (1)
1150 * is redundant.) In this case the condition (1) is equivalent to the
1151 * following condition:
1152 *
1153 * y[n] y[n-1] ... y[2] y[1] <= b[n] b[n-1] ... b[2] b[1], (3)
1154 *
1155 * where "<=" is understood lexicographically.
1156 *
1157 * Algorithmically the condition (3) can be tested as follows:
1158 *
1159 * for (k = n; k >= 1; k--)
1160 * { if (y[k] < b[k])
1161 * y is less than b;
1162 * if (y[k] > b[k])
1163 * y is greater than b;
1164 * }
1165 * y is equal to b;
1166 *
1167 * Thus, y is greater than b iff there exists k, 1 <= k <= n, for which
1168 * the following condition is satisfied:
1169 *
1170 * y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] > b[k]. (4)
1171 *
1172 * Negating the condition (4) we have that y is not greater than b iff
1173 * for all k, 1 <= k <= n, the following condition is satisfied:
1174 *
1175 * y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] <= b[k]. (5)
1176 *
1177 * Note that if b[k] = 1, the literal y[k] <= b[k] is always true, in
1178 * which case the entire clause (5) is true and can be omitted.
1179 *
1180 * RETURNS
1181 *
1182 * Normally the routine returns zero. However, if the constraint (1) is
1183 * infeasible, the routine returns non-zero. */
1184
1185 int npp_sat_encode_leq(NPP *npp, int n, NPPLIT y[], int rhs)
1186 { NPPLIT lit[1+NBIT_MAX];
1187 int j, k, size, temp, b[1+NBIT_MAX];
1188 xassert(0 <= n && n <= NBIT_MAX);
1189 /* check if the constraint (1) is infeasible */
1190 if (rhs < 0)
1191 return 1;
1192 /* determine binary digits of b according to (2) */
1193 for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
1194 b[k] = temp & 1;
1195 if (temp != 0)
1196 { /* b >= 2**n; the constraint (1) is redundant */
1197 return 0;
1198 }
1199 /* main transformation loop */
1200 for (k = 1; k <= n; k++)
1201 { /* build the clause (5) for current k */
1202 size = 0; /* clause size = number of literals */
1203 /* add literal y[k] <= b[k] */
1204 if (b[k] == 1)
1205 { /* b[k] = 1 -> the literal is true */
1206 goto skip;
1207 }
1208 else if (y[k].col == NULL)
1209 { /* y[k] = 0, b[k] = 0 -> the literal is true */
1210 xassert(y[k].neg == 0);
1211 goto skip;
1212 }
1213 else
1214 { /* add literal y[k] = 0 */
1215 lit[++size] = y[k];
1216 lit[size].neg = 1 - lit[size].neg;
1217 }
1218 for (j = k+1; j <= n; j++)
1219 { /* add literal y[j] != b[j] */
1220 if (y[j].col == NULL)
1221 { xassert(y[j].neg == 0);
1222 if (b[j] == 0)
1223 { /* y[j] = 0, b[j] = 0 -> the literal is false */
1224 continue;
1225 }
1226 else
1227 { /* y[j] = 0, b[j] = 1 -> the literal is true */
1228 goto skip;
1229 }
1230 }
1231 else
1232 { lit[++size] = y[j];
1233 if (b[j] != 0)
1234 lit[size].neg = 1 - lit[size].neg;
1235 }
1236 }
1237 /* normalize the clause */
1238 size = npp_sat_normalize_clause(npp, size, lit);
1239 if (size < 0)
1240 { /* the clause is equivalent to the value true */
1241 goto skip;
1242 }
1243 if (size == 0)
1244 { /* the clause is equivalent to the value false; this means
1245 that the constraint (1) is infeasible */
1246 return 2;
1247 }
1248 /* translate the clause to corresponding cover inequality */
1249 npp_sat_encode_clause(npp, size, lit);
1250 skip: ;
1251 }
1252 return 0;
1253 }
1254
1255 /***********************************************************************
1256 * npp_sat_encode_row - encode constraint (row) of general type
1257 *
1258 * PURPOSE
1259 *
1260 * This routine translates to CNF the following constraint (row):
1261 *
1262 * L <= sum a[j] x[j] <= U, (1)
1263 * j
1264 *
1265 * where all x[j] are binary variables.
1266 *
1267 * ALGORITHM
1268 *
1269 * First, the routine performs substitution x[j] = t[j] for j in J+
1270 * and x[j] = 1 - t[j] for j in J-, where J+ = { j : a[j] > 0 } and
1271 * J- = { j : a[j] < 0 }. This gives:
1272 *
1273 * L <= sum a[j] t[j] + sum a[j] (1 - t[j]) <= U ==>
1274 * j in J+ j in J-
1275 *
1276 * L' <= sum |a[j]| t[j] <= U', (2)
1277 * j
1278 *
1279 * where
1280 *
1281 * L' = L - sum a[j], U' = U - sum a[j]. (3)
1282 * j in J- j in J-
1283 *
1284 * (Actually only new bounds L' and U' are computed.)
1285 *
1286 * Then the routine translates to CNF the following equality:
1287 *
1288 * n
1289 * sum |a[j]| t[j] = sum 2**(k-1) * y[k], (4)
1290 * j k=1
1291 *
1292 * where y[k] is either some t[j] or a new literal or a constant zero
1293 * (see the routine npp_sat_encode_sum_ax).
1294 *
1295 * Finally, the routine translates to CNF the following conditions:
1296 *
1297 * n
1298 * sum 2**(k-1) * y[k] >= L' (5)
1299 * k=1
1300 *
1301 * and
1302 *
1303 * n
1304 * sum 2**(k-1) * y[k] <= U' (6)
1305 * k=1
1306 *
1307 * (see the routines npp_sat_encode_geq and npp_sat_encode_leq).
1308 *
1309 * All resulting clauses are encoded as cover inequalities and included
1310 * into the transformed problem.
1311 *
1312 * Note that on exit the routine removes the specified constraint (row)
1313 * from the original problem.
1314 *
1315 * RETURNS
1316 *
1317 * The routine returns one of the following codes:
1318 *
1319 * 0 - translation was successful;
1320 * 1 - constraint (1) was found infeasible;
1321 * 2 - integer arithmetic error occured. */
1322
1323 int npp_sat_encode_row(NPP *npp, NPPROW *row)
1324 { NPPAIJ *aij;
1325 NPPLIT y[1+NBIT_MAX];
1326 int n, rhs;
1327 double lb, ub;
1328 /* the row should not be free */
1329 xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX));
1330 /* compute new bounds L' and U' (3) */
1331 lb = row->lb;
1332 ub = row->ub;
1333 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1334 { if (aij->val < 0.0)
1335 { if (lb != -DBL_MAX)
1336 lb -= aij->val;
1337 if (ub != -DBL_MAX)
1338 ub -= aij->val;
1339 }
1340 }
1341 /* encode the equality (4) */
1342 n = npp_sat_encode_sum_ax(npp, row, y);
1343 if (n < 0)
1344 return 2; /* integer arithmetic error */
1345 /* encode the condition (5) */
1346 if (lb != -DBL_MAX)
1347 { rhs = (int)lb;
1348 if ((double)rhs != lb)
1349 return 2; /* integer arithmetic error */
1350 if (npp_sat_encode_geq(npp, n, y, rhs) != 0)
1351 return 1; /* original constraint is infeasible */
1352 }
1353 /* encode the condition (6) */
1354 if (ub != +DBL_MAX)
1355 { rhs = (int)ub;
1356 if ((double)rhs != ub)
1357 return 2; /* integer arithmetic error */
1358 if (npp_sat_encode_leq(npp, n, y, rhs) != 0)
1359 return 1; /* original constraint is infeasible */
1360 }
1361 /* remove the specified row from the problem */
1362 npp_del_row(npp, row);
1363 return 0;
1364 }
1365
1366 /***********************************************************************
1367 * npp_sat_encode_prob - encode 0-1 feasibility problem
1368 *
1369 * This routine translates the specified 0-1 feasibility problem to an
1370 * equivalent SAT-CNF problem.
1371 *
1372 * N.B. Currently this is a very crude implementation.
1373 *
1374 * RETURNS
1375 *
1376 * 0 success;
1377 *
1378 * GLP_ENOPFS primal/integer infeasibility detected;
1379 *
1380 * GLP_ERANGE integer overflow occured. */
1381
1382 int npp_sat_encode_prob(NPP *npp)
1383 { NPPROW *row, *next_row, *prev_row;
1384 NPPCOL *col, *next_col;
1385 int cover = 0, pack = 0, partn = 0, ret;
1386 /* process and remove free rows */
1387 for (row = npp->r_head; row != NULL; row = next_row)
1388 { next_row = row->next;
1389 if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
1390 npp_sat_free_row(npp, row);
1391 }
1392 /* process and remove fixed columns */
1393 for (col = npp->c_head; col != NULL; col = next_col)
1394 { next_col = col->next;
1395 if (col->lb == col->ub)
1396 xassert(npp_sat_fixed_col(npp, col) == 0);
1397 }
1398 /* only binary variables should remain */
1399 for (col = npp->c_head; col != NULL; col = col->next)
1400 xassert(col->is_int && col->lb == 0.0 && col->ub == 1.0);
1401 /* new rows may be added to the end of the row list, so we walk
1402 from the end to beginning of the list */
1403 for (row = npp->r_tail; row != NULL; row = prev_row)
1404 { prev_row = row->prev;
1405 /* process special cases */
1406 ret = npp_sat_is_cover_ineq(npp, row);
1407 if (ret != 0)
1408 { /* row is covering inequality */
1409 cover++;
1410 /* since it already encodes a clause, just transform it to
1411 canonical form */
1412 if (ret == 2)
1413 { xassert(npp_sat_reverse_row(npp, row) == 0);
1414 ret = npp_sat_is_cover_ineq(npp, row);
1415 }
1416 xassert(ret == 1);
1417 continue;
1418 }
1419 ret = npp_sat_is_partn_eq(npp, row);
1420 if (ret != 0)
1421 { /* row is partitioning equality */
1422 NPPROW *cov;
1423 NPPAIJ *aij;
1424 partn++;
1425 /* transform it to canonical form */
1426 if (ret == 2)
1427 { xassert(npp_sat_reverse_row(npp, row) == 0);
1428 ret = npp_sat_is_partn_eq(npp, row);
1429 }
1430 xassert(ret == 1);
1431 /* and split it into covering and packing inequalities,
1432 both in canonical forms */
1433 cov = npp_add_row(npp);
1434 cov->lb = row->lb, cov->ub = +DBL_MAX;
1435 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1436 npp_add_aij(npp, cov, aij->col, aij->val);
1437 xassert(npp_sat_is_cover_ineq(npp, cov) == 1);
1438 /* the cover inequality already encodes a clause and do
1439 not need any further processing */
1440 row->lb = -DBL_MAX;
1441 xassert(npp_sat_is_pack_ineq(npp, row) == 1);
1442 /* the packing inequality will be processed below */
1443 pack--;
1444 }
1445 ret = npp_sat_is_pack_ineq(npp, row);
1446 if (ret != 0)
1447 { /* row is packing inequality */
1448 NPPROW *rrr;
1449 int nlit, desired_nlit = 4;
1450 pack++;
1451 /* transform it to canonical form */
1452 if (ret == 2)
1453 { xassert(npp_sat_reverse_row(npp, row) == 0);
1454 ret = npp_sat_is_pack_ineq(npp, row);
1455 }
1456 xassert(ret == 1);
1457 /* process the packing inequality */
1458 for (;;)
1459 { /* determine the number of literals in the remaining
1460 inequality */
1461 nlit = npp_row_nnz(npp, row);
1462 if (nlit <= desired_nlit)
1463 break;
1464 /* split the current inequality into one having not more
1465 than desired_nlit literals and remaining one */
1466 rrr = npp_sat_split_pack(npp, row, desired_nlit-1);
1467 /* translate the former inequality to CNF and remove it
1468 from the original problem */
1469 npp_sat_encode_pack(npp, rrr);
1470 }
1471 /* translate the remaining inequality to CNF and remove it
1472 from the original problem */
1473 npp_sat_encode_pack(npp, row);
1474 continue;
1475 }
1476 /* translate row of general type to CNF and remove it from the
1477 original problem */
1478 ret = npp_sat_encode_row(npp, row);
1479 if (ret == 0)
1480 ;
1481 else if (ret == 1)
1482 ret = GLP_ENOPFS;
1483 else if (ret == 2)
1484 ret = GLP_ERANGE;
1485 else
1486 xassert(ret != ret);
1487 if (ret != 0)
1488 goto done;
1489 }
1490 ret = 0;
1491 if (cover != 0)
1492 xprintf("%d covering inequalities\n", cover);
1493 if (pack != 0)
1494 xprintf("%d packing inequalities\n", pack);
1495 if (partn != 0)
1496 xprintf("%d partitioning equalities\n", partn);
1497 done: return ret;
1498 }
1499
1500 /* eof */