COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios03.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: 42.2 KB
Line 
1/* glpios03.c (branch-and-cut driver) */
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*  show_progress - display current progress of the search
29*
30*  This routine displays some information about current progress of the
31*  search.
32*
33*  The information includes:
34*
35*  the current number of iterations performed by the simplex solver;
36*
37*  the objective value for the best known integer feasible solution,
38*  which is upper (minimization) or lower (maximization) global bound
39*  for optimal solution of the original mip problem;
40*
41*  the best local bound for active nodes, which is lower (minimization)
42*  or upper (maximization) global bound for optimal solution of the
43*  original mip problem;
44*
45*  the relative mip gap, in percents;
46*
47*  the number of open (active) subproblems;
48*
49*  the number of completely explored subproblems, i.e. whose nodes have
50*  been removed from the tree. */
51
52static void show_progress(glp_tree *T, int bingo)
53{     int p;
54      double temp;
55      char best_mip[50], best_bound[50], *rho, rel_gap[50];
56      /* format the best known integer feasible solution */
57      if (T->mip->mip_stat == GLP_FEAS)
58         sprintf(best_mip, "%17.9e", T->mip->mip_obj);
59      else
60         sprintf(best_mip, "%17s", "not found yet");
61      /* determine reference number of an active subproblem whose local
62         bound is best */
63      p = ios_best_node(T);
64      /* format the best bound */
65      if (p == 0)
66         sprintf(best_bound, "%17s", "tree is empty");
67      else
68      {  temp = T->slot[p].node->bound;
69         if (temp == -DBL_MAX)
70            sprintf(best_bound, "%17s", "-inf");
71         else if (temp == +DBL_MAX)
72            sprintf(best_bound, "%17s", "+inf");
73         else
74            sprintf(best_bound, "%17.9e", temp);
75      }
76      /* choose the relation sign between global bounds */
77      if (T->mip->dir == GLP_MIN)
78         rho = ">=";
79      else if (T->mip->dir == GLP_MAX)
80         rho = "<=";
81      else
82         xassert(T != T);
83      /* format the relative mip gap */
84      temp = ios_relative_gap(T);
85      if (temp == 0.0)
86         sprintf(rel_gap, "  0.0%%");
87      else if (temp < 0.001)
88         sprintf(rel_gap, "< 0.1%%");
89      else if (temp <= 9.999)
90         sprintf(rel_gap, "%5.1f%%", 100.0 * temp);
91      else
92         sprintf(rel_gap, "%6s", "");
93      /* display progress of the search */
94      xprintf("+%6d: %s %s %s %s %s (%d; %d)\n",
95         T->mip->it_cnt, bingo ? ">>>>>" : "mip =", best_mip, rho,
96         best_bound, rel_gap, T->a_cnt, T->t_cnt - T->n_cnt);
97      T->tm_lag = xtime();
98      return;
99}
100
101/***********************************************************************
102*  is_branch_hopeful - check if specified branch is hopeful
103*
104*  This routine checks if the specified subproblem can have an integer
105*  optimal solution which is better than the best known one.
106*
107*  The check is based on comparison of the local objective bound stored
108*  in the subproblem descriptor and the incumbent objective value which
109*  is the global objective bound.
110*
111*  If there is a chance that the specified subproblem can have a better
112*  integer optimal solution, the routine returns non-zero. Otherwise, if
113*  the corresponding branch can pruned, zero is returned. */
114
115static int is_branch_hopeful(glp_tree *T, int p)
116{     xassert(1 <= p && p <= T->nslots);
117      xassert(T->slot[p].node != NULL);
118      return ios_is_hopeful(T, T->slot[p].node->bound);
119}
120
121/***********************************************************************
122*  check_integrality - check integrality of basic solution
123*
124*  This routine checks if the basic solution of LP relaxation of the
125*  current subproblem satisfies to integrality conditions, i.e. that all
126*  variables of integer kind have integral primal values. (The solution
127*  is assumed to be optimal.)
128*
129*  For each variable of integer kind the routine computes the following
130*  quantity:
131*
132*     ii(x[j]) = min(x[j] - floor(x[j]), ceil(x[j]) - x[j]),         (1)
133*
134*  which is a measure of the integer infeasibility (non-integrality) of
135*  x[j] (for example, ii(2.1) = 0.1, ii(3.7) = 0.3, ii(5.0) = 0). It is
136*  understood that 0 <= ii(x[j]) <= 0.5, and variable x[j] is integer
137*  feasible if ii(x[j]) = 0. However, due to floating-point arithmetic
138*  the routine checks less restrictive condition:
139*
140*     ii(x[j]) <= tol_int,                                           (2)
141*
142*  where tol_int is a given tolerance (small positive number) and marks
143*  each variable which does not satisfy to (2) as integer infeasible by
144*  setting its fractionality flag.
145*
146*  In order to characterize integer infeasibility of the basic solution
147*  in the whole the routine computes two parameters: ii_cnt, which is
148*  the number of variables with the fractionality flag set, and ii_sum,
149*  which is the sum of integer infeasibilities (1). */
150
151static void check_integrality(glp_tree *T)
152{     glp_prob *mip = T->mip;
153      int j, type, ii_cnt = 0;
154      double lb, ub, x, temp1, temp2, ii_sum = 0.0;
155      /* walk through the set of columns (structural variables) */
156      for (j = 1; j <= mip->n; j++)
157      {  GLPCOL *col = mip->col[j];
158         T->non_int[j] = 0;
159         /* if the column is not integer, skip it */
160         if (col->kind != GLP_IV) continue;
161         /* if the column is non-basic, it is integer feasible */
162         if (col->stat != GLP_BS) continue;
163         /* obtain the type and bounds of the column */
164         type = col->type, lb = col->lb, ub = col->ub;
165         /* obtain value of the column in optimal basic solution */
166         x = col->prim;
167         /* if the column's primal value is close to the lower bound,
168            the column is integer feasible within given tolerance */
169         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
170         {  temp1 = lb - T->parm->tol_int;
171            temp2 = lb + T->parm->tol_int;
172            if (temp1 <= x && x <= temp2) continue;
173#if 0
174            /* the lower bound must not be violated */
175            xassert(x >= lb);
176#else
177            if (x < lb) continue;
178#endif
179         }
180         /* if the column's primal value is close to the upper bound,
181            the column is integer feasible within given tolerance */
182         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
183         {  temp1 = ub - T->parm->tol_int;
184            temp2 = ub + T->parm->tol_int;
185            if (temp1 <= x && x <= temp2) continue;
186#if 0
187            /* the upper bound must not be violated */
188            xassert(x <= ub);
189#else
190            if (x > ub) continue;
191#endif
192         }
193         /* if the column's primal value is close to nearest integer,
194            the column is integer feasible within given tolerance */
195         temp1 = floor(x + 0.5) - T->parm->tol_int;
196         temp2 = floor(x + 0.5) + T->parm->tol_int;
197         if (temp1 <= x && x <= temp2) continue;
198         /* otherwise the column is integer infeasible */
199         T->non_int[j] = 1;
200         /* increase the number of fractional-valued columns */
201         ii_cnt++;
202         /* compute the sum of integer infeasibilities */
203         temp1 = x - floor(x);
204         temp2 = ceil(x) - x;
205         xassert(temp1 > 0.0 && temp2 > 0.0);
206         ii_sum += (temp1 <= temp2 ? temp1 : temp2);
207      }
208      /* store ii_cnt and ii_sum to the current problem descriptor */
209      xassert(T->curr != NULL);
210      T->curr->ii_cnt = ii_cnt;
211      T->curr->ii_sum = ii_sum;
212      /* and also display these parameters */
213      if (T->parm->msg_lev >= GLP_MSG_DBG)
214      {  if (ii_cnt == 0)
215            xprintf("There are no fractional columns\n");
216         else if (ii_cnt == 1)
217            xprintf("There is one fractional column, integer infeasibil"
218               "ity is %.3e\n", ii_sum);
219         else
220            xprintf("There are %d fractional columns, integer infeasibi"
221               "lity is %.3e\n", ii_cnt, ii_sum);
222      }
223      return;
224}
225
226/***********************************************************************
227*  record_solution - record better integer feasible solution
228*
229*  This routine records optimal basic solution of LP relaxation of the
230*  current subproblem, which being integer feasible is better than the
231*  best known integer feasible solution. */
232
233static void record_solution(glp_tree *T)
234{     glp_prob *mip = T->mip;
235      int i, j;
236      mip->mip_stat = GLP_FEAS;
237      mip->mip_obj = mip->obj_val;
238      for (i = 1; i <= mip->m; i++)
239      {  GLPROW *row = mip->row[i];
240         row->mipx = row->prim;
241      }
242      for (j = 1; j <= mip->n; j++)
243      {  GLPCOL *col = mip->col[j];
244         if (col->kind == GLP_CV)
245            col->mipx = col->prim;
246         else if (col->kind == GLP_IV)
247         {  /* value of the integer column must be integral */
248            col->mipx = floor(col->prim + 0.5);
249         }
250         else
251            xassert(col != col);
252      }
253      T->sol_cnt++;
254      return;
255}
256
257/***********************************************************************
258*  fix_by_red_cost - fix non-basic integer columns by reduced costs
259*
260*  This routine fixes some non-basic integer columns if their reduced
261*  costs indicate that increasing (decreasing) the column at least by
262*  one involves the objective value becoming worse than the incumbent
263*  objective value. */
264
265static void fix_by_red_cost(glp_tree *T)
266{     glp_prob *mip = T->mip;
267      int j, stat, fixed = 0;
268      double obj, lb, ub, dj;
269      /* the global bound must exist */
270      xassert(T->mip->mip_stat == GLP_FEAS);
271      /* basic solution of LP relaxation must be optimal */
272      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
273      /* determine the objective function value */
274      obj = mip->obj_val;
275      /* walk through the column list */
276      for (j = 1; j <= mip->n; j++)
277      {  GLPCOL *col = mip->col[j];
278         /* if the column is not integer, skip it */
279         if (col->kind != GLP_IV) continue;
280         /* obtain bounds of j-th column */
281         lb = col->lb, ub = col->ub;
282         /* and determine its status and reduced cost */
283         stat = col->stat, dj = col->dual;
284         /* analyze the reduced cost */
285         switch (mip->dir)
286         {  case GLP_MIN:
287               /* minimization */
288               if (stat == GLP_NL)
289               {  /* j-th column is non-basic on its lower bound */
290                  if (dj < 0.0) dj = 0.0;
291                  if (obj + dj >= mip->mip_obj)
292                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
293               }
294               else if (stat == GLP_NU)
295               {  /* j-th column is non-basic on its upper bound */
296                  if (dj > 0.0) dj = 0.0;
297                  if (obj - dj >= mip->mip_obj)
298                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
299               }
300               break;
301            case GLP_MAX:
302               /* maximization */
303               if (stat == GLP_NL)
304               {  /* j-th column is non-basic on its lower bound */
305                  if (dj > 0.0) dj = 0.0;
306                  if (obj + dj <= mip->mip_obj)
307                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
308               }
309               else if (stat == GLP_NU)
310               {  /* j-th column is non-basic on its upper bound */
311                  if (dj < 0.0) dj = 0.0;
312                  if (obj - dj <= mip->mip_obj)
313                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
314               }
315               break;
316            default:
317               xassert(T != T);
318         }
319      }
320      if (T->parm->msg_lev >= GLP_MSG_DBG)
321      {  if (fixed == 0)
322            /* nothing to say */;
323         else if (fixed == 1)
324            xprintf("One column has been fixed by reduced cost\n");
325         else
326            xprintf("%d columns have been fixed by reduced costs\n",
327               fixed);
328      }
329      /* fixing non-basic columns on their current bounds does not
330         change the basic solution */
331      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
332      return;
333}
334
335/***********************************************************************
336*  branch_on - perform branching on specified variable
337*
338*  This routine performs branching on j-th column (structural variable)
339*  of the current subproblem. The specified column must be of integer
340*  kind and must have a fractional value in optimal basic solution of
341*  LP relaxation of the current subproblem (i.e. only columns for which
342*  the flag non_int[j] is set are valid candidates to branch on).
343*
344*  Let x be j-th structural variable, and beta be its primal fractional
345*  value in the current basic solution. Branching on j-th variable is
346*  dividing the current subproblem into two new subproblems, which are
347*  identical to the current subproblem with the following exception: in
348*  the first subproblem that begins the down-branch x has a new upper
349*  bound x <= floor(beta), and in the second subproblem that begins the
350*  up-branch x has a new lower bound x >= ceil(beta).
351*
352*  Depending on estimation of local bounds for down- and up-branches
353*  this routine returns the following:
354*
355*  0 - both branches have been created;
356*  1 - one branch is hopeless and has been pruned, so now the current
357*      subproblem is other branch;
358*  2 - both branches are hopeless and have been pruned; new subproblem
359*      selection is needed to continue the search. */
360
361static int branch_on(glp_tree *T, int j, int next)
362{     glp_prob *mip = T->mip;
363      IOSNPD *node;
364      int m = mip->m;
365      int n = mip->n;
366      int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
367      double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
368      /* determine bounds and value of x[j] in optimal solution to LP
369         relaxation of the current subproblem */
370      xassert(1 <= j && j <= n);
371      type = mip->col[j]->type;
372      lb = mip->col[j]->lb;
373      ub = mip->col[j]->ub;
374      beta = mip->col[j]->prim;
375      /* determine new bounds of x[j] for down- and up-branches */
376      new_ub = floor(beta);
377      new_lb = ceil(beta);
378      switch (type)
379      {  case GLP_FR:
380            dn_type = GLP_UP;
381            up_type = GLP_LO;
382            break;
383         case GLP_LO:
384            xassert(lb <= new_ub);
385            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
386            xassert(lb + 1.0 <= new_lb);
387            up_type = GLP_LO;
388            break;
389         case GLP_UP:
390            xassert(new_ub <= ub - 1.0);
391            dn_type = GLP_UP;
392            xassert(new_lb <= ub);
393            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
394            break;
395         case GLP_DB:
396            xassert(lb <= new_ub && new_ub <= ub - 1.0);
397            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
398            xassert(lb + 1.0 <= new_lb && new_lb <= ub);
399            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
400            break;
401         default:
402            xassert(type != type);
403      }
404      /* compute local bounds to LP relaxation for both branches */
405      ios_eval_degrad(T, j, &dn_lp, &up_lp);
406      /* and improve them by rounding */
407      dn_bnd = ios_round_bound(T, dn_lp);
408      up_bnd = ios_round_bound(T, up_lp);
409      /* check local bounds for down- and up-branches */
410      dn_bad = !ios_is_hopeful(T, dn_bnd);
411      up_bad = !ios_is_hopeful(T, up_bnd);
412      if (dn_bad && up_bad)
413      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
414            xprintf("Both down- and up-branches are hopeless\n");
415         ret = 2;
416         goto done;
417      }
418      else if (up_bad)
419      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
420            xprintf("Up-branch is hopeless\n");
421         glp_set_col_bnds(mip, j, dn_type, lb, new_ub);
422         T->curr->lp_obj = dn_lp;
423         if (mip->dir == GLP_MIN)
424         {  if (T->curr->bound < dn_bnd)
425                T->curr->bound = dn_bnd;
426         }
427         else if (mip->dir == GLP_MAX)
428         {  if (T->curr->bound > dn_bnd)
429                T->curr->bound = dn_bnd;
430         }
431         else
432            xassert(mip != mip);
433         ret = 1;
434         goto done;
435      }
436      else if (dn_bad)
437      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
438            xprintf("Down-branch is hopeless\n");
439         glp_set_col_bnds(mip, j, up_type, new_lb, ub);
440         T->curr->lp_obj = up_lp;
441         if (mip->dir == GLP_MIN)
442         {  if (T->curr->bound < up_bnd)
443                T->curr->bound = up_bnd;
444         }
445         else if (mip->dir == GLP_MAX)
446         {  if (T->curr->bound > up_bnd)
447                T->curr->bound = up_bnd;
448         }
449         else
450            xassert(mip != mip);
451         ret = 1;
452         goto done;
453      }
454      /* both down- and up-branches seem to be hopeful */
455      if (T->parm->msg_lev >= GLP_MSG_DBG)
456         xprintf("Branching on column %d, primal value is %.9e\n",
457            j, beta);
458      /* determine the reference number of the current subproblem */
459      xassert(T->curr != NULL);
460      p = T->curr->p;
461      T->curr->br_var = j;
462      T->curr->br_val = beta;
463      /* freeze the current subproblem */
464      ios_freeze_node(T);
465      /* create two clones of the current subproblem; the first clone
466         begins the down-branch, the second one begins the up-branch */
467      ios_clone_node(T, p, 2, clone);
468      if (T->parm->msg_lev >= GLP_MSG_DBG)
469         xprintf("Node %d begins down branch, node %d begins up branch "
470            "\n", clone[1], clone[2]);
471      /* set new upper bound of j-th column in the down-branch */
472      node = T->slot[clone[1]].node;
473      xassert(node != NULL);
474      xassert(node->up != NULL);
475      xassert(node->b_ptr == NULL);
476      node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
477      node->b_ptr->k = m + j;
478      node->b_ptr->type = (unsigned char)dn_type;
479      node->b_ptr->lb = lb;
480      node->b_ptr->ub = new_ub;
481      node->b_ptr->next = NULL;
482      node->lp_obj = dn_lp;
483      if (mip->dir == GLP_MIN)
484      {  if (node->bound < dn_bnd)
485             node->bound = dn_bnd;
486      }
487      else if (mip->dir == GLP_MAX)
488      {  if (node->bound > dn_bnd)
489             node->bound = dn_bnd;
490      }
491      else
492         xassert(mip != mip);
493      /* set new lower bound of j-th column in the up-branch */
494      node = T->slot[clone[2]].node;
495      xassert(node != NULL);
496      xassert(node->up != NULL);
497      xassert(node->b_ptr == NULL);
498      node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
499      node->b_ptr->k = m + j;
500      node->b_ptr->type = (unsigned char)up_type;
501      node->b_ptr->lb = new_lb;
502      node->b_ptr->ub = ub;
503      node->b_ptr->next = NULL;
504      node->lp_obj = up_lp;
505      if (mip->dir == GLP_MIN)
506      {  if (node->bound < up_bnd)
507             node->bound = up_bnd;
508      }
509      else if (mip->dir == GLP_MAX)
510      {  if (node->bound > up_bnd)
511             node->bound = up_bnd;
512      }
513      else
514         xassert(mip != mip);
515      /* suggest the subproblem to be solved next */
516      xassert(T->child == 0);
517      if (next == GLP_NO_BRNCH)
518         T->child = 0;
519      else if (next == GLP_DN_BRNCH)
520         T->child = clone[1];
521      else if (next == GLP_UP_BRNCH)
522         T->child = clone[2];
523      else
524         xassert(next != next);
525      ret = 0;
526done: return ret;
527}
528
529/***********************************************************************
530*  cleanup_the_tree - prune hopeless branches from the tree
531*
532*  This routine walks through the active list and checks the local
533*  bound for every active subproblem. If the local bound indicates that
534*  the subproblem cannot have integer optimal solution better than the
535*  incumbent objective value, the routine deletes such subproblem that,
536*  in turn, involves pruning the corresponding branch of the tree. */
537
538static void cleanup_the_tree(glp_tree *T)
539{     IOSNPD *node, *next_node;
540      int count = 0;
541      /* the global bound must exist */
542      xassert(T->mip->mip_stat == GLP_FEAS);
543      /* walk through the list of active subproblems */
544      for (node = T->head; node != NULL; node = next_node)
545      {  /* deleting some active problem node may involve deleting its
546            parents recursively; however, all its parents being created
547            *before* it are always *precede* it in the node list, so
548            the next problem node is never affected by such deletion */
549         next_node = node->next;
550         /* if the branch is hopeless, prune it */
551         if (!is_branch_hopeful(T, node->p))
552            ios_delete_node(T, node->p), count++;
553      }
554      if (T->parm->msg_lev >= GLP_MSG_DBG)
555      {  if (count == 1)
556            xprintf("One hopeless branch has been pruned\n");
557         else if (count > 1)
558            xprintf("%d hopeless branches have been pruned\n", count);
559      }
560      return;
561}
562
563/**********************************************************************/
564
565static void generate_cuts(glp_tree *T)
566{     /* generate generic cuts with built-in generators */
567      if (!(T->parm->mir_cuts == GLP_ON ||
568            T->parm->gmi_cuts == GLP_ON ||
569            T->parm->cov_cuts == GLP_ON ||
570            T->parm->clq_cuts == GLP_ON)) goto done;
571#if 1 /* 20/IX-2008 */
572      {  int i, max_cuts, added_cuts;
573         max_cuts = T->n;
574         if (max_cuts < 1000) max_cuts = 1000;
575         added_cuts = 0;
576         for (i = T->orig_m+1; i <= T->mip->m; i++)
577         {  if (T->mip->row[i]->origin == GLP_RF_CUT)
578               added_cuts++;
579         }
580         /* xprintf("added_cuts = %d\n", added_cuts); */
581         if (added_cuts >= max_cuts) goto done;
582      }
583#endif
584      /* generate and add to POOL all cuts violated by x* */
585      if (T->parm->gmi_cuts == GLP_ON)
586      {  if (T->curr->changed < 5)
587            ios_gmi_gen(T);
588      }
589      if (T->parm->mir_cuts == GLP_ON)
590      {  xassert(T->mir_gen != NULL);
591         ios_mir_gen(T, T->mir_gen);
592      }
593      if (T->parm->cov_cuts == GLP_ON)
594      {  /* cover cuts works well along with mir cuts */
595         /*if (T->round <= 5)*/
596            ios_cov_gen(T);
597      }
598      if (T->parm->clq_cuts == GLP_ON)
599      {  if (T->clq_gen != NULL)
600         {  if (T->curr->level == 0 && T->curr->changed < 50 ||
601                T->curr->level >  0 && T->curr->changed < 5)
602               ios_clq_gen(T, T->clq_gen);
603         }
604      }
605done: return;
606}
607
608/**********************************************************************/
609
610static void remove_cuts(glp_tree *T)
611{     /* remove inactive cuts (some valueable globally valid cut might
612         be saved in the global cut pool) */
613      int i, cnt = 0, *num = NULL;
614      xassert(T->curr != NULL);
615      for (i = T->orig_m+1; i <= T->mip->m; i++)
616      {  if (T->mip->row[i]->origin == GLP_RF_CUT &&
617             T->mip->row[i]->level == T->curr->level &&
618             T->mip->row[i]->stat == GLP_BS)
619         {  if (num == NULL)
620               num = xcalloc(1+T->mip->m, sizeof(int));
621            num[++cnt] = i;
622         }
623      }
624      if (cnt > 0)
625      {  glp_del_rows(T->mip, cnt, num);
626#if 0
627         xprintf("%d inactive cut(s) removed\n", cnt);
628#endif
629         xfree(num);
630         xassert(glp_factorize(T->mip) == 0);
631      }
632      return;
633}
634
635/**********************************************************************/
636
637static void display_cut_info(glp_tree *T)
638{     glp_prob *mip = T->mip;
639      int i, gmi = 0, mir = 0, cov = 0, clq = 0, app = 0;
640      for (i = mip->m; i > 0; i--)
641      {  GLPROW *row;
642         row = mip->row[i];
643         /* if (row->level < T->curr->level) break; */
644         if (row->origin == GLP_RF_CUT)
645         {  if (row->klass == GLP_RF_GMI)
646               gmi++;
647            else if (row->klass == GLP_RF_MIR)
648               mir++;
649            else if (row->klass == GLP_RF_COV)
650               cov++;
651            else if (row->klass == GLP_RF_CLQ)
652               clq++;
653            else
654               app++;
655         }
656      }
657      xassert(T->curr != NULL);
658      if (gmi + mir + cov + clq + app > 0)
659      {  xprintf("Cuts on level %d:", T->curr->level);
660         if (gmi > 0) xprintf(" gmi = %d;", gmi);
661         if (mir > 0) xprintf(" mir = %d;", mir);
662         if (cov > 0) xprintf(" cov = %d;", cov);
663         if (clq > 0) xprintf(" clq = %d;", clq);
664         if (app > 0) xprintf(" app = %d;", app);
665         xprintf("\n");
666      }
667      return;
668}
669
670/***********************************************************************
671*  NAME
672*
673*  ios_driver - branch-and-cut driver
674*
675*  SYNOPSIS
676*
677*  #include "glpios.h"
678*  int ios_driver(glp_tree *T);
679*
680*  DESCRIPTION
681*
682*  The routine ios_driver is a branch-and-cut driver. It controls the
683*  MIP solution process.
684*
685*  RETURNS
686*
687*  0  The MIP problem instance has been successfully solved. This code
688*     does not necessarily mean that the solver has found optimal
689*     solution. It only means that the solution process was successful.
690*
691*  GLP_EFAIL
692*     The search was prematurely terminated due to the solver failure.
693*
694*  GLP_EMIPGAP
695*     The search was prematurely terminated, because the relative mip
696*     gap tolerance has been reached.
697*
698*  GLP_ETMLIM
699*     The search was prematurely terminated, because the time limit has
700*     been exceeded.
701*
702*  GLP_ESTOP
703*     The search was prematurely terminated by application. */
704
705int ios_driver(glp_tree *T)
706{     int p, curr_p, p_stat, d_stat, ret;
707#if 1 /* carry out to glp_tree */
708      int pred_p = 0;
709      /* if the current subproblem has been just created due to
710         branching, pred_p is the reference number of its parent
711         subproblem, otherwise pred_p is zero */
712#endif
713      glp_long ttt = T->tm_beg;
714#if 0
715      ((glp_iocp *)T->parm)->msg_lev = GLP_MSG_DBG;
716#endif
717      /* on entry to the B&B driver it is assumed that the active list
718         contains the only active (i.e. root) subproblem, which is the
719         original MIP problem to be solved */
720loop: /* main loop starts here */
721      /* at this point the current subproblem does not exist */
722      xassert(T->curr == NULL);
723      /* if the active list is empty, the search is finished */
724      if (T->head == NULL)
725      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
726            xprintf("Active list is empty!\n");
727         xassert(dmp_in_use(T->pool).lo == 0);
728         ret = 0;
729         goto done;
730      }
731      /* select some active subproblem to continue the search */
732      xassert(T->next_p == 0);
733      /* let the application program select subproblem */
734      if (T->parm->cb_func != NULL)
735      {  xassert(T->reason == 0);
736         T->reason = GLP_ISELECT;
737         T->parm->cb_func(T, T->parm->cb_info);
738         T->reason = 0;
739         if (T->stop)
740         {  ret = GLP_ESTOP;
741            goto done;
742         }
743      }
744      if (T->next_p != 0)
745      {  /* the application program has selected something */
746         ;
747      }
748      else if (T->a_cnt == 1)
749      {  /* the only active subproblem exists, so select it */
750         xassert(T->head->next == NULL);
751         T->next_p = T->head->p;
752      }
753      else if (T->child != 0)
754      {  /* select one of branching childs suggested by the branching
755            heuristic */
756         T->next_p = T->child;
757      }
758      else
759      {  /* select active subproblem as specified by the backtracking
760            technique option */
761         T->next_p = ios_choose_node(T);
762      }
763      /* the active subproblem just selected becomes current */
764      ios_revive_node(T, T->next_p);
765      T->next_p = T->child = 0;
766      /* invalidate pred_p, if it is not the reference number of the
767         parent of the current subproblem */
768      if (T->curr->up != NULL && T->curr->up->p != pred_p) pred_p = 0;
769      /* determine the reference number of the current subproblem */
770      p = T->curr->p;
771      if (T->parm->msg_lev >= GLP_MSG_DBG)
772      {  xprintf("-----------------------------------------------------"
773            "-------------------\n");
774         xprintf("Processing node %d at level %d\n", p, T->curr->level);
775      }
776      /* if it is the root subproblem, initialize cut generators */
777      if (p == 1)
778      {  if (T->parm->gmi_cuts == GLP_ON)
779         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
780               xprintf("Gomory's cuts enabled\n");
781         }
782         if (T->parm->mir_cuts == GLP_ON)
783         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
784               xprintf("MIR cuts enabled\n");
785            xassert(T->mir_gen == NULL);
786            T->mir_gen = ios_mir_init(T);
787         }
788         if (T->parm->cov_cuts == GLP_ON)
789         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
790               xprintf("Cover cuts enabled\n");
791         }
792         if (T->parm->clq_cuts == GLP_ON)
793         {  xassert(T->clq_gen == NULL);
794            if (T->parm->msg_lev >= GLP_MSG_ALL)
795               xprintf("Clique cuts enabled\n");
796            T->clq_gen = ios_clq_init(T);
797         }
798      }
799more: /* minor loop starts here */
800      /* at this point the current subproblem needs either to be solved
801         for the first time or re-optimized due to reformulation */
802      /* display current progress of the search */
803      if (T->parm->msg_lev >= GLP_MSG_DBG ||
804          T->parm->msg_lev >= GLP_MSG_ON &&
805        (double)(T->parm->out_frq - 1) <=
806            1000.0 * xdifftime(xtime(), T->tm_lag))
807         show_progress(T, 0);
808      if (T->parm->msg_lev >= GLP_MSG_ALL &&
809            xdifftime(xtime(), ttt) >= 60.0)
810      {  glp_long total;
811         glp_mem_usage(NULL, NULL, &total, NULL);
812         xprintf("Time used: %.1f secs.  Memory used: %.1f Mb.\n",
813            xdifftime(xtime(), T->tm_beg), xltod(total) / 1048576.0);
814         ttt = xtime();
815      }
816      /* check the mip gap */
817      if (T->parm->mip_gap > 0.0 &&
818          ios_relative_gap(T) <= T->parm->mip_gap)
819      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
820            xprintf("Relative gap tolerance reached; search terminated "
821               "\n");
822         ret = GLP_EMIPGAP;
823         goto done;
824      }
825      /* check if the time limit has been exhausted */
826      if (T->parm->tm_lim < INT_MAX &&
827         (double)(T->parm->tm_lim - 1) <=
828         1000.0 * xdifftime(xtime(), T->tm_beg))
829      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
830            xprintf("Time limit exhausted; search terminated\n");
831         ret = GLP_ETMLIM;
832         goto done;
833      }
834      /* let the application program preprocess the subproblem */
835      if (T->parm->cb_func != NULL)
836      {  xassert(T->reason == 0);
837         T->reason = GLP_IPREPRO;
838         T->parm->cb_func(T, T->parm->cb_info);
839         T->reason = 0;
840         if (T->stop)
841         {  ret = GLP_ESTOP;
842            goto done;
843         }
844      }
845      /* perform basic preprocessing */
846      if (T->parm->pp_tech == GLP_PP_NONE)
847         ;
848      else if (T->parm->pp_tech == GLP_PP_ROOT)
849      {  if (T->curr->level == 0)
850         {  if (ios_preprocess_node(T, 100))
851               goto fath;
852         }
853      }
854      else if (T->parm->pp_tech == GLP_PP_ALL)
855      {  if (ios_preprocess_node(T, T->curr->level == 0 ? 100 : 10))
856            goto fath;
857      }
858      else
859         xassert(T != T);
860      /* preprocessing may improve the global bound */
861      if (!is_branch_hopeful(T, p))
862      {  xprintf("*** not tested yet ***\n");
863         goto fath;
864      }
865      /* solve LP relaxation of the current subproblem */
866      if (T->parm->msg_lev >= GLP_MSG_DBG)
867         xprintf("Solving LP relaxation...\n");
868      ret = ios_solve_node(T);
869      if (!(ret == 0 || ret == GLP_EOBJLL || ret == GLP_EOBJUL))
870      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
871            xprintf("ios_driver: unable to solve current LP relaxation;"
872               " glp_simplex returned %d\n", ret);
873         ret = GLP_EFAIL;
874         goto done;
875      }
876      /* analyze status of the basic solution to LP relaxation found */
877      p_stat = T->mip->pbs_stat;
878      d_stat = T->mip->dbs_stat;
879      if (p_stat == GLP_FEAS && d_stat == GLP_FEAS)
880      {  /* LP relaxation has optimal solution */
881         if (T->parm->msg_lev >= GLP_MSG_DBG)
882            xprintf("Found optimal solution to LP relaxation\n");
883      }
884      else if (d_stat == GLP_NOFEAS)
885      {  /* LP relaxation has no dual feasible solution */
886         /* since the current subproblem cannot have a larger feasible
887            region than its parent, there is something wrong */
888         if (T->parm->msg_lev >= GLP_MSG_ERR)
889            xprintf("ios_driver: current LP relaxation has no dual feas"
890               "ible solution\n");
891         ret = GLP_EFAIL;
892         goto done;
893      }
894      else if (p_stat == GLP_INFEAS && d_stat == GLP_FEAS)
895      {  /* LP relaxation has no primal solution which is better than
896            the incumbent objective value */
897         xassert(T->mip->mip_stat == GLP_FEAS);
898         if (T->parm->msg_lev >= GLP_MSG_DBG)
899            xprintf("LP relaxation has no solution better than incumben"
900               "t objective value\n");
901         /* prune the branch */
902         goto fath;
903      }
904      else if (p_stat == GLP_NOFEAS)
905      {  /* LP relaxation has no primal feasible solution */
906         if (T->parm->msg_lev >= GLP_MSG_DBG)
907            xprintf("LP relaxation has no feasible solution\n");
908         /* prune the branch */
909         goto fath;
910      }
911      else
912      {  /* other cases cannot appear */
913         xassert(T->mip != T->mip);
914      }
915      /* at this point basic solution to LP relaxation of the current
916         subproblem is optimal */
917      xassert(p_stat == GLP_FEAS && d_stat == GLP_FEAS);
918      xassert(T->curr != NULL);
919      T->curr->lp_obj = T->mip->obj_val;
920      /* thus, it defines a local bound to integer optimal solution of
921         the current subproblem */
922      {  double bound = T->mip->obj_val;
923         /* some local bound to the current subproblem could be already
924            set before, so we should only improve it */
925         bound = ios_round_bound(T, bound);
926         if (T->mip->dir == GLP_MIN)
927         {  if (T->curr->bound < bound)
928               T->curr->bound = bound;
929         }
930         else if (T->mip->dir == GLP_MAX)
931         {  if (T->curr->bound > bound)
932               T->curr->bound = bound;
933         }
934         else
935            xassert(T->mip != T->mip);
936         if (T->parm->msg_lev >= GLP_MSG_DBG)
937            xprintf("Local bound is %.9e\n", bound);
938      }
939      /* if the local bound indicates that integer optimal solution of
940         the current subproblem cannot be better than the global bound,
941         prune the branch */
942      if (!is_branch_hopeful(T, p))
943      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
944            xprintf("Current branch is hopeless and can be pruned\n");
945         goto fath;
946      }
947      /* let the application program generate additional rows ("lazy"
948         constraints) */
949      xassert(T->reopt == 0);
950      xassert(T->reinv == 0);
951      if (T->parm->cb_func != NULL)
952      {  xassert(T->reason == 0);
953         T->reason = GLP_IROWGEN;
954         T->parm->cb_func(T, T->parm->cb_info);
955         T->reason = 0;
956         if (T->stop)
957         {  ret = GLP_ESTOP;
958            goto done;
959         }
960         if (T->reopt)
961         {  /* some rows were added; re-optimization is needed */
962            T->reopt = T->reinv = 0;
963            goto more;
964         }
965         if (T->reinv)
966         {  /* no rows were added, however, some inactive rows were
967               removed */
968            T->reinv = 0;
969            xassert(glp_factorize(T->mip) == 0);
970         }
971      }
972      /* check if the basic solution is integer feasible */
973      check_integrality(T);
974      /* if the basic solution satisfies to all integrality conditions,
975         it is a new, better integer feasible solution */
976      if (T->curr->ii_cnt == 0)
977      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
978            xprintf("New integer feasible solution found\n");
979         if (T->parm->msg_lev >= GLP_MSG_ALL)
980            display_cut_info(T);
981         record_solution(T);
982         if (T->parm->msg_lev >= GLP_MSG_ON)
983            show_progress(T, 1);
984         /* make the application program happy */
985         if (T->parm->cb_func != NULL)
986         {  xassert(T->reason == 0);
987            T->reason = GLP_IBINGO;
988            T->parm->cb_func(T, T->parm->cb_info);
989            T->reason = 0;
990            if (T->stop)
991            {  ret = GLP_ESTOP;
992               goto done;
993            }
994         }
995         /* since the current subproblem has been fathomed, prune its
996            branch */
997         goto fath;
998      }
999      /* at this point basic solution to LP relaxation of the current
1000         subproblem is optimal, but integer infeasible */
1001      /* try to fix some non-basic structural variables of integer kind
1002         on their current bounds due to reduced costs */
1003      if (T->mip->mip_stat == GLP_FEAS)
1004         fix_by_red_cost(T);
1005      /* let the application program try to find some solution to the
1006         original MIP with a primal heuristic */
1007      if (T->parm->cb_func != NULL)
1008      {  xassert(T->reason == 0);
1009         T->reason = GLP_IHEUR;
1010         T->parm->cb_func(T, T->parm->cb_info);
1011         T->reason = 0;
1012         if (T->stop)
1013         {  ret = GLP_ESTOP;
1014            goto done;
1015         }
1016         /* check if the current branch became hopeless */
1017         if (!is_branch_hopeful(T, p))
1018         {  if (T->parm->msg_lev >= GLP_MSG_DBG)
1019               xprintf("Current branch became hopeless and can be prune"
1020                  "d\n");
1021            goto fath;
1022         }
1023      }
1024      /* try to find solution with the feasibility pump heuristic */
1025      if (T->parm->fp_heur)
1026      {  xassert(T->reason == 0);
1027         T->reason = GLP_IHEUR;
1028         ios_feas_pump(T);
1029         T->reason = 0;
1030         /* check if the current branch became hopeless */
1031         if (!is_branch_hopeful(T, p))
1032         {  if (T->parm->msg_lev >= GLP_MSG_DBG)
1033               xprintf("Current branch became hopeless and can be prune"
1034                  "d\n");
1035            goto fath;
1036         }
1037      }
1038      /* it's time to generate cutting planes */
1039      xassert(T->local != NULL);
1040      xassert(T->local->size == 0);
1041      /* let the application program generate some cuts; note that it
1042         can add cuts either to the local cut pool or directly to the
1043         current subproblem */
1044      if (T->parm->cb_func != NULL)
1045      {  xassert(T->reason == 0);
1046         T->reason = GLP_ICUTGEN;
1047         T->parm->cb_func(T, T->parm->cb_info);
1048         T->reason = 0;
1049         if (T->stop)
1050         {  ret = GLP_ESTOP;
1051            goto done;
1052         }
1053      }
1054      /* try to generate generic cuts with built-in generators
1055         (as suggested by Matteo Fischetti et al. the built-in cuts
1056         are not generated at each branching node; an intense attempt
1057         of generating new cuts is only made at the root node, and then
1058         a moderate effort is spent after each backtracking step) */
1059      if (T->curr->level == 0 || pred_p == 0)
1060      {  xassert(T->reason == 0);
1061         T->reason = GLP_ICUTGEN;
1062         generate_cuts(T);
1063         T->reason = 0;
1064      }
1065      /* if the local cut pool is not empty, select useful cuts and add
1066         them to the current subproblem */
1067      if (T->local->size > 0)
1068      {  xassert(T->reason == 0);
1069         T->reason = GLP_ICUTGEN;
1070         ios_process_cuts(T);
1071         T->reason = 0;
1072      }
1073      /* clear the local cut pool */
1074      ios_clear_pool(T, T->local);
1075      /* perform re-optimization, if necessary */
1076      if (T->reopt)
1077      {  T->reopt = 0;
1078         T->curr->changed++;
1079         goto more;
1080      }
1081      /* no cuts were generated; remove inactive cuts */
1082      remove_cuts(T);
1083      if (T->parm->msg_lev >= GLP_MSG_ALL && T->curr->level == 0)
1084         display_cut_info(T);
1085      /* update history information used on pseudocost branching */
1086      if (T->pcost != NULL) ios_pcost_update(T);
1087      /* it's time to perform branching */
1088      xassert(T->br_var == 0);
1089      xassert(T->br_sel == 0);
1090      /* let the application program choose variable to branch on */
1091      if (T->parm->cb_func != NULL)
1092      {  xassert(T->reason == 0);
1093         xassert(T->br_var == 0);
1094         xassert(T->br_sel == 0);
1095         T->reason = GLP_IBRANCH;
1096         T->parm->cb_func(T, T->parm->cb_info);
1097         T->reason = 0;
1098         if (T->stop)
1099         {  ret = GLP_ESTOP;
1100            goto done;
1101         }
1102      }
1103      /* if nothing has been chosen, choose some variable as specified
1104         by the branching technique option */
1105      if (T->br_var == 0)
1106         T->br_var = ios_choose_var(T, &T->br_sel);
1107      /* perform actual branching */
1108      curr_p = T->curr->p;
1109      ret = branch_on(T, T->br_var, T->br_sel);
1110      T->br_var = T->br_sel = 0;
1111      if (ret == 0)
1112      {  /* both branches have been created */
1113         pred_p = curr_p;
1114         goto loop;
1115      }
1116      else if (ret == 1)
1117      {  /* one branch is hopeless and has been pruned, so now the
1118            current subproblem is other branch */
1119         /* the current subproblem should be considered as a new one,
1120            since one bound of the branching variable was changed */
1121         T->curr->solved = T->curr->changed = 0;
1122         goto more;
1123      }
1124      else if (ret == 2)
1125      {  /* both branches are hopeless and have been pruned; new
1126            subproblem selection is needed to continue the search */
1127         goto fath;
1128      }
1129      else
1130         xassert(ret != ret);
1131fath: /* the current subproblem has been fathomed */
1132      if (T->parm->msg_lev >= GLP_MSG_DBG)
1133         xprintf("Node %d fathomed\n", p);
1134      /* freeze the current subproblem */
1135      ios_freeze_node(T);
1136      /* and prune the corresponding branch of the tree */
1137      ios_delete_node(T, p);
1138      /* if a new integer feasible solution has just been found, other
1139         branches may become hopeless and therefore must be pruned */
1140      if (T->mip->mip_stat == GLP_FEAS) cleanup_the_tree(T);
1141      /* new subproblem selection is needed due to backtracking */
1142      pred_p = 0;
1143      goto loop;
1144done: /* display progress of the search on exit from the solver */
1145      if (T->parm->msg_lev >= GLP_MSG_ON)
1146         show_progress(T, 0);
1147      if (T->mir_gen != NULL)
1148         ios_mir_term(T->mir_gen), T->mir_gen = NULL;
1149      if (T->clq_gen != NULL)
1150         ios_clq_term(T->clq_gen), T->clq_gen = NULL;
1151      /* return to the calling program */
1152      return ret;
1153}
1154
1155/* eof */
Note: See TracBrowser for help on using the repository browser.