COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios11.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: 10.5 KB
Line 
1/* glpios11.c (process cuts stored in the local cut pool) */
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 "glpios.h"
26
27/***********************************************************************
28*  NAME
29*
30*  ios_process_cuts - process cuts stored in the local cut pool
31*
32*  SYNOPSIS
33*
34*  #include "glpios.h"
35*  void ios_process_cuts(glp_tree *T);
36*
37*  DESCRIPTION
38*
39*  The routine ios_process_cuts analyzes each cut currently stored in
40*  the local cut pool, which must be non-empty, and either adds the cut
41*  to the current subproblem or just discards it. All cuts are assumed
42*  to be locally valid. On exit the local cut pool remains unchanged.
43*
44*  REFERENCES
45*
46*  1. E.Balas, S.Ceria, G.Cornuejols, "Mixed 0-1 Programming by
47*     Lift-and-Project in a Branch-and-Cut Framework", Management Sc.,
48*     42 (1996) 1229-1246.
49*
50*  2. G.Andreello, A.Caprara, and M.Fischetti, "Embedding Cuts in
51*     a Branch&Cut Framework: a Computational Study with {0,1/2}-Cuts",
52*     Preliminary Draft, October 28, 2003, pp.6-8. */
53
54struct info
55{     /* estimated cut efficiency */
56      IOSCUT *cut;
57      /* pointer to cut in the cut pool */
58      char flag;
59      /* if this flag is set, the cut is included into the current
60         subproblem */
61      double eff;
62      /* cut efficacy (normalized residual) */
63      double deg;
64      /* lower bound to objective degradation */
65};
66
67static int fcmp(const void *arg1, const void *arg2)
68{     const struct info *info1 = arg1, *info2 = arg2;
69      if (info1->deg == 0.0 && info2->deg == 0.0)
70      {  if (info1->eff > info2->eff) return -1;
71         if (info1->eff < info2->eff) return +1;
72      }
73      else
74      {  if (info1->deg > info2->deg) return -1;
75         if (info1->deg < info2->deg) return +1;
76      }
77      return 0;
78}
79
80static double parallel(IOSCUT *a, IOSCUT *b, double work[]);
81
82void ios_process_cuts(glp_tree *T)
83{     IOSPOOL *pool;
84      IOSCUT *cut;
85      IOSAIJ *aij;
86      struct info *info;
87      int k, kk, max_cuts, len, ret, *ind;
88      double *val, *work;
89      /* the current subproblem must exist */
90      xassert(T->curr != NULL);
91      /* the pool must exist and be non-empty */
92      pool = T->local;
93      xassert(pool != NULL);
94      xassert(pool->size > 0);
95      /* allocate working arrays */
96      info = xcalloc(1+pool->size, sizeof(struct info));
97      ind = xcalloc(1+T->n, sizeof(int));
98      val = xcalloc(1+T->n, sizeof(double));
99      work = xcalloc(1+T->n, sizeof(double));
100      for (k = 1; k <= T->n; k++) work[k] = 0.0;
101      /* build the list of cuts stored in the cut pool */
102      for (k = 0, cut = pool->head; cut != NULL; cut = cut->next)
103         k++, info[k].cut = cut, info[k].flag = 0;
104      xassert(k == pool->size);
105      /* estimate efficiency of all cuts in the cut pool */
106      for (k = 1; k <= pool->size; k++)
107      {  double temp, dy, dz;
108         cut = info[k].cut;
109         /* build the vector of cut coefficients and compute its
110            Euclidean norm */
111         len = 0; temp = 0.0;
112         for (aij = cut->ptr; aij != NULL; aij = aij->next)
113         {  xassert(1 <= aij->j && aij->j <= T->n);
114            len++, ind[len] = aij->j, val[len] = aij->val;
115            temp += aij->val * aij->val;
116         }
117         if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
118         /* transform the cut to express it only through non-basic
119            (auxiliary and structural) variables */
120         len = glp_transform_row(T->mip, len, ind, val);
121         /* determine change in the cut value and in the objective
122            value for the adjacent basis by simulating one step of the
123            dual simplex */
124         ret = _glp_analyze_row(T->mip, len, ind, val, cut->type,
125            cut->rhs, 1e-9, NULL, NULL, NULL, NULL, &dy, &dz);
126         /* determine normalized residual and lower bound to objective
127            degradation */
128         if (ret == 0)
129         {  info[k].eff = fabs(dy) / sqrt(temp);
130            /* if some reduced costs violates (slightly) their zero
131               bounds (i.e. have wrong signs) due to round-off errors,
132               dz also may have wrong sign being close to zero */
133            if (T->mip->dir == GLP_MIN)
134            {  if (dz < 0.0) dz = 0.0;
135               info[k].deg = + dz;
136            }
137            else /* GLP_MAX */
138            {  if (dz > 0.0) dz = 0.0;
139               info[k].deg = - dz;
140            }
141         }
142         else if (ret == 1)
143         {  /* the constraint is not violated at the current point */
144            info[k].eff = info[k].deg = 0.0;
145         }
146         else if (ret == 2)
147         {  /* no dual feasible adjacent basis exists */
148            info[k].eff = 1.0;
149            info[k].deg = DBL_MAX;
150         }
151         else
152            xassert(ret != ret);
153         /* if the degradation is too small, just ignore it */
154         if (info[k].deg < 0.01) info[k].deg = 0.0;
155      }
156      /* sort the list of cuts by decreasing objective degradation and
157         then by decreasing efficacy */
158      qsort(&info[1], pool->size, sizeof(struct info), fcmp);
159      /* only first (most efficient) max_cuts in the list are qualified
160         as candidates to be added to the current subproblem */
161      max_cuts = (T->curr->level == 0 ? 90 : 10);
162      if (max_cuts > pool->size) max_cuts = pool->size;
163      /* add cuts to the current subproblem */
164#if 0
165      xprintf("*** adding cuts ***\n");
166#endif
167      for (k = 1; k <= max_cuts; k++)
168      {  int i, len;
169         /* if this cut seems to be inefficient, skip it */
170         if (info[k].deg < 0.01 && info[k].eff < 0.01) continue;
171         /* if the angle between this cut and every other cut included
172            in the current subproblem is small, skip this cut */
173         for (kk = 1; kk < k; kk++)
174         {  if (info[kk].flag)
175            {  if (parallel(info[k].cut, info[kk].cut, work) > 0.90)
176                  break;
177            }
178         }
179         if (kk < k) continue;
180         /* add this cut to the current subproblem */
181#if 0
182         xprintf("eff = %g; deg = %g\n", info[k].eff, info[k].deg);
183#endif
184         cut = info[k].cut, info[k].flag = 1;
185         i = glp_add_rows(T->mip, 1);
186         if (cut->name != NULL)
187            glp_set_row_name(T->mip, i, cut->name);
188         xassert(T->mip->row[i]->origin == GLP_RF_CUT);
189         T->mip->row[i]->klass = cut->klass;
190         len = 0;
191         for (aij = cut->ptr; aij != NULL; aij = aij->next)
192            len++, ind[len] = aij->j, val[len] = aij->val;
193         glp_set_mat_row(T->mip, i, len, ind, val);
194         xassert(cut->type == GLP_LO || cut->type == GLP_UP);
195         glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);
196      }
197      /* free working arrays */
198      xfree(info);
199      xfree(ind);
200      xfree(val);
201      xfree(work);
202      return;
203}
204
205#if 0
206/***********************************************************************
207*  Given a cut a * x >= b (<= b) the routine efficacy computes the cut
208*  efficacy as follows:
209*
210*     eff = d * (a * x~ - b) / ||a||,
211*
212*  where d is -1 (in case of '>= b') or +1 (in case of '<= b'), x~ is
213*  the vector of values of structural variables in optimal solution to
214*  LP relaxation of the current subproblem, ||a|| is the Euclidean norm
215*  of the vector of cut coefficients.
216*
217*  If the cut is violated at point x~, the efficacy eff is positive,
218*  and its value is the Euclidean distance between x~ and the cut plane
219*  a * x = b in the space of structural variables.
220*
221*  Following geometrical intuition, it is quite natural to consider
222*  this distance as a first-order measure of the expected efficacy of
223*  the cut: the larger the distance the better the cut [1]. */
224
225static double efficacy(glp_tree *T, IOSCUT *cut)
226{     glp_prob *mip = T->mip;
227      IOSAIJ *aij;
228      double s = 0.0, t = 0.0, temp;
229      for (aij = cut->ptr; aij != NULL; aij = aij->next)
230      {  xassert(1 <= aij->j && aij->j <= mip->n);
231         s += aij->val * mip->col[aij->j]->prim;
232         t += aij->val * aij->val;
233      }
234      temp = sqrt(t);
235      if (temp < DBL_EPSILON) temp = DBL_EPSILON;
236      if (cut->type == GLP_LO)
237         temp = (s >= cut->rhs ? 0.0 : (cut->rhs - s) / temp);
238      else if (cut->type == GLP_UP)
239         temp = (s <= cut->rhs ? 0.0 : (s - cut->rhs) / temp);
240      else
241         xassert(cut != cut);
242      return temp;
243}
244#endif
245
246/***********************************************************************
247*  Given two cuts a1 * x >= b1 (<= b1) and a2 * x >= b2 (<= b2) the
248*  routine parallel computes the cosine of angle between the cut planes
249*  a1 * x = b1 and a2 * x = b2 (which is the acute angle between two
250*  normals to these planes) in the space of structural variables as
251*  follows:
252*
253*     cos phi = (a1' * a2) / (||a1|| * ||a2||),
254*
255*  where (a1' * a2) is a dot product of vectors of cut coefficients,
256*  ||a1|| and ||a2|| are Euclidean norms of vectors a1 and a2.
257*
258*  Note that requirement cos phi = 0 forces the cuts to be orthogonal,
259*  i.e. with disjoint support, while requirement cos phi <= 0.999 means
260*  only avoiding duplicate (parallel) cuts [1]. */
261
262static double parallel(IOSCUT *a, IOSCUT *b, double work[])
263{     IOSAIJ *aij;
264      double s = 0.0, sa = 0.0, sb = 0.0, temp;
265      for (aij = a->ptr; aij != NULL; aij = aij->next)
266      {  work[aij->j] = aij->val;
267         sa += aij->val * aij->val;
268      }
269      for (aij = b->ptr; aij != NULL; aij = aij->next)
270      {  s += work[aij->j] * aij->val;
271         sb += aij->val * aij->val;
272      }
273      for (aij = a->ptr; aij != NULL; aij = aij->next)
274         work[aij->j] = 0.0;
275      temp = sqrt(sa) * sqrt(sb);
276      if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
277      return s / temp;
278}
279
280/* eof */
Note: See TracBrowser for help on using the repository browser.