Scippy

SCIP

Solving Constraint Integer Programs

prop_obbt.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file prop_obbt.c
17  * @ingroup PROPAGATORS
18  * @brief optimization-based bound tightening propagator
19  * @author Stefan Weltge
20  */
21 
22 /**@todo if bound tightenings of other propagators are the reason for lpsolstat != SCIP_LPSOLSTAT_OPTIMAL, resolve LP */
23 /**@todo only run more than once in root node if primal bound improved or many cuts were added to the LP */
24 /**@todo filter bounds of a variable already if SCIPisLbBetter()/SCIPisUbBetter() would return FALSE */
25 /**@todo improve warmstarting of LP solving */
26 /**@todo include bound value (finite/infinite) into getScore() function */
27 /**@todo use unbounded ray in filtering */
28 /**@todo do we want to run if the LP is unbounded, maybe for infinite variable bounds? */
29 /**@todo add first filter round in direction of objective function */
30 /**@todo implement conflict resolving callback by calling public method of genvbounds propagator, since the reason are
31  * exactly the variable bounds with nonnegative reduced costs stored in the right-hand side of the generated
32  * generalized variable bound (however, this only makes sense if we run locally)
33  */
34 
35 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
36 
37 #include <assert.h>
38 #include <string.h>
39 
40 #include "scip/prop_obbt.h"
41 #include "scip/prop_genvbounds.h"
42 #include "scip/debug.h"
43 
44 #define PROP_NAME "obbt"
45 #define PROP_DESC "optimization-based bound tightening propagator"
46 #define PROP_TIMING SCIP_PROPTIMING_AFTERLPLOOP
47 #define PROP_PRIORITY -1000000 /**< propagator priority */
48 #define PROP_FREQ 0 /**< propagator frequency */
49 #define PROP_DELAY TRUE /**< should propagation method be delayed, if other propagators
50  * found reductions? */
51 #define DEFAULT_CREATE_GENVBOUNDS TRUE /**< should obbt try to provide genvbounds if possible? */
52 #define DEFAULT_FILTERING_NORM TRUE /**< should coefficients in filtering be normalized w.r.t. the
53  * domains sizes? */
54 #define DEFAULT_APPLY_FILTERROUNDS FALSE /**< try to filter bounds in so-called filter rounds by solving
55  * auxiliary LPs? */
56 #define DEFAULT_DUALFEASTOL 1e-9 /**< feasibility tolerance for reduced costs used in obbt; this value
57  * is used if SCIP's dual feastol is greater */
58 #define DEFAULT_CONDITIONLIMIT -1.0 /**< maximum condition limit used in LP solver (-1.0: no limit) */
59 #define DEFAULT_FILTERING_MIN 2 /**< minimal number of filtered bounds to apply another filter
60  * round */
61 #define DEFAULT_ITLIMITFACTOR 5.0 /**< multiple of root node LP iterations used as total LP iteration
62  * limit for obbt (<= 0: no limit ) */
63 #define DEFAULT_MAXLOOKAHEAD 3 /**< maximal number of bounds evaluated without success per group
64  * (-1: no limit) */
65 
66 #define OBBT_SCOREBASE 5 /**< base that is used to calculate a bounds score value */
67 #define GENVBOUND_PROP_NAME "genvbounds"
68 
69 
70 /*
71  * Data structures
72  */
73 
74 /** bound data */
75 struct Bound
76 {
77  SCIP_VAR* var; /**< variable */
78  SCIP_Real newval; /**< stores a probably tighter value for this bound */
79  SCIP_BOUNDTYPE boundtype; /**< type of bound */
80  unsigned int score; /**< score value that is used to group bounds */
81  unsigned int filtered:1; /**< thrown out during pre-filtering step */
82  unsigned int found:1; /**< stores whether a probably tighter value for this bound was found */
83 };
84 typedef struct Bound BOUND;
85 
86 /** bound group */
87 struct BoundGroup
88 {
89  int firstbdindex; /**< index of the first bound of this group in propdata->bounds array */
90  int nbounds; /**< number of bounds in this group */
91 };
92 typedef struct BoundGroup BOUNDGROUP;
93 
94 /** propagator data */
95 struct SCIP_PropData
96 {
97  BOUND** bounds; /**< array of interesting bounds */
98  BOUNDGROUP* boundgroups; /**< array of bound groups */
99  SCIP_ROW* cutoffrow; /**< pointer to current objective cutoff row */
100  SCIP_PROP* genvboundprop; /**< pointer to genvbound propagator */
101  SCIP_Longint lastnode; /**< number of last node where obbt was performed */
102  SCIP_Real dualfeastol; /**< feasibility tolerance for reduced costs used in obbt; this value is
103  * used if SCIP's dual feastol is greater */
104  SCIP_Real conditionlimit; /**< maximum condition limit used in LP solver (-1.0: no limit) */
105  SCIP_Real itlimitfactor; /**< LP iteration limit for obbt will be this factor times total LP
106  * iterations in root node */
107  SCIP_Bool applyfilterrounds; /**< apply filter rounds? */
108  SCIP_Bool creategenvbounds; /**< should obbt try to provide genvbounds if possible? */
109  SCIP_Bool normalize; /**< should coefficients in filtering be normalized w.r.t. the domains
110  * sizes? */
111  int maxlookahead; /**< maximal number of bounds evaluated without success per group
112  * (-1: no limit) */
113  int nbounds; /**< length of interesting bounds array */
114  int nboundgroups; /**< length of boundgroups array */
115  int nminfilter; /**< minimal number of filtered bounds to apply another filter round */
116 };
117 
118 
119 /*
120  * Local methods
121  */
122 
123 /** solves the LP and handles errors */
124 static
125 SCIP_RETCODE solveLP(
126  SCIP* scip, /**< SCIP data structure */
127  int itlimit, /**< maximal number of LP iterations to perform, or -1 for no limit */
128  SCIP_Bool* error, /**< pointer to store whether an unresolved LP error occurred */
129  SCIP_Bool* optimal /**< was the LP solved to optimalilty? */
130  )
131 {
132  SCIP_LPSOLSTAT lpsolstat;
133  SCIP_RETCODE retcode;
134 
135  assert(scip != NULL);
136  assert(itlimit == -1 || itlimit >= 0);
137  assert(error != NULL);
138  assert(optimal != NULL);
139 
140  *optimal = FALSE;
141  *error = FALSE;
142 
143  retcode = SCIPsolveDiveLP(scip, itlimit, error, NULL);
144  lpsolstat = SCIPgetLPSolstat(scip);
145 
146  /* an error should not kill the overall solving process */
147  if( retcode != SCIP_OKAY )
148  {
149  SCIPwarningMessage(scip, " error while solving LP in obbt propagator; LP solve terminated with code <%d>\n", retcode);
150  SCIPwarningMessage(scip, " this does not affect the remaining solution procedure --> continue\n");
151 
152  *error = TRUE;
153 
154  return SCIP_OKAY;
155  }
156 
157  if( lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
158  {
159  assert(!*error);
160  *optimal = TRUE;
161  }
162 #ifdef SCIP_DEBUG
163  else
164  {
165  switch( lpsolstat )
166  {
167  case SCIP_LPSOLSTAT_ITERLIMIT:
168  SCIPdebugMessage(" reached lp iteration limit\n");
169  break;
170  case SCIP_LPSOLSTAT_TIMELIMIT:
171  SCIPdebugMessage(" reached time limit while solving lp\n");
172  break;
173  case SCIP_LPSOLSTAT_UNBOUNDEDRAY:
174  SCIPdebugMessage(" lp was unbounded\n");
175  break;
176  case SCIP_LPSOLSTAT_NOTSOLVED:
177  SCIPdebugMessage(" lp was not solved\n");
178  break;
179  case SCIP_LPSOLSTAT_ERROR:
180  SCIPdebugMessage(" an error occured during solving lp\n");
181  break;
182  case SCIP_LPSOLSTAT_INFEASIBLE:
183  case SCIP_LPSOLSTAT_OBJLIMIT:
184  case SCIP_LPSOLSTAT_OPTIMAL: /* should not appear because it is handled earlier */
185  default:
186  SCIPdebugMessage(" received an unexpected solstat during solving lp: %d\n", lpsolstat);
187  }
188  }
189 #endif
190 
191  return SCIP_OKAY;
192 }
193 
194 /** adds the objective cutoff to the LP; must be in diving mode */
195 static
196 SCIP_RETCODE addObjCutoff(
197  SCIP* scip, /**< SCIP data structure */
198  SCIP_PROPDATA* propdata /**< data of the obbt propagator */
199  )
200 {
201  SCIP_ROW* row;
202  SCIP_VAR** vars;
203  char rowname[SCIP_MAXSTRLEN];
204 
205  int nvars;
206  int i;
207 
208  assert(scip != NULL);
209  assert(SCIPinDive(scip));
210  assert(propdata != NULL);
211  assert(propdata->cutoffrow == NULL);
212 
213  if( SCIPisInfinity(scip, SCIPgetCutoffbound(scip)) )
214  {
215  SCIPdebugMessage("no objective cutoff since there is no cutoff bound\n");
216  return SCIP_OKAY;
217  }
218 
219  SCIPdebugMessage("create objective cutoff and add it to the LP\n");
220 
221  /* get variables data */
222  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
223 
224  /* create objective cutoff row; set local flag to FALSE since primal cutoff is globally valid */
225  (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "obbt_objcutoff");
226  SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, rowname, -SCIPinfinity(scip), SCIPgetCutoffbound(scip), FALSE, FALSE, FALSE) );
227  SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
228 
229  for( i = 0; i < nvars; i++ )
230  {
231  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i], SCIPvarGetObj(vars[i])) );
232  }
233  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
234 
235  /* add row to the LP */
236  SCIP_CALL( SCIPaddRowDive(scip, row) );
237 
238  propdata->cutoffrow = row;
239  assert(SCIProwIsInLP(propdata->cutoffrow));
240 
241  return SCIP_OKAY;
242 }
243 
244 /** determines, whether a variable is already locally fixed */
245 static
246 SCIP_Bool varIsFixedLocal(
247  SCIP* scip, /**< SCIP data structure */
248  SCIP_VAR* var /**< variable to check */
249  )
250 {
251  return SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
252 }
253 
254 /** returns number of LP iterations left (-1: no limit ) */
255 static
256 int getIterationsLeft(
257  SCIP* scip, /**< SCIP data structure */
258  SCIP_Longint nolditerations, /**< iterations count at the beginning of the corresponding function */
259  SCIP_Longint itlimit /**< LP iteration limit (-1: no limit) */
260  )
261 {
262  SCIP_Longint itsleft;
263 
264  assert(scip != NULL);
265  assert(nolditerations >= 0);
266  assert(itlimit == -1 || itlimit >= 0);
267 
268  if( itlimit == -1 )
269  {
270  SCIPdebugMessage("iterations left: unlimited\n");
271  return -1;
272  }
273  else
274  {
275  itsleft = itlimit - ( SCIPgetNLPIterations(scip) - nolditerations );
276  itsleft = MAX(itsleft, 0);
277  itsleft = MIN(itsleft, INT_MAX);
278 
279  SCIPdebugMessage("iterations left: %d\n", (int) itsleft);
280  return (int) itsleft;
281  }
282 }
283 
284 /** returns the objective coefficient for a variable's bound that will be chosen during filtering */
285 static
286 SCIP_Real getFilterCoef(
287  SCIP* scip, /**< SCIP data structure */
288  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
289  SCIP_VAR* var, /**< variable */
290  SCIP_BOUNDTYPE boundtype /**< boundtype to be filtered? */
291  )
292 {
293  SCIP_Real lb;
294  SCIP_Real ub;
295 
296  assert(scip != NULL);
297  assert(propdata != NULL);
298  assert(var != NULL);
299 
300  lb = SCIPvarGetLbLocal(var);
301  ub = SCIPvarGetUbLocal(var);
302 
303  /* this function should not be called for fixed variables */
304  assert(!varIsFixedLocal(scip, var));
305 
306  /* infinite bounds will not be reached */
307  if( boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisInfinity(scip, -lb) )
308  return 0.0;
309  if( boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisInfinity(scip, ub) )
310  return 0.0;
311 
312  if( propdata->normalize )
313  {
314  /* if the length of the domain is too large then the coefficient should be set to +/- 1.0 */
315  if( boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisInfinity(scip, ub) )
316  return 1.0;
317  if( boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisInfinity(scip, -lb) )
318  return -1.0;
319 
320  /* otherwise the coefficient is +/- 1.0 / ( ub - lb ) */
321  return boundtype == SCIP_BOUNDTYPE_LOWER ? 1.0 / (ub - lb) : -1.0 / (ub - lb);
322  }
323  else
324  {
325  return boundtype == SCIP_BOUNDTYPE_LOWER ? 1.0 : -1.0;
326  }
327 }
328 
329 /** trying to filter some bounds using the existing LP solution */
330 static
331 SCIP_RETCODE filterExistingLP(
332  SCIP* scip, /**< original SCIP data structure */
333  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
334  int* nfiltered /**< how many bounds were filtered this round? */
335  )
336 {
337  int i;
338 
339  assert(scip != NULL);
340  assert(propdata != NULL);
341  assert(nfiltered != NULL);
342 
343  *nfiltered = 0;
344 
345  /* only apply filtering if an LP solution is at hand */
346  if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
347  {
348  SCIPdebugMessage("can't filter using existing lp solution since it was not solved to optimality\n");
349  return SCIP_OKAY;
350  }
351 
352  /* check if a bound is tight */
353  for( i = 0; i < propdata->nbounds; i++ )
354  {
355  BOUND* bound; /* shortcut for current bound */
356 
357  SCIP_Real solval; /* the variables value in the current solution */
358  SCIP_Real boundval; /* current local bound for the variable */
359 
360  bound = propdata->bounds[i];
361  if( bound->filtered )
362  continue;
363 
364  boundval = bound->boundtype == SCIP_BOUNDTYPE_UPPER ?
365  SCIPvarGetUbLocal(bound->var) : SCIPvarGetLbLocal(bound->var);
366  solval = SCIPvarGetLPSol(bound->var);
367 
368  /* bound is tight; since this holds for all fixed variables, those are filtered here automatically */
369  if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisFeasGE(scip, solval, boundval))
370  || (bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisFeasLE(scip, solval, boundval)) )
371  {
372  /* mark bound as filtered */
373  bound->filtered = TRUE;
374 
375  /* increase number of filtered variables */
376  (*nfiltered)++;
377  }
378  }
379 
380  return SCIP_OKAY;
381 }
382 
383 /** enforces one round of filtering */
384 static
385 SCIP_RETCODE filterRound(
386  SCIP* scip, /**< SCIP data structure */
387  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
388  int itlimit, /**< LP iteration limit (-1: no limit) */
389  int* nfiltered /**< how many bounds were filtered this round */
390  )
391 {
392  SCIP_VAR** vars; /* array of the problems variables */
393  SCIP_Bool error;
394  SCIP_Bool optimal;
395 
396  int nvars; /* number of the problems variables */
397  int i;
398 
399  assert(scip != NULL);
400  assert(SCIPinDive(scip));
401  assert(propdata != NULL);
402  assert(itlimit == -1 || itlimit >= 0);
403  assert(nfiltered != NULL);
404 
405  *nfiltered = 0;
406 
407  /* get variable data */
408  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
409 
410  /* solve LP */
411  SCIP_CALL( solveLP(scip, itlimit, &error, &optimal) );
412 
413  if( !optimal )
414  {
415  SCIPdebugMessage("skipping filter round since the LP was not solved to optimality\n");
416  return SCIP_OKAY;
417  }
418 
419  assert(!error);
420 
421  /* check if a bound is tight */
422  for( i = 0; i < propdata->nbounds; i++ )
423  {
424  BOUND* bound; /* shortcut for current bound */
425 
426  SCIP_Real solval; /* the variables value in the current solution */
427  SCIP_Real boundval; /* current local bound for the variable */
428 
429  bound = propdata->bounds[i];
430 
431  /* if bound is filtered it was handled already before */
432  if( bound->filtered )
433  continue;
434 
435  boundval = bound->boundtype == SCIP_BOUNDTYPE_UPPER ?
436  SCIPvarGetUbLocal(bound->var) : SCIPvarGetLbLocal(bound->var);
437  solval = SCIPvarGetLPSol(bound->var);
438 
439  /* bound is tight */
440  if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisFeasGE(scip, solval, boundval))
441  || (bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisFeasLE(scip, solval, boundval)) )
442  {
443  SCIP_Real objcoef;
444 
445  /* mark bound as filtered */
446  bound->filtered = TRUE;
447 
448  /* increase number of filtered variables */
449  (*nfiltered)++;
450 
451  /* get the corresponding variable's objective coeffient */
452  objcoef = SCIPgetVarObjDive(scip, bound->var);
453 
454  /* change objective coefficient if it was set up for this bound */
455  if( (bound->boundtype == SCIP_BOUNDTYPE_UPPER && SCIPisNegative(scip, objcoef))
456  || (bound->boundtype == SCIP_BOUNDTYPE_LOWER && SCIPisPositive(scip, objcoef)) )
457  {
458  SCIP_CALL( SCIPchgVarObjDive(scip, bound->var, 0.0) );
459  }
460  }
461  }
462 
463  return SCIP_OKAY;
464 }
465 
466 /** filter some bounds that are not improvable by solving auxiliary LPs */
467 static
468 SCIP_RETCODE filterBounds(
469  SCIP* scip, /**< SCIP data structure */
470  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
471  SCIP_Longint itlimit /**< LP iteration limit (-1: no limit) */
472  )
473 {
474  SCIP_VAR** vars;
475  SCIP_Longint nolditerations;
476  int nleftiterations;
477  int i;
478  int nfiltered;
479  int nvars;
480 
481  assert(scip != NULL);
482  assert(SCIPinDive(scip));
483  assert(propdata != NULL);
484  assert(itlimit == -1 || itlimit >= 0);
485 
486  nolditerations = SCIPgetNLPIterations(scip);
487  nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
488 
489  /* get variable data */
490  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
491 
492  SCIPdebugMessage("start filter rounds\n");
493 
494  /*
495  * 1.) Try first to filter lower bounds of interesting variables, whose bounds are not already filtered
496  */
497 
498  for( i = 0; i < nvars; i++ )
499  {
500  SCIP_CALL( SCIPchgVarObjDive(scip, vars[i], 0.0) );
501  }
502 
503  for( i = 0; i < propdata->nbounds; i++ )
504  {
505  if( propdata->bounds[i]->boundtype == SCIP_BOUNDTYPE_LOWER && !propdata->bounds[i]->filtered )
506  {
507  SCIP_CALL( SCIPchgVarObjDive(scip, propdata->bounds[i]->var,
508  getFilterCoef(scip, propdata, propdata->bounds[i]->var, SCIP_BOUNDTYPE_LOWER)) );
509  }
510  }
511 
512  do
513  {
514  SCIPdebugMessage("doing a lower bounds round\n");
515  SCIP_CALL( filterRound(scip, propdata, nleftiterations, &nfiltered) );
516  SCIPdebugMessage("filtered %d more bounds in lower bounds round\n", nfiltered);
517 
518  /* update iterations left */
519  nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
520  }
521  while( nfiltered >= propdata->nminfilter && ( nleftiterations == -1 || nleftiterations > 0 ) );
522 
523  /*
524  * 2.) Now try to filter the remaining upper bounds of interesting variables, whose bounds are not already filtered
525  */
526 
527  for( i = 0; i < nvars; i++ )
528  {
529  SCIP_CALL( SCIPchgVarObjDive(scip, vars[i], 0.0) );
530  }
531 
532  for( i = 0; i < propdata->nbounds; i++ )
533  {
534  if( propdata->bounds[i]->boundtype == SCIP_BOUNDTYPE_UPPER && !propdata->bounds[i]->filtered )
535  {
536  SCIP_CALL( SCIPchgVarObjDive(scip, propdata->bounds[i]->var,
537  getFilterCoef(scip, propdata, propdata->bounds[i]->var, SCIP_BOUNDTYPE_UPPER)) );
538  }
539  }
540 
541  do
542  {
543  SCIPdebugMessage("doing an upper bounds round\n");
544  SCIP_CALL( filterRound(scip, propdata, nleftiterations, &nfiltered) );
545  SCIPdebugMessage("filtered %d more bounds in upper bounds round\n", nfiltered);
546 
547  /* update iterations left */
548  nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
549  }
550  while( nfiltered >= propdata->nminfilter && ( nleftiterations == -1 || nleftiterations > 0 ) );
551 
552  return SCIP_OKAY;
553 }
554 
555 /** applies possible bound changes that were found */
556 static
557 SCIP_RETCODE applyBoundChgs(
558  SCIP* scip, /**< SCIP data structure */
559  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
560  SCIP_RESULT* result /**< result pointer */
561  )
562 {
563 #ifdef SCIP_DEBUG
564  int ntightened; /* stores the number of successful bound changes */
565 #endif
566  int i;
567 
568  assert(scip != NULL);
569  assert(!SCIPinDive(scip));
570  assert(propdata != NULL);
571  assert(result != NULL);
572  assert(*result == SCIP_DIDNOTFIND);
573 
574  SCIPdebug( ntightened = 0 );
575 
576  for( i = 0; i < propdata->nbounds; i++ )
577  {
578  BOUND* bound; /* shortcut to the current bound */
579  SCIP_Bool infeas; /* stores wether a tightening approach forced an infeasibilty */
580  SCIP_Bool tightened; /* stores wether a tightening approach was successful */
581 
582  bound = propdata->bounds[i];
583 
584  if( bound->found )
585  {
586  if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
587  {
588  SCIP_CALL( SCIPtightenVarLb(scip, bound->var, bound->newval, FALSE, &infeas, &tightened) );
589  }
590  else
591  {
592  SCIP_CALL( SCIPtightenVarUb(scip, bound->var, bound->newval, FALSE, &infeas, &tightened) );
593  }
594 
595  /* handle information about the success */
596  if( infeas )
597  {
598  *result = SCIP_CUTOFF;
599  SCIPdebugMessage("cut off\n");
600  break;
601  }
602 
603  if( tightened )
604  {
605  *result = SCIP_REDUCEDDOM;
606  SCIPdebug( ntightened++ );
607  }
608  }
609  }
610 
611 #ifdef SCIP_DEBUG
612  SCIPdebugMessage("tightened %d bounds\n", ntightened);
613 #endif
614 
615  return SCIP_OKAY;
616 }
617 
618 /** tries to tighten a bound in diving mode */
619 static
620 SCIP_RETCODE tightenBoundDive(
621  SCIP* scip, /**< SCIP data structure */
622  BOUND* bound, /**< bound that could be tightened */
623  SCIP_Real newval, /**< new bound value */
624  SCIP_Bool* tightened /**< was tightening successful? */
625  )
626 {
627  SCIP_Real lb;
628  SCIP_Real ub;
629 
630  assert(scip != NULL);
631  assert(SCIPinDive(scip));
632  assert(bound != NULL);
633  assert(tightened != NULL);
634 
635  *tightened = FALSE;
636 
637  /* get old bounds */
638  lb = SCIPgetVarLbDive(scip, bound->var);
639  ub = SCIPgetVarUbDive(scip, bound->var);
640 
641  if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
642  {
643  /* round bounds new value if variable is integral */
644  if( SCIPvarIsIntegral(bound->var) )
645  newval = SCIPceil(scip, newval);
646 
647  /* ensure that we give consistent bounds to the LP solver */
648  if( newval > ub )
649  newval = ub;
650 
651  /* tighten if really better */
652  if( SCIPisLbBetter(scip, newval, lb, ub) )
653  {
654  SCIP_CALL( SCIPchgVarLbDive(scip, bound->var, newval) );
655  *tightened = TRUE;
656  }
657  }
658  else
659  {
660  /* round bounds new value if variable is integral */
661  if( SCIPvarIsIntegral(bound->var) )
662  newval = SCIPfloor(scip, newval);
663 
664  /* ensure that we give consistent bounds to the LP solver */
665  if( newval < lb )
666  newval = lb;
667 
668  /* tighten if really better */
669  if( SCIPisUbBetter(scip, newval, lb, ub) )
670  {
671  SCIP_CALL( SCIPchgVarUbDive(scip, bound->var, newval) );
672  *tightened = TRUE;
673  }
674  }
675 
676  return SCIP_OKAY;
677 }
678 
679 /** creates a genvbound if the dual LP solution provides such information
680  *
681  * Consider the problem
682  *
683  * min { +/- x_i : obj * x <= z, lb <= Ax <= ub, l <= x <= u },
684  *
685  * where z is the current cutoff bound. Let (mu, nu, gamma, alpha, beta) >= 0 be the optimal solution of the dual of
686  * problem (P), where the variables correspond to the primal inequalities in the following way:
687  *
688  * Ax >= lb <-> mu
689  * -Ax >= -ub <-> nu
690  * -obj * x >= -z <-> gamma
691  * x >= l <-> alpha
692  * -x >= -u <-> beta
693  *
694  * Fixing these multipliers, by weak duality, we obtain the inequality
695  *
696  * +/- x_i >= lb*mu - ub*nu - z*gamma + l*alpha - u*beta
697  *
698  * that holds for all primal feasible points x with objective value at least z. Setting
699  *
700  * c = lb*mu - ub*nu, redcost_k = alpha_k - beta_k
701  *
702  * we obtain the inequality
703  *
704  * +/- x_i >= sum ( redcost_k * x_k ) + (-gamma) * cutoff_bound + c,
705  *
706  * that holds for all primal feasible points with objective value at least cutoff_bound. Therefore, the latter
707  * inequality can be added as a generalized variable bound.
708  */
709 static
710 SCIP_RETCODE createGenVBound(
711  SCIP* scip, /**< SCIP data structure */
712  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
713  BOUND* bound /**< bound of x_i */
714  )
715 {
716  assert(scip != NULL);
717  assert(bound != NULL);
718  assert(propdata != NULL);
719  assert(propdata->genvboundprop != NULL);
720 
721  /* make sure we are in dive mode having an optimal LP solution */
722  assert(SCIPinDive(scip));
723  assert(SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL);
724 
725  /* only genvbounds created in the root node are globally valid */
726  assert(SCIPgetDepth(scip) == 0);
727 
728  SCIPdebugMessage(" try to create a genvbound for <%s>...\n", SCIPvarGetName(bound->var));
729 
730  /* a genvbound with a multiplier for x_i would not help us */
731  if( SCIPisZero(scip, SCIPgetVarRedcost(scip, bound->var)) )
732  {
733  SCIP_VAR** vars; /* global variables array */
734  SCIP_VAR** genvboundvars; /* genvbound variables array */
735 
736  SCIP_VAR* xi; /* variable x_i */
737 
738  SCIP_Real* genvboundcoefs; /* genvbound coefficients array */
739 
740  SCIP_Real gamma_dual; /* dual multiplier of objective cutoff */
741 
742  int k; /* variable for indexing global variables array */
743  int ncoefs; /* number of nonzero coefficients in genvbound */
744  int nvars; /* number of global variables */
745 
746  /* set x_i */
747  xi = bound->var;
748 
749  /* get variable data */
750  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
751 
752  /* count nonzero coefficients in genvbound */
753  ncoefs = 0;
754  for( k = 0; k < nvars; k++ )
755  {
756  if( !SCIPisZero(scip, SCIPgetVarRedcost(scip, vars[k])) )
757  {
758  assert(vars[k] != xi);
759  ncoefs++;
760  }
761  }
762 
763  /* get dual multiplier for the objective cutoff (set to zero if there is no) */
764  if( propdata->cutoffrow == NULL )
765  {
766  gamma_dual = 0.0;
767  }
768  else
769  {
770  assert(!SCIPisInfinity(scip, SCIPgetCutoffbound(scip)));
771 
772  /* note that the objective cutoff is of the form
773  * -inf <= obj * x <= cutoff_bound
774  * but we want the positive dual multiplier!
775  */
776  gamma_dual = -SCIProwGetDualsol(propdata->cutoffrow);
777  }
778 
779  /* we need at least one nonzero coefficient or a nonzero dual multiplier for the objective cutoff */
780  if( ncoefs > 0 || !SCIPisZero(scip, gamma_dual) )
781  {
782  SCIP_Real c; /* helper variable to calculate constant term in genvbound */
783  int idx; /* variable for indexing genvbound's coefficients array */
784 
785  /* there should be no coefficient for x_i */
786  assert(SCIPisZero(scip, SCIPgetVarRedcost(scip, xi)));
787 
788  /* allocate memory for storing the genvbounds right-hand side variables and coefficients */
789  SCIP_CALL( SCIPallocMemoryArray(scip, &(genvboundvars), ncoefs) );
790  SCIP_CALL( SCIPallocMemoryArray(scip, &(genvboundcoefs), ncoefs) );
791 
792  /* set c = lb*mu - ub*nu - z*gamma + l*alpha - u*beta */
793  c = SCIPgetLPObjval(scip);
794 
795  /* subtract ( - z * gamma ) from c */
796  c += SCIPgetCutoffbound(scip) * gamma_dual;
797 
798  /* subtract ( l*alpha - u*beta ) from c and set the coefficients of the variables */
799  idx = 0;
800  for( k = 0; k < nvars; k++ )
801  {
802  SCIP_VAR* xk;
803  SCIP_Real redcost;
804 
805  xk = vars[k];
806  redcost = SCIPgetVarRedcost(scip, xk);
807 
808  if( !SCIPisZero(scip, redcost) )
809  {
810  assert(xk != xi);
811 
812  /* store coefficients */
813  assert(idx < ncoefs);
814  genvboundvars[idx] = xk;
815  genvboundcoefs[idx] = redcost;
816  idx++;
817 
818  /* if redcost > 0, then redcost = alpha_k, otherwise redcost = - beta_k */
819  assert(redcost <= 0 || !SCIPisInfinity(scip, -SCIPgetVarLbDive(scip, xk)));
820  assert(redcost >= 0 || !SCIPisInfinity(scip, SCIPgetVarUbDive(scip, xk)));
821  c -= redcost > 0 ? redcost * SCIPgetVarLbDive(scip, xk) : redcost * SCIPgetVarUbDive(scip, xk);
822  }
823  }
824 
825  /* add genvbound */
826  SCIPdebugMessage(" adding genvbound\n");
827  SCIP_CALL( SCIPgenVBoundAdd(scip, propdata->genvboundprop, genvboundvars, xi, genvboundcoefs, ncoefs,
828  !SCIPisPositive(scip, gamma_dual) ? 0.0 : -gamma_dual, c, bound->boundtype) );
829 
830  /* free arrays */
831  SCIPfreeMemoryArray(scip, &genvboundcoefs);
832  SCIPfreeMemoryArray(scip, &genvboundvars);
833  }
834  else
835  {
836  SCIPdebugMessage(" trivial genvbound, skipping\n");
837  }
838  }
839  else
840  {
841  SCIPdebugMessage(" found multiplier for <%s>: %g, skipping\n",
842  SCIPvarGetName(bound->var), SCIPgetVarRedcost(scip, bound->var));
843  }
844 
845  return SCIP_OKAY;
846 }
847 
848 /** tries to find tighter values for bounds and stores them in the bound data structure */
849 static
850 SCIP_RETCODE findNewBounds(
851  SCIP* scip, /**< SCIP data structure */
852  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
853  SCIP_Longint itlimit /**< LP iteration limit (-1: no limit) */
854  )
855 {
856  SCIP_VAR** vars; /* array of the problems variables */
857  int* nextindices; /* next relative index to try per bound group */
858 
859  SCIP_Bool boundsleft;
860  SCIP_Bool iterationsleft;
861  SCIP_Longint nolditerations;
862  int nleftiterations;
863  int i;
864  int obbtround;
865  int nvars; /* number of the problems variables */
866 
867  assert(scip != NULL);
868  assert(SCIPinDive(scip));
869  assert(propdata != NULL);
870  assert(itlimit == -1 || itlimit >= 0);
871 
872  SCIPdebugMessage("solve LPs...\n");
873 
874  nolditerations = SCIPgetNLPIterations(scip);
875  nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
876 
877  /* get variable data */
878  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
879 
880  /* set objective coefficients to zero */
881  for( i = 0; i < nvars; i++ )
882  {
883  SCIP_CALL( SCIPchgVarObjDive(scip, vars[i], 0.0) );
884  }
885 
886  SCIP_CALL( SCIPallocBufferArray(scip, &nextindices, propdata->nboundgroups) );
887  BMSclearMemoryArray(nextindices, propdata->nboundgroups);
888 
889  boundsleft = TRUE;
890  iterationsleft = nleftiterations == -1 || nleftiterations > 0;
891 
892  /* as long as bounds are left and maxlookahead is not exceeded, iterate over all groups */
893  for( obbtround = 1; boundsleft && (propdata->maxlookahead == -1 || obbtround <= propdata->maxlookahead)
894  && iterationsleft; obbtround++ )
895  {
896  SCIPdebugMessage("obbt round %d\n", obbtround);
897 
898  boundsleft = FALSE;
899 
900  for( i = 0; i < propdata->nboundgroups && iterationsleft; i++ )
901  {
902  BOUNDGROUP* bdgroup;
903  SCIP_Bool lastboundsuccessful;
904 
905  bdgroup = &(propdata->boundgroups[i]);
906  lastboundsuccessful = TRUE;
907 
908  while( lastboundsuccessful && nextindices[i] < bdgroup->nbounds && iterationsleft )
909  {
910  BOUND* bound;
911  SCIP_VAR* var;
912 
913  SCIP_Bool error;
914  SCIP_Bool optimal; /* was the LP solved to optimalilty? */
915 
916  bound = propdata->bounds[bdgroup->firstbdindex + nextindices[i]];
917  var = bound->var;
918  nextindices[i]++;
919 
920  SCIPdebugMessage(" applying obbt on %s bound of <%s> (local bounds: [%f,%f])\n",
921  bound->boundtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper", SCIPvarGetName(var),
922  SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
923 
924  /* ignore filtered bounds */
925  if( bound->filtered )
926  {
927  SCIPdebugMessage(" bound already filtered\n");
928  continue;
929  }
930 
931  assert(!varIsFixedLocal(scip, var));
932 
933  lastboundsuccessful = FALSE;
934 
935  /* set objective coefficient to +/- 1 (note that we minimize) */
936  SCIP_CALL( SCIPchgVarObjDive(scip, var, (bound->boundtype == SCIP_BOUNDTYPE_LOWER) ? 1.0 : -1.0 ) );
937 
938  /* solve LP */
939  SCIP_CALL( solveLP(scip, nleftiterations, &error, &optimal) );
940 
941  /* update iterations left */
942  nleftiterations = getIterationsLeft(scip, nolditerations, itlimit);
943  iterationsleft = nleftiterations == -1 || nleftiterations > 0;
944 
945  /* stop this procedure if an error occured */
946  if( error )
947  {
948  SCIPfreeBufferArray(scip, &nextindices);
949  return SCIP_OKAY;
950  }
951 
952  /* this variable should also be ignored in further calls to filterExistingLP() */
953  bound->filtered = TRUE;
954 
955  if( optimal )
956  {
957  int nfiltered;
958 
959  /* try to filter further bounds using the lp solution */
960  nfiltered = 0;
961  SCIP_CALL( filterExistingLP(scip, propdata, &nfiltered));
962  if( nfiltered > 0 )
963  {
964  SCIPdebugMessage(" filtered %d more bounds using existing lp solution\n", nfiltered);
965  }
966 
967  /* store this value in the bound data structure */
968  bound->newval = SCIPvarGetLPSol(var);
969  bound->found = TRUE;
970 
971  SCIPdebugMessage(" var <%s>, LP value: %f\n", SCIPvarGetName(var), bound->newval);
972 
973 #ifdef SCIP_DEBUG_SOLUTION
974  if( bound->boundtype == SCIP_BOUNDTYPE_LOWER )
975  {
976  SCIP_CALL( SCIPdebugCheckLbGlobal(scip, var, bound->newval) );
977  }
978  else
979  {
980  SCIP_CALL( SCIPdebugCheckUbGlobal(scip, var, bound->newval) );
981  }
982 #endif
983 
984  /* in root node we may want to create a genvbound (independent of tightening success) */
985  if( SCIPgetDepth(scip) == 0 && propdata->genvboundprop != NULL )
986  {
987  SCIP_CALL( createGenVBound(scip, propdata, bound) );
988  }
989 
990  /* try to tighten bound in dive mode */
991  SCIP_CALL( tightenBoundDive(scip, bound, bound->newval, &lastboundsuccessful) );
992  }
993 
994  /* set objective coefficient back to zero */
995  SCIP_CALL( SCIPchgVarObjDive(scip, var, 0.0 ) );
996  }
997 
998  if( nextindices[i] < bdgroup->nbounds )
999  boundsleft = TRUE;
1000  }
1001  }
1002 
1003  SCIPfreeBufferArray(scip, &nextindices);
1004  return SCIP_OKAY;
1005 }
1006 
1007 /** main function of obbt */
1008 static
1009 SCIP_RETCODE applyObbt(
1010  SCIP* scip, /**< SCIP data structure */
1011  SCIP_PROPDATA* propdata, /**< data of the obbt propagator */
1012  SCIP_Longint itlimit, /**< LP iteration limit (-1: no limit) */
1013  SCIP_RESULT* result /**< result pointer */
1014  )
1015 {
1016  SCIP_Longint nolditerations;
1017  SCIP_Bool hasconditionlimit;
1018  SCIP_Real oldconditionlimit;
1019  SCIP_Real olddualfeastol;
1020  int nfiltered;
1021  int i;
1022 
1023  assert(scip != NULL);
1024  assert(propdata != NULL);
1025  assert(itlimit == -1 || itlimit >= 0);
1026 
1027  *result = SCIP_DIDNOTFIND;
1028  nolditerations = SCIPgetNLPIterations(scip);
1029  nfiltered = 0;
1030 
1031  /* reset bound data structure flags; fixed variables are marked as filtered */
1032  for( i = 0; i < propdata->nbounds; i++ )
1033  {
1034  propdata->bounds[i]->filtered = varIsFixedLocal(scip, propdata->bounds[i]->var);
1035  propdata->bounds[i]->found = FALSE;
1036  }
1037 
1038  /* filter variables via inspecting present LP solution */
1039  SCIP_CALL( filterExistingLP(scip, propdata, &nfiltered) );
1040  SCIPdebugMessage("filtered %d bounds via inspecting present LP solution\n", nfiltered);
1041 
1042  /* start diving */
1043  SCIP_CALL( SCIPstartDive(scip) );
1044 
1045  /* tighten dual feastol */
1046  olddualfeastol = SCIPdualfeastol(scip);
1047  if( propdata->dualfeastol < olddualfeastol )
1048  {
1049  SCIP_CALL( SCIPchgDualfeastol(scip, propdata->dualfeastol) );
1050  }
1051 
1052  /* tighten condition limit */
1053  hasconditionlimit = (SCIPgetRealParam(scip, "lp/conditionlimit", &oldconditionlimit) == SCIP_OKAY);
1054  if( !hasconditionlimit )
1055  {
1056  SCIPwarningMessage(scip, "obbt propagator could not set condition limit in LP solver - running without\n");
1057  }
1058  else if( propdata->conditionlimit > 0.0 && (oldconditionlimit < 0.0 || propdata->conditionlimit < oldconditionlimit) )
1059  {
1060  SCIP_CALL( SCIPsetRealParam(scip, "lp/conditionlimit", propdata->conditionlimit) );
1061  }
1062 
1063  /* add objective cutoff */
1064  SCIP_CALL( addObjCutoff(scip, propdata) );
1065 
1066  /* apply filtering */
1067  if( propdata->applyfilterrounds )
1068  {
1069  SCIP_CALL( filterBounds(scip, propdata, itlimit) );
1070  }
1071 
1072  /**@todo maybe endDive, startDive, addObjCutoff to restore old LP basis information here */
1073 
1074  if( itlimit > 0 )
1075  {
1076  itlimit = itlimit - ( SCIPgetNLPIterations(scip) - nolditerations );
1077  itlimit = MAX(itlimit, 0);
1078  }
1079 
1080  /* try to find new bounds and store them in the bound data structure */
1081  SCIP_CALL( findNewBounds(scip, propdata, itlimit) );
1082 
1083  /* reset dual feastol and condition limit */
1084  SCIP_CALL( SCIPchgDualfeastol(scip, olddualfeastol) );
1085  if( hasconditionlimit )
1086  {
1087  SCIP_CALL( SCIPsetRealParam(scip, "lp/conditionlimit", oldconditionlimit) );
1088  }
1089 
1090  /* end diving */
1091  SCIP_CALL( SCIPendDive(scip) );
1092 
1093  /* release cutoff row if there is one */
1094  if( propdata->cutoffrow != NULL )
1095  {
1096  assert(!SCIProwIsInLP(propdata->cutoffrow));
1097  SCIP_CALL( SCIPreleaseRow(scip, &(propdata->cutoffrow)) );
1098  }
1099 
1100  /* apply buffered bound changes */
1101  SCIP_CALL( applyBoundChgs(scip, propdata, result) );
1102 
1103  return SCIP_OKAY;
1104 }
1105 
1106 /** computes the score of a bound */
1107 static
1108 unsigned int getScore(
1109  SCIP* scip, /**< SCIP data structure */
1110  BOUND* bound, /**< pointer of bound */
1111  int nlcount, /**< number of nonlinear constraints containing the bounds variable */
1112  int maxnlcount /**< maximal number of nonlinear constraints a variable appears in */
1113  )
1114 {
1115  unsigned int score; /* score to be computed */
1116 
1117  assert(scip != NULL);
1118  assert(bound != NULL);
1119  assert(nlcount >= 0);
1120  assert(maxnlcount >= nlcount);
1121 
1122  /* score = ( nlcount * ( BASE - 1 ) / maxnlcount ) * BASE^2 + vartype * BASE + boundtype */
1123  score = (unsigned int) ( nlcount > 0 ? (OBBT_SCOREBASE * nlcount * ( OBBT_SCOREBASE - 1 )) / maxnlcount : 0 );
1124  switch( SCIPvarGetType(bound->var) )
1125  {
1126  case SCIP_VARTYPE_INTEGER:
1127  score += 1;
1128  break;
1129  case SCIP_VARTYPE_IMPLINT:
1130  score += 2;
1131  break;
1132  case SCIP_VARTYPE_CONTINUOUS:
1133  score += 3;
1134  break;
1135  case SCIP_VARTYPE_BINARY:
1136  score += 4;
1137  break;
1138  default:
1139  break;
1140  }
1141 
1142  score *= OBBT_SCOREBASE;
1143  if( bound->boundtype == SCIP_BOUNDTYPE_UPPER )
1144  score += 1;
1145 
1146  return score;
1147 }
1148 
1149 /** comparison method for two bounds w.r.t. their scores */
1150 static
1151 SCIP_DECL_SORTPTRCOMP(compBounds)
1152 {
1153  BOUND* bound1 = (BOUND*) elem1;
1154  BOUND* bound2 = (BOUND*) elem2;
1155 
1156  return bound1->score == bound2->score ? 0 : ( bound1->score > bound2->score ? 1 : -1 );
1157 }
1158 
1159 #ifdef SCIP_DEBUG
1160 /** prints groups of variables */
1161 static
1162 void printGroups(
1163  SCIP_PROPDATA* propdata /**< data of the obbt propagator */
1164  )
1165 {
1166  int i;
1167 
1168  assert(propdata != NULL);
1169  assert(propdata->nbounds > 0);
1170 
1171  SCIPdebugPrintf("groups={\n");
1172 
1173  for( i = 0; i < propdata->nboundgroups; i++ )
1174  {
1175  int j;
1176 
1177  SCIPdebugPrintf(" {\n");
1178 
1179  for( j = 0; j < propdata->boundgroups[i].nbounds; j++ )
1180  {
1181  BOUND* bound;
1182 
1183  bound = propdata->bounds[propdata->boundgroups[i].firstbdindex + j];
1184 
1185  SCIPdebugPrintf(" %s bound of <%s>, scoreval=%u\n", bound->boundtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper",
1186  SCIPvarGetName(bound->var), bound->score);
1187  }
1188 
1189  SCIPdebugPrintf(" }\n");
1190  }
1191 
1192  SCIPdebugPrintf("}\n");
1193 }
1194 #endif
1195 
1196 /** creates groups for the bounds array */
1197 static
1198 SCIP_RETCODE createGroups(
1199  SCIP_PROPDATA* propdata /**< data of the obbt propagator */
1200  )
1201 {
1202  unsigned int oldscoreval; /* old scores value */
1203 
1204  int i;
1205  int j;
1206 
1207  assert(propdata != NULL);
1208  assert(propdata->nbounds > 0);
1209 
1210  /* count boundgroups */
1211  propdata->nboundgroups = 0;
1212  oldscoreval = propdata->bounds[0]->score + 1;
1213  for( i = 0; i < propdata->nbounds; i++ )
1214  {
1215  if( propdata->bounds[i]->score != oldscoreval )
1216  {
1217  propdata->nboundgroups++;
1218  oldscoreval = propdata->bounds[i]->score;
1219  }
1220  }
1221 
1222  /* allocate bound groups array */
1223  SCIP_CALL( SCIPallocMemoryArray(scip, &(propdata->boundgroups), propdata->nboundgroups) );
1224 
1225  oldscoreval = propdata->bounds[0]->score + 1;
1226  j = 0;
1227  for( i = 0; i < propdata->nbounds; i++ )
1228  {
1229  /* if a new score value appears, create a new group */
1230  if( propdata->bounds[i]->score != oldscoreval )
1231  {
1232  BOUNDGROUP* newgroup;
1233 
1234  newgroup = &(propdata->boundgroups[j]);
1235  newgroup->firstbdindex = i;
1236  newgroup->nbounds = 1;
1237 
1238  j++;
1239 
1240  oldscoreval = propdata->bounds[i]->score;
1241  }
1242  /* otherwise the group has another member */
1243  else
1244  {
1245  assert(j >= 1);
1246  propdata->boundgroups[j - 1].nbounds++;
1247  }
1248  }
1249  assert(j == propdata->nboundgroups);
1250 
1251  return SCIP_OKAY;
1252 }
1253 
1254 /** determines whether a variable is interesting */
1255 static
1256 SCIP_Bool varIsInteresting(
1257  SCIP* scip, /**< SCIP data structure */
1258  SCIP_VAR* var, /**< variable to check */
1259  int nlcount /**< number of nonlinear constraints containing the variable */
1260  )
1261 {
1262  assert(SCIPgetDepth(scip) == 0);
1263 
1264  return !SCIPvarIsBinary(var) && !varIsFixedLocal(scip, var) && nlcount > 0;
1265 }
1266 
1267 /** initializes interesting bounds */
1268 static
1269 SCIP_RETCODE initBounds(
1270  SCIP* scip, /**< SCIP data structure */
1271  SCIP_PROPDATA* propdata /**< data of the obbt propagator */
1272  )
1273 {
1274  SCIP_VAR** vars; /* array of the problems variables */
1275  int* nlcount; /* array that stores in how many nonlinear constraints each variable
1276  * appears */
1277 
1278  int bdidx; /* bound index inside propdata->bounds */
1279  int maxnlcount; /* maximal number of nonlinear constraints a variable appears in */
1280  int nvars; /* number of the problems variables */
1281  int i;
1282 
1283  assert(scip != NULL);
1284  assert(propdata != NULL);
1285 
1286  SCIPdebugMessage("initialize bounds\n");
1287 
1288  /* get variable data */
1289  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
1290 
1291  /* count nonlinearities */
1292  SCIP_CALL( SCIPallocBufferArray(scip, &nlcount, nvars) );
1293  maxnlcount = 0;
1294  if( SCIPisNLPConstructed(scip) )
1295  {
1296  assert(SCIPgetNNLPVars(scip) == nvars);
1297  SCIP_CALL( SCIPgetNLPVarsNonlinearity(scip, nlcount) );
1298 
1299  for( i = 0; i < nvars; i++ )
1300  if( maxnlcount < nlcount[i] )
1301  maxnlcount = nlcount[i];
1302  }
1303  else
1304  {
1305  BMSclearMemoryArray(nlcount, nvars);
1306  }
1307 
1308  /* allocate interesting bounds array */
1309  SCIP_CALL( SCIPallocMemoryArray(scip, &(propdata->bounds), 2 * nvars) );
1310 
1311  /* get all interesting variables and their bounds */
1312  bdidx = 0;
1313  for( i = 0; i < nvars; i++ )
1314  {
1315  if( varIsInteresting(scip, vars[i], nlcount[i]) )
1316  {
1317  BOUND** bdaddress;
1318 
1319  /* create lower bound */
1320  bdaddress = &(propdata->bounds[bdidx]);
1321  SCIP_CALL( SCIPallocMemory(scip, bdaddress) );
1322  propdata->bounds[bdidx]->boundtype = SCIP_BOUNDTYPE_LOWER;
1323  propdata->bounds[bdidx]->var = vars[i];
1324  propdata->bounds[bdidx]->found = FALSE;
1325  propdata->bounds[bdidx]->filtered = FALSE;
1326  propdata->bounds[bdidx]->newval = 0.0;
1327  propdata->bounds[bdidx]->score = getScore(scip, propdata->bounds[bdidx], nlcount[i], maxnlcount);
1328  bdidx++;
1329 
1330  /* create upper bound */
1331  bdaddress = &(propdata->bounds[bdidx]);
1332  SCIP_CALL( SCIPallocMemory(scip, bdaddress) );
1333  propdata->bounds[bdidx]->boundtype = SCIP_BOUNDTYPE_UPPER;
1334  propdata->bounds[bdidx]->var = vars[i];
1335  propdata->bounds[bdidx]->found = FALSE;
1336  propdata->bounds[bdidx]->filtered = FALSE;
1337  propdata->bounds[bdidx]->newval = 0.0;
1338  propdata->bounds[bdidx]->score = getScore(scip, propdata->bounds[bdidx], nlcount[i], maxnlcount);
1339  bdidx++;
1340  }
1341  }
1342 
1343  /* free memory for buffering nonlinearities */
1344  assert(nlcount != NULL);
1345  SCIPfreeBufferArray(scip, &nlcount);
1346 
1347  /* set number of interesting bounds */
1348  propdata->nbounds = bdidx;
1349 
1350  /* resize propdata->bounds array */
1351  if( propdata->nbounds > 0 )
1352  {
1353  SCIP_CALL( SCIPreallocMemoryArray(scip, &(propdata->bounds), propdata->nbounds) );
1354  }
1355  else
1356  {
1357  assert(propdata->nbounds == 0);
1358  SCIPfreeMemoryArray(scip, &(propdata->bounds));
1359  }
1360 
1361  SCIPdebugMessage("problem has %d/%d interesting bounds\n", propdata->nbounds, 2 * nvars);
1362 
1363  if( propdata->nbounds > 0 )
1364  {
1365  /* sort bounds according to decreasing score */
1366  SCIPsortDownPtr((void**) propdata->bounds, compBounds, propdata->nbounds);
1367 
1368  /* create groups */
1369  SCIP_CALL( createGroups(propdata) );
1370  SCIPdebug( printGroups(propdata) );
1371  }
1372 
1373  return SCIP_OKAY;
1374 }
1375 
1376 
1377 /*
1378  * Callback methods of propagator
1379  */
1380 
1381 /** solving process initialization method of propagator (called when branch and bound process is about to begin) */
1382 static
1383 SCIP_DECL_PROPINITSOL(propInitsolObbt)
1384 { /*lint --e{715}*/
1385  SCIP_PROPDATA* propdata;
1386 
1387  assert(scip != NULL);
1388  assert(prop != NULL);
1389  assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
1390 
1391  /* get propagator data */
1392  propdata = SCIPpropGetData(prop);
1393  assert(propdata != NULL);
1394 
1395  propdata->bounds = NULL;
1396  propdata->nbounds = -1;
1397  propdata->boundgroups = NULL;
1398  propdata->nboundgroups = -1;
1399  propdata->cutoffrow = NULL;
1400  propdata->lastnode = -1;
1401 
1402  /* if genvbounds propagator is not available, we cannot create genvbounds */
1403  propdata->genvboundprop = propdata->creategenvbounds ? SCIPfindProp(scip, GENVBOUND_PROP_NAME) : NULL;
1404 
1405  SCIPdebugMessage("creating genvbounds: %s\n", propdata->genvboundprop != NULL ? "true" : "false");
1406 
1407  return SCIP_OKAY;
1408 }
1409 
1410 /** execution method of propagator */
1411 static
1412 SCIP_DECL_PROPEXEC(propExecObbt)
1413 { /*lint --e{715}*/
1414  SCIP_PROPDATA* propdata;
1415 
1416  assert(scip != NULL);
1417  assert(prop != NULL);
1418  assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
1419 
1420  *result = SCIP_DIDNOTRUN;
1421 
1422  /* do not run in: presolving, repropagation, probing mode */
1423  if( SCIPgetStage(scip) != SCIP_STAGE_SOLVING || SCIPinRepropagation(scip) || SCIPinProbing(scip) )
1424  return SCIP_OKAY;
1425 
1426  /* only run for nonlinear problems, i.e., if NLP is constructed */
1427  if( !SCIPisNLPConstructed(scip) )
1428  {
1429  SCIPdebugMessage("NLP not constructed, skipping obbt\n");
1430  return SCIP_OKAY;
1431  }
1432 
1433  /* only run if LP all columns are in the LP, i.e., the LP is a relaxation; e.g., do not run if pricers are active
1434  * since pricing is not performed in diving mode
1435  */
1436  if( !SCIPallColsInLP(scip) )
1437  {
1438  SCIPdebugMessage("not all columns in LP, skipping obbt\n");
1439  return SCIP_OKAY;
1440  }
1441 
1442  /* get propagator data */
1443  propdata = SCIPpropGetData(prop);
1444  assert(propdata != NULL);
1445 
1446  /* ensure that bounds are initialized */
1447  if( propdata->nbounds == -1 )
1448  {
1449  /* bounds must be initialized at root node */
1450  if( SCIPgetDepth(scip) == 0 )
1451  {
1452  SCIP_CALL( initBounds(scip, propdata) );
1453  }
1454  else
1455  {
1456  assert(!SCIPinProbing(scip));
1457  return SCIP_OKAY;
1458  }
1459  }
1460  assert(propdata->nbounds >= 0);
1461 
1462  /* do not run if there are no interesting bounds */
1463  /**@todo disable */
1464  if( propdata->nbounds <= 0 )
1465  {
1466  SCIPdebugMessage("there are no interesting bounds\n");
1467  return SCIP_OKAY;
1468  }
1469 
1470  /* only run once in a node != root */
1471  if( SCIPgetDepth(scip) > 0 && SCIPnodeGetNumber(SCIPgetCurrentNode(scip)) == propdata->lastnode )
1472  {
1473  return SCIP_OKAY;
1474  }
1475 
1476  SCIPdebugMessage("applying obbt for problem <%s> at depth %d\n", SCIPgetProbName(scip), SCIPgetDepth(scip));
1477 
1478  /* without an optimal LP solution we don't want to run; this may be because propagators with higher priority have
1479  * already found reductions or numerical troubles occured during LP solving
1480  */
1481  if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL && SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_UNBOUNDEDRAY )
1482  {
1483  SCIPdebugMessage("aborting since no optimal LP solution is at hand\n");
1484  return SCIP_OKAY;
1485  }
1486 
1487  /* apply obbt */
1488  SCIP_CALL( applyObbt(scip, propdata, propdata->itlimitfactor > 0.0 ?
1489  (SCIP_Longint) (propdata->itlimitfactor * SCIPgetNRootLPIterations(scip)) : -1, result) );
1490  assert(*result != SCIP_DIDNOTRUN);
1491 
1492  /* set current node as last node */
1493  propdata->lastnode = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
1494 
1495  return SCIP_OKAY;
1496 }
1497 
1498 /** propagation conflict resolving method of propagator */
1499 static
1500 SCIP_DECL_PROPRESPROP(propRespropObbt)
1501 { /*lint --e{715}*/
1502  *result = SCIP_DIDNOTFIND;
1503 
1504  return SCIP_OKAY;
1505 }
1506 
1507 /** solving process deinitialization method of propagator (called before branch and bound process data is freed) */
1508 static
1509 SCIP_DECL_PROPEXITSOL(propExitsolObbt)
1510 { /*lint --e{715}*/
1511  SCIP_PROPDATA* propdata;
1512  int i;
1513 
1514  assert(scip != NULL);
1515  assert(prop != NULL);
1516  assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
1517 
1518  /* get propagator data */
1519  propdata = SCIPpropGetData(prop);
1520  assert(propdata != NULL);
1521 
1522  /* free memory allocated for the bounds */
1523  if( propdata->nbounds > 0 )
1524  {
1525  assert(propdata->bounds != NULL);
1526  assert(propdata->boundgroups != NULL);
1527 
1528  /* free bound groups */
1529  SCIPfreeMemoryArray(scip, &(propdata->boundgroups));
1530 
1531  /* free bounds */
1532  for( i = propdata->nbounds - 1; i >= 0; i-- )
1533  {
1534  SCIPfreeMemory(scip, &(propdata->bounds[i]));
1535  }
1536  SCIPfreeMemoryArray(scip, &(propdata->bounds));
1537  }
1538 
1539  propdata->nbounds = -1;
1540  propdata->nboundgroups = -1;
1541 
1542  return SCIP_OKAY;
1543 }
1544 
1545 /** destructor of propagator to free user data (called when SCIP is exiting) */
1546 static
1547 SCIP_DECL_PROPFREE(propFreeObbt)
1548 { /*lint --e{715}*/
1549  SCIP_PROPDATA* propdata;
1550 
1551  assert(strcmp(SCIPpropGetName(prop), PROP_NAME) == 0);
1552 
1553  /* free propagator data */
1554  propdata = SCIPpropGetData(prop);
1555  assert(propdata != NULL);
1556 
1557  SCIPfreeMemory(scip, &propdata);
1558 
1559  SCIPpropSetData(prop, NULL);
1560 
1561  return SCIP_OKAY;
1562 }
1563 
1564 
1565 /*
1566  * propagator specific interface methods
1567  */
1568 
1569 /** creates the obbt propagator and includes it in SCIP */
1570 SCIP_RETCODE SCIPincludePropObbt(
1571  SCIP* scip /**< SCIP data structure */
1572  )
1573 {
1574  SCIP_PROPDATA* propdata;
1575  SCIP_PROP* prop;
1576 
1577  /* create obbt propagator data */
1578  SCIP_CALL( SCIPallocMemory(scip, &propdata) );
1579 
1580  /* include propagator */
1581  SCIP_CALL( SCIPincludePropBasic(scip, &prop, PROP_NAME, PROP_DESC, PROP_PRIORITY, PROP_FREQ, PROP_DELAY, PROP_TIMING,
1582  propExecObbt, propdata) );
1583 
1584  SCIP_CALL( SCIPsetPropFree(scip, prop, propFreeObbt) );
1585  SCIP_CALL( SCIPsetPropExitsol(scip, prop, propExitsolObbt) );
1586  SCIP_CALL( SCIPsetPropInitsol(scip, prop, propInitsolObbt) );
1587  SCIP_CALL( SCIPsetPropResprop(scip, prop, propRespropObbt) );
1588 
1589  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/"PROP_NAME"/creategenvbounds",
1590  "should obbt try to provide genvbounds if possible?",
1591  &propdata->creategenvbounds, TRUE, DEFAULT_CREATE_GENVBOUNDS, NULL, NULL) );
1592 
1593  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/"PROP_NAME"/normalize",
1594  "should coefficients in filtering be normalized w.r.t. the domains sizes?",
1595  &propdata->normalize, TRUE, DEFAULT_FILTERING_NORM, NULL, NULL) );
1596 
1597  SCIP_CALL( SCIPaddBoolParam(scip, "propagating/"PROP_NAME"/applyfilterrounds",
1598  "try to filter bounds in so-called filter rounds by solving auxiliary LPs?",
1599  &propdata->applyfilterrounds, TRUE, DEFAULT_APPLY_FILTERROUNDS, NULL, NULL) );
1600 
1601  SCIP_CALL( SCIPaddIntParam(scip, "propagating/"PROP_NAME"/minfilter",
1602  "minimal number of filtered bounds to apply another filter round",
1603  &propdata->nminfilter, TRUE, DEFAULT_FILTERING_MIN, 1, INT_MAX, NULL, NULL) );
1604 
1605  SCIP_CALL( SCIPaddIntParam(scip, "propagating/"PROP_NAME"/maxlookahead",
1606  "maximal number of bounds evaluated without success per group (-1: no limit)",
1607  &propdata->maxlookahead, FALSE, DEFAULT_MAXLOOKAHEAD, -1, INT_MAX, NULL, NULL) );
1608 
1609  SCIP_CALL( SCIPaddRealParam(scip, "propagating/"PROP_NAME"/itlimitfactor",
1610  "multiple of root node LP iterations used as total LP iteration limit for obbt (<= 0: no limit )",
1611  &propdata->itlimitfactor, FALSE, DEFAULT_ITLIMITFACTOR, SCIP_REAL_MIN, SCIP_REAL_MAX, NULL, NULL) );
1612 
1613  SCIP_CALL( SCIPaddRealParam(scip, "propagating/"PROP_NAME"/dualfeastol",
1614  "feasibility tolerance for reduced costs used in obbt; this value is used if SCIP's dual feastol is greater",
1615  &propdata->dualfeastol, FALSE, DEFAULT_DUALFEASTOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1616 
1617  SCIP_CALL( SCIPaddRealParam(scip, "propagating/"PROP_NAME"/conditionlimit",
1618  "maximum condition limit used in LP solver (-1.0: no limit)",
1619  &propdata->conditionlimit, FALSE, DEFAULT_CONDITIONLIMIT, -1.0, SCIP_REAL_MAX, NULL, NULL) );
1620 
1621  return SCIP_OKAY;
1622 }
1623