COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnpp06.c @ 9:33de93886c88

subpack-glpk
Last change on this file since 9:33de93886c88 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 12 years ago

Import GLPK 4.47

File size: 50.3 KB
Line 
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
38void 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
83struct 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
91static int rcv_sat_fixed_col(NPP *, void *);
92
93int 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      }
124skip: /* remove the column from the problem */
125      npp_del_col(npp, q);
126      return 0;
127}
128
129static 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
143int 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
163int 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
180int 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
224int 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
283int 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
342int 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
377int 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
435NPPROW *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
494void 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
541void 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
640void 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
822static 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
836int 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
931int 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];
954skip:    ;
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
979NPPROW *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
1058int 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);
1121skip:    ;
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
1185int 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);
1250skip:    ;
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
1323int 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
1382int 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);
1497done: return ret;
1498}
1499
1500/* eof */
Note: See TracBrowser for help on using the repository browser.