COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpnpp04.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 49.3 KB
Line 
1/* glpnpp04.c */
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 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*  NAME
29*
30*  npp_binarize_prob - binarize MIP problem
31*
32*  SYNOPSIS
33*
34*  #include "glpnpp.h"
35*  int npp_binarize_prob(NPP *npp);
36*
37*  DESCRIPTION
38*
39*  The routine npp_binarize_prob replaces in the original MIP problem
40*  every integer variable:
41*
42*     l[q] <= x[q] <= u[q],                                          (1)
43*
44*  where l[q] < u[q], by an equivalent sum of binary variables.
45*
46*  RETURNS
47*
48*  The routine returns the number of integer variables for which the
49*  transformation failed, because u[q] - l[q] > d_max.
50*
51*  PROBLEM TRANSFORMATION
52*
53*  If variable x[q] has non-zero lower bound, it is first processed
54*  with the routine npp_lbnd_col. Thus, we can assume that:
55*
56*     0 <= x[q] <= u[q].                                             (2)
57*
58*  If u[q] = 1, variable x[q] is already binary, so further processing
59*  is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a
60*  smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2).
61*  Then variable x[q] can be replaced by the following sum:
62*
63*            n-1
64*     x[q] = sum 2^k x[k],                                           (3)
65*            k=0
66*
67*  where x[k] are binary columns (variables). If u[q] < 2^n - 1, the
68*  following additional inequality constraint must be also included in
69*  the transformed problem:
70*
71*     n-1
72*     sum 2^k x[k] <= u[q].                                          (4)
73*     k=0
74*
75*  Note: Assuming that in the transformed problem x[q] becomes binary
76*  variable x[0], this transformation causes new n-1 binary variables
77*  to appear.
78*
79*  Substituting x[q] from (3) to the objective row gives:
80*
81*     z = sum c[j] x[j] + c[0] =
82*          j
83*
84*       = sum c[j] x[j] + c[q] x[q] + c[0] =
85*         j!=q
86*                              n-1
87*       = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] =
88*         j!=q                 k=0
89*                         n-1
90*       = sum c[j] x[j] + sum c[k] x[k] + c[0],
91*         j!=q            k=0
92*
93*  where:
94*
95*     c[k] = 2^k c[q],  k = 0, ..., n-1.                             (5)
96*
97*  And substituting x[q] from (3) to i-th constraint row i gives:
98*
99*     L[i] <= sum a[i,j] x[j] <= U[i]  ==>
100*              j
101*
102*     L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i]  ==>
103*             j!=q
104*                                      n-1
105*     L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i]  ==>
106*             j!=q                     k=0
107*                               n-1
108*     L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i],
109*             j!=q              k=0
110*
111*  where:
112*
113*     a[i,k] = 2^k a[i,q],  k = 0, ..., n-1.                         (6)
114*
115*  RECOVERING SOLUTION
116*
117*  Value of variable x[q] is computed with formula (3). */
118
119struct binarize
120{     int q;
121      /* column reference number for x[q] = x[0] */
122      int j;
123      /* column reference number for x[1]; x[2] has reference number
124         j+1, x[3] - j+2, etc. */
125      int n;
126      /* total number of binary variables, n >= 2 */
127};
128
129static int rcv_binarize_prob(NPP *npp, void *info);
130
131int npp_binarize_prob(NPP *npp)
132{     /* binarize MIP problem */
133      struct binarize *info;
134      NPPROW *row;
135      NPPCOL *col, *bin;
136      NPPAIJ *aij;
137      int u, n, k, temp, nfails, nvars, nbins, nrows;
138      /* new variables will be added to the end of the column list, so
139         we go from the end to beginning of the column list */
140      nfails = nvars = nbins = nrows = 0;
141      for (col = npp->c_tail; col != NULL; col = col->prev)
142      {  /* skip continuous variable */
143         if (!col->is_int) continue;
144         /* skip fixed variable */
145         if (col->lb == col->ub) continue;
146         /* skip binary variable */
147         if (col->lb == 0.0 && col->ub == 1.0) continue;
148         /* check if the transformation is applicable */
149         if (col->lb < -1e6 || col->ub > +1e6 ||
150             col->ub - col->lb > 4095.0)
151         {  /* unfortunately, not */
152            nfails++;
153            continue;
154         }
155         /* process integer non-binary variable x[q] */
156         nvars++;
157         /* make x[q] non-negative, if its lower bound is non-zero */
158         if (col->lb != 0.0)
159            npp_lbnd_col(npp, col);
160         /* now 0 <= x[q] <= u[q] */
161         xassert(col->lb == 0.0);
162         u = (int)col->ub;
163         xassert(col->ub == (double)u);
164         /* if x[q] is binary, further processing is not needed */
165         if (u == 1) continue;
166         /* determine smallest n such that u <= 2^n - 1 (thus, n is the
167            number of binary variables needed) */
168         n = 2, temp = 4;
169         while (u >= temp)
170            n++, temp += temp;
171         nbins += n;
172         /* create transformation stack entry */
173         info = npp_push_tse(npp,
174            rcv_binarize_prob, sizeof(struct binarize));
175         info->q = col->j;
176         info->j = 0; /* will be set below */
177         info->n = n;
178         /* if u < 2^n - 1, we need one additional row for (4) */
179         if (u < temp - 1)
180         {  row = npp_add_row(npp), nrows++;
181            row->lb = -DBL_MAX, row->ub = u;
182         }
183         else
184            row = NULL;
185         /* in the transformed problem variable x[q] becomes binary
186            variable x[0], so its objective and constraint coefficients
187            are not changed */
188         col->ub = 1.0;
189         /* include x[0] into constraint (4) */
190         if (row != NULL)
191            npp_add_aij(npp, row, col, 1.0);
192         /* add other binary variables x[1], ..., x[n-1] */
193         for (k = 1, temp = 2; k < n; k++, temp += temp)
194         {  /* add new binary variable x[k] */
195            bin = npp_add_col(npp);
196            bin->is_int = 1;
197            bin->lb = 0.0, bin->ub = 1.0;
198            bin->coef = (double)temp * col->coef;
199            /* store column reference number for x[1] */
200            if (info->j == 0)
201               info->j = bin->j;
202            else
203               xassert(info->j + (k-1) == bin->j);
204            /* duplicate constraint coefficients for x[k]; this also
205               automatically includes x[k] into constraint (4) */
206            for (aij = col->ptr; aij != NULL; aij = aij->c_next)
207               npp_add_aij(npp, aij->row, bin, (double)temp * aij->val);
208         }
209      }
210      if (nvars > 0)
211         xprintf("%d integer variable(s) were replaced by %d binary one"
212            "s\n", nvars, nbins);
213      if (nrows > 0)
214         xprintf("%d row(s) were added due to binarization\n", nrows);
215      if (nfails > 0)
216         xprintf("Binarization failed for %d integer variable(s)\n",
217            nfails);
218      return nfails;
219}
220
221static int rcv_binarize_prob(NPP *npp, void *_info)
222{     /* recovery binarized variable */
223      struct binarize *info = _info;
224      int k, temp;
225      double sum;
226      /* compute value of x[q]; see formula (3) */
227      sum = npp->c_value[info->q];
228      for (k = 1, temp = 2; k < info->n; k++, temp += temp)
229         sum += (double)temp * npp->c_value[info->j + (k-1)];
230      npp->c_value[info->q] = sum;
231      return 0;
232}
233
234/**********************************************************************/
235
236struct elem
237{     /* linear form element a[j] x[j] */
238      double aj;
239      /* non-zero coefficient value */
240      NPPCOL *xj;
241      /* pointer to variable (column) */
242      struct elem *next;
243      /* pointer to another term */
244};
245
246static struct elem *copy_form(NPP *npp, NPPROW *row, double s)
247{     /* copy linear form */
248      NPPAIJ *aij;
249      struct elem *ptr, *e;
250      ptr = NULL;
251      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
252      {  e = dmp_get_atom(npp->pool, sizeof(struct elem));
253         e->aj = s * aij->val;
254         e->xj = aij->col;
255         e->next = ptr;
256         ptr = e;
257      }
258      return ptr;
259}
260
261static void drop_form(NPP *npp, struct elem *ptr)
262{     /* drop linear form */
263      struct elem *e;
264      while (ptr != NULL)
265      {  e = ptr;
266         ptr = e->next;
267         dmp_free_atom(npp->pool, e, sizeof(struct elem));
268      }
269      return;
270}
271
272/***********************************************************************
273*  NAME
274*
275*  npp_is_packing - test if constraint is packing inequality
276*
277*  SYNOPSIS
278*
279*  #include "glpnpp.h"
280*  int npp_is_packing(NPP *npp, NPPROW *row);
281*
282*  RETURNS
283*
284*  If the specified row (constraint) is packing inequality (see below),
285*  the routine npp_is_packing returns non-zero. Otherwise, it returns
286*  zero.
287*
288*  PACKING INEQUALITIES
289*
290*  In canonical format the packing inequality is the following:
291*
292*     sum  x[j] <= 1,                                                (1)
293*    j in J
294*
295*  where all variables x[j] are binary. This inequality expresses the
296*  condition that in any integer feasible solution at most one variable
297*  from set J can take non-zero (unity) value while other variables
298*  must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because
299*  if J is empty or |J| = 1, the inequality (1) is redundant.
300*
301*  In general case the packing inequality may include original variables
302*  x[j] as well as their complements x~[j]:
303*
304*     sum   x[j] + sum   x~[j] <= 1,                                 (2)
305*    j in Jp      j in Jn
306*
307*  where Jp and Jn are not intersected. Therefore, using substitution
308*  x~[j] = 1 - x[j] gives the packing inequality in generalized format:
309*
310*     sum   x[j] - sum   x[j] <= 1 - |Jn|.                           (3)
311*    j in Jp      j in Jn */
312
313int npp_is_packing(NPP *npp, NPPROW *row)
314{     /* test if constraint is packing inequality */
315      NPPCOL *col;
316      NPPAIJ *aij;
317      int b;
318      xassert(npp == npp);
319      if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX))
320         return 0;
321      b = 1;
322      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
323      {  col = aij->col;
324         if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
325            return 0;
326         if (aij->val == +1.0)
327            ;
328         else if (aij->val == -1.0)
329            b--;
330         else
331            return 0;
332      }
333      if (row->ub != (double)b) return 0;
334      return 1;
335}
336
337/***********************************************************************
338*  NAME
339*
340*  npp_hidden_packing - identify hidden packing inequality
341*
342*  SYNOPSIS
343*
344*  #include "glpnpp.h"
345*  int npp_hidden_packing(NPP *npp, NPPROW *row);
346*
347*  DESCRIPTION
348*
349*  The routine npp_hidden_packing processes specified inequality
350*  constraint, which includes only binary variables, and the number of
351*  the variables is not less than two. If the original inequality is
352*  equivalent to a packing inequality, the routine replaces it by this
353*  equivalent inequality. If the original constraint is double-sided
354*  inequality, it is replaced by a pair of single-sided inequalities,
355*  if necessary.
356*
357*  RETURNS
358*
359*  If the original inequality constraint was replaced by equivalent
360*  packing inequality, the routine npp_hidden_packing returns non-zero.
361*  Otherwise, it returns zero.
362*
363*  PROBLEM TRANSFORMATION
364*
365*  Consider an inequality constraint:
366*
367*     sum  a[j] x[j] <= b,                                           (1)
368*    j in J
369*
370*  where all variables x[j] are binary, and |J| >= 2. (In case of '>='
371*  inequality it can be transformed to '<=' format by multiplying both
372*  its sides by -1.)
373*
374*  Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
375*  x[j] = 1 - x~[j] for all j in Jn, we have:
376*
377*     sum   a[j] x[j] <= b  ==>
378*    j in J
379*
380*     sum   a[j] x[j] + sum   a[j] x[j] <= b  ==>
381*    j in Jp           j in Jn
382*
383*     sum   a[j] x[j] + sum   a[j] (1 - x~[j]) <= b  ==>
384*    j in Jp           j in Jn
385*
386*     sum   a[j] x[j] - sum   a[j] x~[j] <= b - sum   a[j].
387*    j in Jp           j in Jn                 j in Jn
388*
389*  Thus, meaning the transformation above, we can assume that in
390*  inequality (1) all coefficients a[j] are positive. Moreover, we can
391*  assume that a[j] <= b. In fact, let a[j] > b; then the following
392*  three cases are possible:
393*
394*  1) b < 0. In this case inequality (1) is infeasible, so the problem
395*     has no feasible solution (see the routine npp_analyze_row);
396*
397*  2) b = 0. In this case inequality (1) is a forcing inequality on its
398*     upper bound (see the routine npp_forcing row), from which it
399*     follows that all variables x[j] should be fixed at zero;
400*
401*  3) b > 0. In this case inequality (1) defines an implied zero upper
402*     bound for variable x[j] (see the routine npp_implied_bounds), from
403*     which it follows that x[j] should be fixed at zero.
404*
405*  It is assumed that all three cases listed above have been recognized
406*  by the routine npp_process_prob, which performs basic MIP processing
407*  prior to a call the routine npp_hidden_packing. So, if one of these
408*  cases occurs, we should just skip processing such constraint.
409*
410*  Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is
411*  equivalent to packing inquality only if:
412*
413*     a[j] + a[k] > b + eps                                          (2)
414*
415*  for all j, k in J, j != k, where eps is an absolute tolerance for
416*  row (linear form) value. Checking the condition (2) for all j and k,
417*  j != k, requires time O(|J|^2). However, this time can be reduced to
418*  O(|J|), if use minimal a[j] and a[k], in which case it is sufficient
419*  to check the condition (2) only once.
420*
421*  Once the original inequality (1) is replaced by equivalent packing
422*  inequality, we need to perform back substitution x~[j] = 1 - x[j] for
423*  all j in Jn (see above).
424*
425*  RECOVERING SOLUTION
426*
427*  None needed. */
428
429static int hidden_packing(NPP *npp, struct elem *ptr, double *_b)
430{     /* process inequality constraint: sum a[j] x[j] <= b;
431         0 - specified row is NOT hidden packing inequality;
432         1 - specified row is packing inequality;
433         2 - specified row is hidden packing inequality. */
434      struct elem *e, *ej, *ek;
435      int neg;
436      double b = *_b, eps;
437      xassert(npp == npp);
438      /* a[j] must be non-zero, x[j] must be binary, for all j in J */
439      for (e = ptr; e != NULL; e = e->next)
440      {  xassert(e->aj != 0.0);
441         xassert(e->xj->is_int);
442         xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
443      }
444      /* check if the specified inequality constraint already has the
445         form of packing inequality */
446      neg = 0; /* neg is |Jn| */
447      for (e = ptr; e != NULL; e = e->next)
448      {  if (e->aj == +1.0)
449            ;
450         else if (e->aj == -1.0)
451            neg++;
452         else
453            break;
454      }
455      if (e == NULL)
456      {  /* all coefficients a[j] are +1 or -1; check rhs b */
457         if (b == (double)(1 - neg))
458         {  /* it is packing inequality; no processing is needed */
459            return 1;
460         }
461      }
462      /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
463         positive; the result is a~[j] = |a[j]| and new rhs b */
464      for (e = ptr; e != NULL; e = e->next)
465         if (e->aj < 0) b -= e->aj;
466      /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
467      /* if a[j] > b, skip processing--this case must not appear */
468      for (e = ptr; e != NULL; e = e->next)
469         if (fabs(e->aj) > b) return 0;
470      /* now 0 < a[j] <= b for all j in J */
471      /* find two minimal coefficients a[j] and a[k], j != k */
472      ej = NULL;
473      for (e = ptr; e != NULL; e = e->next)
474         if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e;
475      xassert(ej != NULL);
476      ek = NULL;
477      for (e = ptr; e != NULL; e = e->next)
478         if (e != ej)
479            if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e;
480      xassert(ek != NULL);
481      /* the specified constraint is equivalent to packing inequality
482         iff a[j] + a[k] > b + eps */
483      eps = 1e-3 + 1e-6 * fabs(b);
484      if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0;
485      /* perform back substitution x~[j] = 1 - x[j] and construct the
486         final equivalent packing inequality in generalized format */
487      b = 1.0;
488      for (e = ptr; e != NULL; e = e->next)
489      {  if (e->aj > 0.0)
490            e->aj = +1.0;
491         else /* e->aj < 0.0 */
492            e->aj = -1.0, b -= 1.0;
493      }
494      *_b = b;
495      return 2;
496}
497
498int npp_hidden_packing(NPP *npp, NPPROW *row)
499{     /* identify hidden packing inequality */
500      NPPROW *copy;
501      NPPAIJ *aij;
502      struct elem *ptr, *e;
503      int kase, ret, count = 0;
504      double b;
505      /* the row must be inequality constraint */
506      xassert(row->lb < row->ub);
507      for (kase = 0; kase <= 1; kase++)
508      {  if (kase == 0)
509         {  /* process row upper bound */
510            if (row->ub == +DBL_MAX) continue;
511            ptr = copy_form(npp, row, +1.0);
512            b = + row->ub;
513         }
514         else
515         {  /* process row lower bound */
516            if (row->lb == -DBL_MAX) continue;
517            ptr = copy_form(npp, row, -1.0);
518            b = - row->lb;
519         }
520         /* now the inequality has the form "sum a[j] x[j] <= b" */
521         ret = hidden_packing(npp, ptr, &b);
522         xassert(0 <= ret && ret <= 2);
523         if (kase == 1 && ret == 1 || ret == 2)
524         {  /* the original inequality has been identified as hidden
525               packing inequality */
526            count++;
527#ifdef GLP_DEBUG
528            xprintf("Original constraint:\n");
529            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
530               xprintf(" %+g x%d", aij->val, aij->col->j);
531            if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
532            if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
533            xprintf("\n");
534            xprintf("Equivalent packing inequality:\n");
535            for (e = ptr; e != NULL; e = e->next)
536               xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
537            xprintf(", <= %g\n", b);
538#endif
539            if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
540            {  /* the original row is single-sided inequality; no copy
541                  is needed */
542               copy = NULL;
543            }
544            else
545            {  /* the original row is double-sided inequality; we need
546                  to create its copy for other bound before replacing it
547                  with the equivalent inequality */
548               copy = npp_add_row(npp);
549               if (kase == 0)
550               {  /* the copy is for lower bound */
551                  copy->lb = row->lb, copy->ub = +DBL_MAX;
552               }
553               else
554               {  /* the copy is for upper bound */
555                  copy->lb = -DBL_MAX, copy->ub = row->ub;
556               }
557               /* copy original row coefficients */
558               for (aij = row->ptr; aij != NULL; aij = aij->r_next)
559                  npp_add_aij(npp, copy, aij->col, aij->val);
560            }
561            /* replace the original inequality by equivalent one */
562            npp_erase_row(npp, row);
563            row->lb = -DBL_MAX, row->ub = b;
564            for (e = ptr; e != NULL; e = e->next)
565               npp_add_aij(npp, row, e->xj, e->aj);
566            /* continue processing lower bound for the copy */
567            if (copy != NULL) row = copy;
568         }
569         drop_form(npp, ptr);
570      }
571      return count;
572}
573
574/***********************************************************************
575*  NAME
576*
577*  npp_implied_packing - identify implied packing inequality
578*
579*  SYNOPSIS
580*
581*  #include "glpnpp.h"
582*  int npp_implied_packing(NPP *npp, NPPROW *row, int which,
583*     NPPCOL *var[], char set[]);
584*
585*  DESCRIPTION
586*
587*  The routine npp_implied_packing processes specified row (constraint)
588*  of general format:
589*
590*     L <= sum a[j] x[j] <= U.                                       (1)
591*           j
592*
593*  If which = 0, only lower bound L, which must exist, is considered,
594*  while upper bound U is ignored. Similarly, if which = 1, only upper
595*  bound U, which must exist, is considered, while lower bound L is
596*  ignored. Thus, if the specified row is a double-sided inequality or
597*  equality constraint, this routine should be called twice for both
598*  lower and upper bounds.
599*
600*  The routine npp_implied_packing attempts to find a non-trivial (i.e.
601*  having not less than two binary variables) packing inequality:
602*
603*     sum   x[j] - sum   x[j] <= 1 - |Jn|,                           (2)
604*    j in Jp      j in Jn
605*
606*  which is relaxation of the constraint (1) in the sense that any
607*  solution satisfying to that constraint also satisfies to the packing
608*  inequality (2). If such relaxation exists, the routine stores
609*  pointers to descriptors of corresponding binary variables and their
610*  flags, resp., to locations var[1], var[2], ..., var[len] and set[1],
611*  set[2], ..., set[len], where set[j] = 0 means that j in Jp and
612*  set[j] = 1 means that j in Jn.
613*
614*  RETURNS
615*
616*  The routine npp_implied_packing returns len, which is the total
617*  number of binary variables in the packing inequality found, len >= 2.
618*  However, if the relaxation does not exist, the routine returns zero.
619*
620*  ALGORITHM
621*
622*  If which = 0, the constraint coefficients (1) are multiplied by -1
623*  and b is assigned -L; if which = 1, the constraint coefficients (1)
624*  are not changed and b is assigned +U. In both cases the specified
625*  constraint gets the following format:
626*
627*     sum a[j] x[j] <= b.                                            (3)
628*      j
629*
630*  (Note that (3) is a relaxation of (1), because one of bounds L or U
631*  is ignored.)
632*
633*  Let J be set of binary variables, Kp be set of non-binary (integer
634*  or continuous) variables with a[j] > 0, and Kn be set of non-binary
635*  variables with a[j] < 0. Then the inequality (3) can be written as
636*  follows:
637*
638*     sum  a[j] x[j] <= b - sum   a[j] x[j] - sum   a[j] x[j].       (4)
639*    j in J                j in Kp           j in Kn
640*
641*  To get rid of non-binary variables we can replace the inequality (4)
642*  by the following relaxed inequality:
643*
644*     sum  a[j] x[j] <= b~,                                          (5)
645*    j in J
646*
647*  where:
648*
649*     b~ = sup(b - sum   a[j] x[j] - sum   a[j] x[j]) =
650*                 j in Kp           j in Kn
651*
652*        = b - inf sum   a[j] x[j] - inf sum   a[j] x[j] =           (6)
653*                 j in Kp               j in Kn
654*
655*        = b - sum   a[j] l[j] - sum   a[j] u[j].
656*             j in Kp           j in Kn
657*
658*  Note that if lower bound l[j] (if j in Kp) or upper bound u[j]
659*  (if j in Kn) of some non-binary variable x[j] does not exist, then
660*  formally b = +oo, in which case further analysis is not performed.
661*
662*  Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all
663*  the inequality coefficients in (5) positive, we replace all x[j] in
664*  Bn by their complementaries, substituting x[j] = 1 - x~[j] for all
665*  j in Bn, that gives:
666*
667*     sum   a[j] x[j] - sum   a[j] x~[j] <= b~ - sum   a[j].         (7)
668*    j in Bp           j in Bn                  j in Bn
669*
670*  This inequality is a relaxation of the original constraint (1), and
671*  it is a binary knapsack inequality. Writing it in the standard format
672*  we have:
673*
674*     sum  alfa[j] z[j] <= beta,                                     (8)
675*    j in J
676*
677*  where:
678*               ( + a[j],   if j in Bp,
679*     alfa[j] = <                                                    (9)
680*               ( - a[j],   if j in Bn,
681*
682*               ( x[j],     if j in Bp,
683*        z[j] = <                                                   (10)
684*               ( 1 - x[j], if j in Bn,
685*
686*        beta = b~ - sum   a[j].                                    (11)
687*                   j in Bn
688*
689*  In the inequality (8) all coefficients are positive, therefore, the
690*  packing relaxation to be found for this inequality is the following:
691*
692*     sum  z[j] <= 1.                                               (12)
693*    j in P
694*
695*  It is obvious that set P within J, which we would like to find, must
696*  satisfy to the following condition:
697*
698*     alfa[j] + alfa[k] > beta + eps  for all j, k in P, j != k,    (13)
699*
700*  where eps is an absolute tolerance for value of the linear form.
701*  Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}.
702*  Moreover, if in the equality (8) there exist coefficients alfa[k],
703*  for which alfa[k] <= (beta + eps) / 2, but which, nevertheless,
704*  satisfies to the condition (13) for all j in P, *one* corresponding
705*  variable z[k] (having, for example, maximal coefficient alfa[k]) can
706*  be included in set P, that allows increasing the number of binary
707*  variables in (12) by one.
708*
709*  Once the set P has been built, for the inequality (12) we need to
710*  perform back substitution according to (10) in order to express it
711*  through the original binary variables. As the result of such back
712*  substitution the relaxed packing inequality get its final format (2),
713*  where Jp = J intersect Bp, and Jn = J intersect Bn. */
714
715int npp_implied_packing(NPP *npp, NPPROW *row, int which,
716      NPPCOL *var[], char set[])
717{     struct elem *ptr, *e, *i, *k;
718      int len = 0;
719      double b, eps;
720      /* build inequality (3) */
721      if (which == 0)
722      {  ptr = copy_form(npp, row, -1.0);
723         xassert(row->lb != -DBL_MAX);
724         b = - row->lb;
725      }
726      else if (which == 1)
727      {  ptr = copy_form(npp, row, +1.0);
728         xassert(row->ub != +DBL_MAX);
729         b = + row->ub;
730      }
731      /* remove non-binary variables to build relaxed inequality (5);
732         compute its right-hand side b~ with formula (6) */
733      for (e = ptr; e != NULL; e = e->next)
734      {  if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
735         {  /* x[j] is non-binary variable */
736            if (e->aj > 0.0)
737            {  if (e->xj->lb == -DBL_MAX) goto done;
738               b -= e->aj * e->xj->lb;
739            }
740            else /* e->aj < 0.0 */
741            {  if (e->xj->ub == +DBL_MAX) goto done;
742               b -= e->aj * e->xj->ub;
743            }
744            /* a[j] = 0 means that variable x[j] is removed */
745            e->aj = 0.0;
746         }
747      }
748      /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8);
749         compute its right-hand side beta with formula (11) */
750      for (e = ptr; e != NULL; e = e->next)
751         if (e->aj < 0.0) b -= e->aj;
752      /* if beta is close to zero, the knapsack inequality is either
753         infeasible or forcing inequality; this must never happen, so
754         we skip further analysis */
755      if (b < 1e-3) goto done;
756      /* build set P as well as sets Jp and Jn, and determine x[k] as
757         explained above in comments to the routine */
758      eps = 1e-3 + 1e-6 * b;
759      i = k = NULL;
760      for (e = ptr; e != NULL; e = e->next)
761      {  /* note that alfa[j] = |a[j]| */
762         if (fabs(e->aj) > 0.5 * (b + eps))
763         {  /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in
764               set Jp or Jn */
765            var[++len] = e->xj;
766            set[len] = (char)(e->aj > 0.0 ? 0 : 1);
767            /* alfa[i] = min alfa[j] over all j included in set P */
768            if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e;
769         }
770         else if (fabs(e->aj) >= 1e-3)
771         {  /* alfa[k] = max alfa[j] over all j not included in set P;
772               we skip coefficient a[j] if it is close to zero to avoid
773               numerically unreliable results */
774            if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e;
775         }
776      }
777      /* if alfa[k] satisfies to condition (13) for all j in P, include
778         x[k] in P */
779      if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps)
780      {  var[++len] = k->xj;
781         set[len] = (char)(k->aj > 0.0 ? 0 : 1);
782      }
783      /* trivial packing inequality being redundant must never appear,
784         so we just ignore it */
785      if (len < 2) len = 0;
786done: drop_form(npp, ptr);
787      return len;
788}
789
790/***********************************************************************
791*  NAME
792*
793*  npp_is_covering - test if constraint is covering inequality
794*
795*  SYNOPSIS
796*
797*  #include "glpnpp.h"
798*  int npp_is_covering(NPP *npp, NPPROW *row);
799*
800*  RETURNS
801*
802*  If the specified row (constraint) is covering inequality (see below),
803*  the routine npp_is_covering returns non-zero. Otherwise, it returns
804*  zero.
805*
806*  COVERING INEQUALITIES
807*
808*  In canonical format the covering inequality is the following:
809*
810*     sum  x[j] >= 1,                                                (1)
811*    j in J
812*
813*  where all variables x[j] are binary. This inequality expresses the
814*  condition that in any integer feasible solution variables in set J
815*  cannot be all equal to zero at the same time, i.e. at least one
816*  variable must take non-zero (unity) value. W.l.o.g. it is assumed
817*  that |J| >= 2, because if J is empty, the inequality (1) is
818*  infeasible, and if |J| = 1, the inequality (1) is a forcing row.
819*
820*  In general case the covering inequality may include original
821*  variables x[j] as well as their complements x~[j]:
822*
823*     sum   x[j] + sum   x~[j] >= 1,                                 (2)
824*    j in Jp      j in Jn
825*
826*  where Jp and Jn are not intersected. Therefore, using substitution
827*  x~[j] = 1 - x[j] gives the packing inequality in generalized format:
828*
829*     sum   x[j] - sum   x[j] >= 1 - |Jn|.                           (3)
830*    j in Jp      j in Jn
831*
832*  (May note that the inequality (3) cuts off infeasible solutions,
833*  where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.)
834*
835*  NOTE: If |J| = 2, the inequality (3) is equivalent to packing
836*        inequality (see the routine npp_is_packing). */
837
838int npp_is_covering(NPP *npp, NPPROW *row)
839{     /* test if constraint is covering inequality */
840      NPPCOL *col;
841      NPPAIJ *aij;
842      int b;
843      xassert(npp == npp);
844      if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX))
845         return 0;
846      b = 1;
847      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
848      {  col = aij->col;
849         if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
850            return 0;
851         if (aij->val == +1.0)
852            ;
853         else if (aij->val == -1.0)
854            b--;
855         else
856            return 0;
857      }
858      if (row->lb != (double)b) return 0;
859      return 1;
860}
861
862/***********************************************************************
863*  NAME
864*
865*  npp_hidden_covering - identify hidden covering inequality
866*
867*  SYNOPSIS
868*
869*  #include "glpnpp.h"
870*  int npp_hidden_covering(NPP *npp, NPPROW *row);
871*
872*  DESCRIPTION
873*
874*  The routine npp_hidden_covering processes specified inequality
875*  constraint, which includes only binary variables, and the number of
876*  the variables is not less than three. If the original inequality is
877*  equivalent to a covering inequality (see below), the routine
878*  replaces it by the equivalent inequality. If the original constraint
879*  is double-sided inequality, it is replaced by a pair of single-sided
880*  inequalities, if necessary.
881*
882*  RETURNS
883*
884*  If the original inequality constraint was replaced by equivalent
885*  covering inequality, the routine npp_hidden_covering returns
886*  non-zero. Otherwise, it returns zero.
887*
888*  PROBLEM TRANSFORMATION
889*
890*  Consider an inequality constraint:
891*
892*     sum  a[j] x[j] >= b,                                           (1)
893*    j in J
894*
895*  where all variables x[j] are binary, and |J| >= 3. (In case of '<='
896*  inequality it can be transformed to '>=' format by multiplying both
897*  its sides by -1.)
898*
899*  Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
900*  x[j] = 1 - x~[j] for all j in Jn, we have:
901*
902*     sum   a[j] x[j] >= b  ==>
903*    j in J
904*
905*     sum   a[j] x[j] + sum   a[j] x[j] >= b  ==>
906*    j in Jp           j in Jn
907*
908*     sum   a[j] x[j] + sum   a[j] (1 - x~[j]) >= b  ==>
909*    j in Jp           j in Jn
910*
911*     sum  m   a[j] x[j] - sum   a[j] x~[j] >= b - sum   a[j].
912*    j in Jp              j in Jn                 j in Jn
913*
914*  Thus, meaning the transformation above, we can assume that in
915*  inequality (1) all coefficients a[j] are positive. Moreover, we can
916*  assume that b > 0, because otherwise the inequality (1) would be
917*  redundant (see the routine npp_analyze_row). It is then obvious that
918*  constraint (1) is equivalent to covering inequality only if:
919*
920*     a[j] >= b,                                                     (2)
921*
922*  for all j in J.
923*
924*  Once the original inequality (1) is replaced by equivalent covering
925*  inequality, we need to perform back substitution x~[j] = 1 - x[j] for
926*  all j in Jn (see above).
927*
928*  RECOVERING SOLUTION
929*
930*  None needed. */
931
932static int hidden_covering(NPP *npp, struct elem *ptr, double *_b)
933{     /* process inequality constraint: sum a[j] x[j] >= b;
934         0 - specified row is NOT hidden covering inequality;
935         1 - specified row is covering inequality;
936         2 - specified row is hidden covering inequality. */
937      struct elem *e;
938      int neg;
939      double b = *_b, eps;
940      xassert(npp == npp);
941      /* a[j] must be non-zero, x[j] must be binary, for all j in J */
942      for (e = ptr; e != NULL; e = e->next)
943      {  xassert(e->aj != 0.0);
944         xassert(e->xj->is_int);
945         xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
946      }
947      /* check if the specified inequality constraint already has the
948         form of covering inequality */
949      neg = 0; /* neg is |Jn| */
950      for (e = ptr; e != NULL; e = e->next)
951      {  if (e->aj == +1.0)
952            ;
953         else if (e->aj == -1.0)
954            neg++;
955         else
956            break;
957      }
958      if (e == NULL)
959      {  /* all coefficients a[j] are +1 or -1; check rhs b */
960         if (b == (double)(1 - neg))
961         {  /* it is covering inequality; no processing is needed */
962            return 1;
963         }
964      }
965      /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
966         positive; the result is a~[j] = |a[j]| and new rhs b */
967      for (e = ptr; e != NULL; e = e->next)
968         if (e->aj < 0) b -= e->aj;
969      /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
970      /* if b <= 0, skip processing--this case must not appear */
971      if (b < 1e-3) return 0;
972      /* now a[j] > 0 for all j in J, and b > 0 */
973      /* the specified constraint is equivalent to covering inequality
974         iff a[j] >= b for all j in J */
975      eps = 1e-9 + 1e-12 * fabs(b);
976      for (e = ptr; e != NULL; e = e->next)
977         if (fabs(e->aj) < b - eps) return 0;
978      /* perform back substitution x~[j] = 1 - x[j] and construct the
979         final equivalent covering inequality in generalized format */
980      b = 1.0;
981      for (e = ptr; e != NULL; e = e->next)
982      {  if (e->aj > 0.0)
983            e->aj = +1.0;
984         else /* e->aj < 0.0 */
985            e->aj = -1.0, b -= 1.0;
986      }
987      *_b = b;
988      return 2;
989}
990
991int npp_hidden_covering(NPP *npp, NPPROW *row)
992{     /* identify hidden covering inequality */
993      NPPROW *copy;
994      NPPAIJ *aij;
995      struct elem *ptr, *e;
996      int kase, ret, count = 0;
997      double b;
998      /* the row must be inequality constraint */
999      xassert(row->lb < row->ub);
1000      for (kase = 0; kase <= 1; kase++)
1001      {  if (kase == 0)
1002         {  /* process row lower bound */
1003            if (row->lb == -DBL_MAX) continue;
1004            ptr = copy_form(npp, row, +1.0);
1005            b = + row->lb;
1006         }
1007         else
1008         {  /* process row upper bound */
1009            if (row->ub == +DBL_MAX) continue;
1010            ptr = copy_form(npp, row, -1.0);
1011            b = - row->ub;
1012         }
1013         /* now the inequality has the form "sum a[j] x[j] >= b" */
1014         ret = hidden_covering(npp, ptr, &b);
1015         xassert(0 <= ret && ret <= 2);
1016         if (kase == 1 && ret == 1 || ret == 2)
1017         {  /* the original inequality has been identified as hidden
1018               covering inequality */
1019            count++;
1020#ifdef GLP_DEBUG
1021            xprintf("Original constraint:\n");
1022            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1023               xprintf(" %+g x%d", aij->val, aij->col->j);
1024            if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
1025            if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
1026            xprintf("\n");
1027            xprintf("Equivalent covering inequality:\n");
1028            for (e = ptr; e != NULL; e = e->next)
1029               xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
1030            xprintf(", >= %g\n", b);
1031#endif
1032            if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
1033            {  /* the original row is single-sided inequality; no copy
1034                  is needed */
1035               copy = NULL;
1036            }
1037            else
1038            {  /* the original row is double-sided inequality; we need
1039                  to create its copy for other bound before replacing it
1040                  with the equivalent inequality */
1041               copy = npp_add_row(npp);
1042               if (kase == 0)
1043               {  /* the copy is for upper bound */
1044                  copy->lb = -DBL_MAX, copy->ub = row->ub;
1045               }
1046               else
1047               {  /* the copy is for lower bound */
1048                  copy->lb = row->lb, copy->ub = +DBL_MAX;
1049               }
1050               /* copy original row coefficients */
1051               for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1052                  npp_add_aij(npp, copy, aij->col, aij->val);
1053            }
1054            /* replace the original inequality by equivalent one */
1055            npp_erase_row(npp, row);
1056            row->lb = b, row->ub = +DBL_MAX;
1057            for (e = ptr; e != NULL; e = e->next)
1058               npp_add_aij(npp, row, e->xj, e->aj);
1059            /* continue processing upper bound for the copy */
1060            if (copy != NULL) row = copy;
1061         }
1062         drop_form(npp, ptr);
1063      }
1064      return count;
1065}
1066
1067/***********************************************************************
1068*  NAME
1069*
1070*  npp_is_partitioning - test if constraint is partitioning equality
1071*
1072*  SYNOPSIS
1073*
1074*  #include "glpnpp.h"
1075*  int npp_is_partitioning(NPP *npp, NPPROW *row);
1076*
1077*  RETURNS
1078*
1079*  If the specified row (constraint) is partitioning equality (see
1080*  below), the routine npp_is_partitioning returns non-zero. Otherwise,
1081*  it returns zero.
1082*
1083*  PARTITIONING EQUALITIES
1084*
1085*  In canonical format the partitioning equality is the following:
1086*
1087*     sum  x[j] = 1,                                                 (1)
1088*    j in J
1089*
1090*  where all variables x[j] are binary. This equality expresses the
1091*  condition that in any integer feasible solution exactly one variable
1092*  in set J must take non-zero (unity) value while other variables must
1093*  be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if
1094*  J is empty, the inequality (1) is infeasible, and if |J| = 1, the
1095*  inequality (1) is a fixing row.
1096*
1097*  In general case the partitioning equality may include original
1098*  variables x[j] as well as their complements x~[j]:
1099*
1100*     sum   x[j] + sum   x~[j] = 1,                                  (2)
1101*    j in Jp      j in Jn
1102*
1103*  where Jp and Jn are not intersected. Therefore, using substitution
1104*  x~[j] = 1 - x[j] leads to the partitioning equality in generalized
1105*  format:
1106*
1107*     sum   x[j] - sum   x[j] = 1 - |Jn|.                            (3)
1108*    j in Jp      j in Jn */
1109
1110int npp_is_partitioning(NPP *npp, NPPROW *row)
1111{     /* test if constraint is partitioning equality */
1112      NPPCOL *col;
1113      NPPAIJ *aij;
1114      int b;
1115      xassert(npp == npp);
1116      if (row->lb != row->ub) return 0;
1117      b = 1;
1118      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1119      {  col = aij->col;
1120         if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
1121            return 0;
1122         if (aij->val == +1.0)
1123            ;
1124         else if (aij->val == -1.0)
1125            b--;
1126         else
1127            return 0;
1128      }
1129      if (row->lb != (double)b) return 0;
1130      return 1;
1131}
1132
1133/***********************************************************************
1134*  NAME
1135*
1136*  npp_reduce_ineq_coef - reduce inequality constraint coefficients
1137*
1138*  SYNOPSIS
1139*
1140*  #include "glpnpp.h"
1141*  int npp_reduce_ineq_coef(NPP *npp, NPPROW *row);
1142*
1143*  DESCRIPTION
1144*
1145*  The routine npp_reduce_ineq_coef processes specified inequality
1146*  constraint attempting to replace it by an equivalent constraint,
1147*  where magnitude of coefficients at binary variables is smaller than
1148*  in the original constraint. If the inequality is double-sided, it is
1149*  replaced by a pair of single-sided inequalities, if necessary.
1150*
1151*  RETURNS
1152*
1153*  The routine npp_reduce_ineq_coef returns the number of coefficients
1154*  reduced.
1155*
1156*  BACKGROUND
1157*
1158*  Consider an inequality constraint:
1159*
1160*     sum  a[j] x[j] >= b.                                           (1)
1161*    j in J
1162*
1163*  (In case of '<=' inequality it can be transformed to '>=' format by
1164*  multiplying both its sides by -1.) Let x[k] be a binary variable;
1165*  other variables can be integer as well as continuous. We can write
1166*  constraint (1) as follows:
1167*
1168*     a[k] x[k] + t[k] >= b,                                         (2)
1169*
1170*  where:
1171*
1172*     t[k] = sum      a[j] x[j].                                     (3)
1173*           j in J\{k}
1174*
1175*  Since x[k] is binary, constraint (2) is equivalent to disjunction of
1176*  the following two constraints:
1177*
1178*     x[k] = 0,  t[k] >= b                                           (4)
1179*
1180*        OR
1181*
1182*     x[k] = 1,  t[k] >= b - a[k].                                   (5)
1183*
1184*  Let also that for the partial sum t[k] be known some its implied
1185*  lower bound inf t[k].
1186*
1187*  Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints
1188*  (4) and (5) and therefore constraint (2) are redundant.
1189*  If inf t[k] > b - a[k], only constraint (5) is redundant, in which
1190*  case it can be replaced with the following redundant and therefore
1191*  equivalent constraint:
1192*
1193*     t[k] >= b - a'[k] = inf t[k],                                  (6)
1194*
1195*  where:
1196*
1197*     a'[k] = b - inf t[k].                                          (7)
1198*
1199*  Thus, the original constraint (2) is equivalent to the following
1200*  constraint with coefficient at variable x[k] changed:
1201*
1202*     a'[k] x[k] + t[k] >= b.                                        (8)
1203*
1204*  From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient
1205*  at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that
1206*  a'[k] < a[k], i.e. the coefficient reduces in magnitude.
1207*
1208*  Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both
1209*  constraints (4) and (5) and therefore constraint (2) are redundant.
1210*  If inf t[k] > b, only constraint (4) is redundant, in which case it
1211*  can be replaced with the following redundant and therefore equivalent
1212*  constraint:
1213*
1214*     t[k] >= b' = inf t[k].                                         (9)
1215*
1216*  Rewriting constraint (5) as follows:
1217*
1218*     t[k] >= b - a[k] = b' - a'[k],                                (10)
1219*
1220*  where:
1221*
1222*     a'[k] = a[k] + b' - b = a[k] + inf t[k] - b,                  (11)
1223*
1224*  we can see that disjunction of constraint (9) and (10) is equivalent
1225*  to disjunction of constraint (4) and (5), from which it follows that
1226*  the original constraint (2) is equivalent to the following constraint
1227*  with both coefficient at variable x[k] and right-hand side changed:
1228*
1229*     a'[k] x[k] + t[k] >= b'.                                      (12)
1230*
1231*  From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the
1232*  coefficient at x[k] keeps its sign. And from inf t[k] > b it follows
1233*  that a'[k] > a[k], i.e. the coefficient reduces in magnitude.
1234*
1235*  PROBLEM TRANSFORMATION
1236*
1237*  In the routine npp_reduce_ineq_coef the following implied lower
1238*  bound of the partial sum (3) is used:
1239*
1240*     inf t[k] = sum       a[j] l[j] + sum       a[j] u[j],         (13)
1241*               j in Jp\{k}           k in Jn\{k}
1242*
1243*  where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are
1244*  lower and upper bounds, resp., of variable x[j].
1245*
1246*  In order to compute inf t[k] more efficiently, the following formula,
1247*  which is equivalent to (13), is actually used:
1248*
1249*                ( h - a[k] l[k] = h,        if a[k] > 0,
1250*     inf t[k] = <                                                  (14)
1251*                ( h - a[k] u[k] = h - a[k], if a[k] < 0,
1252*
1253*  where:
1254*
1255*     h = sum   a[j] l[j] + sum   a[j] u[j]                         (15)
1256*        j in Jp           j in Jn
1257*
1258*  is the implied lower bound of row (1).
1259*
1260*  Reduction of positive coefficient (a[k] > 0) does not change value
1261*  of h, since l[k] = 0. In case of reduction of negative coefficient
1262*  (a[k] < 0) from (11) it follows that:
1263*
1264*     delta a[k] = a'[k] - a[k] = inf t[k] - b  (> 0),              (16)
1265*
1266*  so new value of h (accounting that u[k] = 1) can be computed as
1267*  follows:
1268*
1269*     h := h + delta a[k] = h + (inf t[k] - b).                     (17)
1270*
1271*  RECOVERING SOLUTION
1272*
1273*  None needed. */
1274
1275static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b)
1276{     /* process inequality constraint: sum a[j] x[j] >= b */
1277      /* returns: the number of coefficients reduced */
1278      struct elem *e;
1279      int count = 0;
1280      double h, inf_t, new_a, b = *_b;
1281      xassert(npp == npp);
1282      /* compute h; see (15) */
1283      h = 0.0;
1284      for (e = ptr; e != NULL; e = e->next)
1285      {  if (e->aj > 0.0)
1286         {  if (e->xj->lb == -DBL_MAX) goto done;
1287            h += e->aj * e->xj->lb;
1288         }
1289         else /* e->aj < 0.0 */
1290         {  if (e->xj->ub == +DBL_MAX) goto done;
1291            h += e->aj * e->xj->ub;
1292         }
1293      }
1294      /* perform reduction of coefficients at binary variables */
1295      for (e = ptr; e != NULL; e = e->next)
1296      {  /* skip non-binary variable */
1297         if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
1298            continue;
1299         if (e->aj > 0.0)
1300         {  /* compute inf t[k]; see (14) */
1301            inf_t = h;
1302            if (b - e->aj < inf_t && inf_t < b)
1303            {  /* compute reduced coefficient a'[k]; see (7) */
1304               new_a = b - inf_t;
1305               if (new_a >= +1e-3 &&
1306                   e->aj - new_a >= 0.01 * (1.0 + e->aj))
1307               {  /* accept a'[k] */
1308#ifdef GLP_DEBUG
1309                  xprintf("+");
1310#endif
1311                  e->aj = new_a;
1312                  count++;
1313               }
1314            }
1315         }
1316         else /* e->aj < 0.0 */
1317         {  /* compute inf t[k]; see (14) */
1318            inf_t = h - e->aj;
1319            if (b < inf_t && inf_t < b - e->aj)
1320            {  /* compute reduced coefficient a'[k]; see (11) */
1321               new_a = e->aj + (inf_t - b);
1322               if (new_a <= -1e-3 &&
1323                   new_a - e->aj >= 0.01 * (1.0 - e->aj))
1324               {  /* accept a'[k] */
1325#ifdef GLP_DEBUG
1326                  xprintf("-");
1327#endif
1328                  e->aj = new_a;
1329                  /* update h; see (17) */
1330                  h += (inf_t - b);
1331                  /* compute b'; see (9) */
1332                  b = inf_t;
1333                  count++;
1334               }
1335            }
1336         }
1337      }
1338      *_b = b;
1339done: return count;
1340}
1341
1342int npp_reduce_ineq_coef(NPP *npp, NPPROW *row)
1343{     /* reduce inequality constraint coefficients */
1344      NPPROW *copy;
1345      NPPAIJ *aij;
1346      struct elem *ptr, *e;
1347      int kase, count[2];
1348      double b;
1349      /* the row must be inequality constraint */
1350      xassert(row->lb < row->ub);
1351      count[0] = count[1] = 0;
1352      for (kase = 0; kase <= 1; kase++)
1353      {  if (kase == 0)
1354         {  /* process row lower bound */
1355            if (row->lb == -DBL_MAX) continue;
1356#ifdef GLP_DEBUG
1357            xprintf("L");
1358#endif
1359            ptr = copy_form(npp, row, +1.0);
1360            b = + row->lb;
1361         }
1362         else
1363         {  /* process row upper bound */
1364            if (row->ub == +DBL_MAX) continue;
1365#ifdef GLP_DEBUG
1366            xprintf("U");
1367#endif
1368            ptr = copy_form(npp, row, -1.0);
1369            b = - row->ub;
1370         }
1371         /* now the inequality has the form "sum a[j] x[j] >= b" */
1372         count[kase] = reduce_ineq_coef(npp, ptr, &b);
1373         if (count[kase] > 0)
1374         {  /* the original inequality has been replaced by equivalent
1375               one with coefficients reduced */
1376            if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
1377            {  /* the original row is single-sided inequality; no copy
1378                  is needed */
1379               copy = NULL;
1380            }
1381            else
1382            {  /* the original row is double-sided inequality; we need
1383                  to create its copy for other bound before replacing it
1384                  with the equivalent inequality */
1385#ifdef GLP_DEBUG
1386               xprintf("*");
1387#endif
1388               copy = npp_add_row(npp);
1389               if (kase == 0)
1390               {  /* the copy is for upper bound */
1391                  copy->lb = -DBL_MAX, copy->ub = row->ub;
1392               }
1393               else
1394               {  /* the copy is for lower bound */
1395                  copy->lb = row->lb, copy->ub = +DBL_MAX;
1396               }
1397               /* copy original row coefficients */
1398               for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1399                  npp_add_aij(npp, copy, aij->col, aij->val);
1400            }
1401            /* replace the original inequality by equivalent one */
1402            npp_erase_row(npp, row);
1403            row->lb = b, row->ub = +DBL_MAX;
1404            for (e = ptr; e != NULL; e = e->next)
1405               npp_add_aij(npp, row, e->xj, e->aj);
1406            /* continue processing upper bound for the copy */
1407            if (copy != NULL) row = copy;
1408         }
1409         drop_form(npp, ptr);
1410      }
1411      return count[0] + count[1];
1412}
1413
1414/* eof */
Note: See TracBrowser for help on using the repository browser.