source: trunk/Couenne/src/branch/CouenneObject.cpp @ 123

Last change on this file since 123 was 123, checked in by stefan, 11 years ago

use CoinIsnan?

File size: 17.4 KB
Line 
1/*
2 * Name:    CouenneObject.cpp
3 * Authors: Pierre Bonami, IBM Corp.
4 *          Pietro Belotti, Carnegie Mellon University
5 * Purpose: Base object for variables (to be used in branching)
6 *
7 * (C) Carnegie-Mellon University, 2006-08.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CoinHelperFunctions.hpp"
12#include "CoinFinite.hpp"
13
14#include "CouenneProblem.hpp"
15#include "CouenneObject.hpp"
16#include "CouenneBranchingObject.hpp"
17
18const CouNumber default_alpha  = 0.2;
19const CouNumber default_clamp  = 0.2;
20const CouNumber max_pseudocost = 1000.;
21
22/// Empty constructor
23CouenneObject::CouenneObject ():
24
25  OsiObject (),
26  problem_         (NULL),
27  reference_       (NULL),
28  strategy_        (MID_INTERVAL),
29  jnlst_           (NULL),
30  alpha_           (default_alpha),
31  lp_clamp_        (default_clamp),
32  feas_tolerance_  (feas_tolerance_default),
33  doFBBT_          (true),
34  doConvCuts_      (true),
35  downEstimate_    (0.),
36  upEstimate_      (0.),
37  pseudoMultType_  (INFEASIBILITY) {}
38
39
40/// Constructor with information for branching point selection strategy
41CouenneObject::CouenneObject (CouenneProblem *p, 
42                              exprVar *ref, 
43                              Bonmin::BabSetupBase *base, 
44                              JnlstPtr jnlst):
45  OsiObject (),
46  problem_        (p), 
47  reference_      (ref),
48  strategy_       (MID_INTERVAL),
49  jnlst_          (jnlst),
50  alpha_          (default_alpha),
51  lp_clamp_       (default_clamp),
52  feas_tolerance_ (feas_tolerance_default),
53  doFBBT_         (true),
54  doConvCuts_     (true),
55  downEstimate_   (0.),
56  upEstimate_     (0.),
57  pseudoMultType_ (INFEASIBILITY) {
58
59  // read options
60  setParameters (base);
61
62  // debug output ////////////////////////////////////////////
63
64  if (reference_ &&
65      (reference_ -> Type () == AUX) && 
66      jnlst_ -> ProduceOutput (J_SUMMARY, J_BRANCHING)) {
67
68    printf ("created Expression Object: "); reference_ -> print (); 
69    if (reference_ -> Image ()) {
70      printf (" := "); 
71      reference_ -> Image () -> print ();
72    }
73
74    printf (" with %s strategy [clamp=%g, alpha=%g]\n", 
75            (strategy_ == LP_CLAMPED)   ? "lp-clamped" : 
76            (strategy_ == LP_CENTRAL)   ? "lp-central" : 
77            (strategy_ == BALANCED)     ? "balanced"   : 
78            (strategy_ == MIN_AREA)     ? "min-area"   : 
79            (strategy_ == MID_INTERVAL) ? "mid-point"  : 
80            (strategy_ == NO_BRANCH)    ? "no-branching (null infeasibility)" : 
81                                          "no strategy",
82            lp_clamp_, alpha_);
83  }
84}
85
86
87/// Constructor with lesser information, used for infeasibility only
88CouenneObject::CouenneObject (exprVar *ref, 
89                              Bonmin::BabSetupBase *base, 
90                              JnlstPtr jnlst):
91
92  OsiObject (),
93  problem_        (NULL), 
94  reference_      (ref),
95  strategy_       (MID_INTERVAL),
96  jnlst_          (jnlst),
97  alpha_          (default_alpha),
98  lp_clamp_       (default_clamp),
99  feas_tolerance_ (feas_tolerance_default),
100  doFBBT_         (true),
101  doConvCuts_     (true),
102  downEstimate_   (0.),
103  upEstimate_     (0.),
104  pseudoMultType_ (INFEASIBILITY) {
105
106  // read options
107  setParameters (base);
108}
109
110
111/// Copy constructor
112CouenneObject::CouenneObject (const CouenneObject &src):
113  OsiObject       (src),
114  problem_        (src.problem_),
115  reference_      (src.reference_),
116  strategy_       (src.strategy_),
117  jnlst_          (src.jnlst_),
118  alpha_          (src.alpha_),
119  lp_clamp_       (src.lp_clamp_),
120  feas_tolerance_ (src.feas_tolerance_),
121  doFBBT_         (src.doFBBT_),
122  doConvCuts_     (src.doConvCuts_),
123  downEstimate_   (src.downEstimate_),
124  upEstimate_     (src.upEstimate_),
125  pseudoMultType_ (src.pseudoMultType_) {}
126
127
128/// apply the branching rule
129OsiBranchingObject *CouenneObject::createBranch (OsiSolverInterface *si,
130                                                 const OsiBranchingInformation *info,
131                                                 int way) const {
132
133  // a nonlinear constraint w = f(x) is violated. The infeasibility
134  // should be given by something more elaborate than |w-f(x)|, that
135  // is, it is the minimum, among the two branching nodes, of the
136  // distance from the current optimum (w,x) and the optimum obtained
137  // after convexifying the two subproblems. We call selectBranch for
138  // the purpose, and save the output parameter into the branching
139  // point that should be used later in createBranch.
140
141  if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BRANCHING)) {
142    printf ("CouObj::createBranch on ");
143    reference_ -> print (); printf ("\n");
144  }
145
146  // copy current point into Couenne
147  problem_ -> domain () -> push
148    (problem_ -> nVars (),
149     info -> solution_,
150     info -> lower_,
151     info -> upper_);
152
153  CouNumber 
154    *brPts  = NULL,         // branching point(s)
155    *brDist = NULL;         // distances from current LP point to each
156                            // new convexification (usually two)
157  expression *brVar = NULL; // branching variable
158  int whichWay = 0;
159
160  CouNumber improv;
161
162  if (reference_ -> Image ())  // not used out of debug
163    improv = reference_ -> Image () -> 
164      selectBranch (this, info,                      // input parameters
165                    brVar, brPts, brDist, whichWay); // result: who, where, distance, and direction
166  else {
167    brVar = reference_;
168    brPts  = (double *) realloc (brPts, sizeof (double)); 
169    brDist = (double *) realloc (brDist, 2 * sizeof (double)); 
170
171    double point = info -> solution_ [reference_ -> Index ()];
172
173    *brPts = point;
174    improv = 0.;
175
176    if (point > floor (point)) {improv =                  brDist [0] = point - floor (point);}
177    if (point < ceil  (point)) {improv = CoinMin (improv, brDist [1] = ceil (point) - point);}
178
179    point -= floor (point);
180    whichWay = (point < 0.45) ? TWO_LEFT : (point > 0.55) ? TWO_RIGHT : TWO_RAND;
181  }
182
183  assert (brVar); // MUST have a branching variable
184
185  /*  if        (pseudoMultType_ == INTERVAL) {
186    int index = brVar -> Index ();
187    downEstimate_ = *brPts - info -> lower_ [index];
188    upEstimate_   =          info -> upper_ [index] - *brPts;
189  } else if (pseudoMultType_ == REV_INTERVAL) {
190    int index = brVar -> Index ();
191    downEstimate_ = *brPts - info -> upper_ [index];          // notice upper&lower inverted
192    upEstimate_   =          info -> lower_ [index] - *brPts;
193  } else
194  */
195
196  if (pseudoMultType_ == PROJECTDIST) {
197    downEstimate_ = brDist [0];
198    upEstimate_   = brDist [1];
199  } else setEstimates (info, NULL, brPts);
200
201  /// Debug output
202  if (jnlst_ -> ProduceOutput (J_MOREMATRIX, J_BRANCHING)) {
203    printf ("brpts for "); reference_ -> print (); 
204    if (reference_ -> Image ()) {printf (" := "); reference_ -> Image () -> print ();}
205    printf (" is on "); brVar -> print ();
206    printf (" @ %.12g [%.12g,%.12g]\n", *brPts, 
207            problem_ -> Lb (brVar -> Index ()), 
208            problem_ -> Ub (brVar -> Index ()));
209
210    if (brVar) {
211
212      if (improv <= COUENNE_EPS) {
213        printf ("### warning, infeas = %g for ", improv);
214        reference_ -> print (); 
215        if (reference_ -> Image ()) {printf (":="); reference_ -> Image () -> print ();}
216        printf ("\n");
217      }
218
219      int index = brVar -> Index ();
220      if (info -> lower_ [index] >= 
221          info -> upper_ [index] - COUENNE_EPS) {
222        printf ("### warning, tiny bounding box [%g,%g] for x_%d\n", 
223                info -> lower_ [index],
224                info -> upper_ [index], index);
225      }
226    }
227  }
228
229  // create branching object ///////////////////////////////////////
230  OsiBranchingObject *brObj = new CouenneBranchingObject
231    (si, this, jnlst_, brVar, way, *brPts, doFBBT_, doConvCuts_);
232
233  problem_ -> domain () -> pop (); // Couenne discards current point
234
235  if (brPts)  free (brPts);
236  if (brDist) free (brDist);
237
238  return brObj;
239}
240
241
242/// computes a not-too-bad point where to branch, in the "middle" of an interval
243CouNumber CouenneObject::midInterval (CouNumber x, CouNumber l, CouNumber u) const {
244
245  if (u < l + COUENNE_EPS)
246    return (0.5 * (l + u));
247
248  if      (x<l) x = l;
249  else if (x>u) x = u;
250
251  if   (l < -COUENNE_INFINITY / 10)
252    if (u >  COUENNE_INFINITY / 10) return x; // 0.                                    // ]-inf,+inf[
253    else                            return ((x < -COUENNE_EPS) ? (AGGR_MUL * (-1+x)) : // ]-inf,u]
254                                            (x >  COUENNE_EPS) ? 0. : -AGGR_MUL);
255  else
256    if (u >  COUENNE_INFINITY / 10) return ((x >  COUENNE_EPS) ? (AGGR_MUL *  (1+x)) : // [l,+inf[
257                                            (x < -COUENNE_EPS) ? 0. :  AGGR_MUL);
258    else {                                                                             // [l,u]
259      CouNumber point = alpha_ * x + (1. - alpha_) * (l + u) / 2.;
260      if      ((point-l) / (u-l) < closeToBounds) point = l + (u-l) * closeToBounds;
261      else if ((u-point) / (u-l) < closeToBounds) point = u + (l-u) * closeToBounds;
262      return point;
263    }
264}
265
266
267/// pick branching point based on current strategy
268CouNumber CouenneObject::getBrPoint (funtriplet *ft, CouNumber x0, CouNumber l, CouNumber u) const {
269
270  if ((l < -COUENNE_EPS) && 
271      (u >  COUENNE_EPS) && 
272      (-l/u >= THRES_ZERO_SYMM) &&
273      (-u/l >= THRES_ZERO_SYMM))
274    return 0.;
275
276  CouNumber width = lp_clamp_ * (u-l);
277
278  switch (strategy_) {
279
280  case CouenneObject::MIN_AREA:     return maxHeight   (ft, l, u);
281  case CouenneObject::BALANCED:     return minMaxDelta (ft, l, u);
282  case CouenneObject::LP_CLAMPED:   return CoinMax (l + width, CoinMin (x0, u - width));
283  case CouenneObject::LP_CENTRAL:   return ((x0 < l + width) || (x0 > u - width)) ? (l+u)/2 : x0;
284  case CouenneObject::MID_INTERVAL: return midInterval (x0, l, u);
285  default:
286    printf ("Couenne: unknown branching point selection strategy\n");
287    exit (-1);
288  }
289}
290
291
292/// non-linear infeasibility -- not called by independent's CouenneVarObject
293double CouenneObject::infeasibility (const OsiBranchingInformation *info, int &way) const {
294
295  if (strategy_ == NO_BRANCH) 
296    return (upEstimate_ = downEstimate_ = 0.);
297
298  problem_ -> domain () -> push
299    (problem_ -> nVars (),
300     info -> solution_, 
301     info -> lower_, 
302     info -> upper_);
303
304  double retval = checkInfeasibility (info);
305
306  problem_ -> domain () -> pop ();
307
308  if (pseudoMultType_ == INFEASIBILITY) {
309
310    double point = info -> solution_ [reference_ -> Index ()];
311
312    if (reference_ -> isInteger ()) {
313      if (retval < intInfeasibility (point)) {
314        if (downEstimate_ <       point  - floor (point)) downEstimate_ =       point  - floor (point);
315        if (upEstimate_   < ceil (point) -        point)  upEstimate_   = ceil (point) -        point;
316        retval = intInfeasibility (point);
317      }
318    }
319    else upEstimate_ = downEstimate_ = retval;
320  }
321  else setEstimates (info, &retval, NULL);
322
323  return (reference_ -> isInteger ()) ? 
324    CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
325    retval;
326}
327
328
329/// non-linear infeasibility -- no need for the domain's push
330/// instruction as this is called from
331/// CouenneVarObject::infeasibility()
332double CouenneObject::checkInfeasibility (const OsiBranchingInformation *info) const {
333
334  if (reference_ -> Type () == VAR)
335    return (reference_ -> isInteger ()) ? 
336      intInfeasibility (info -> solution_ [reference_ -> Index ()]) : 0.;
337
338  double 
339    vval = info -> solution_ [reference_ -> Index ()],
340    fval = (*(reference_ -> Image ())) (),
341    denom  = CoinMax (1., reference_ -> Image () -> gradientNorm (info -> solution_));
342
343  // check if fval is a number (happens with e.g. w13 = w12/w5 and w5=0, see test/harker.nl)
344  if (CoinIsnan(fval)) {
345    fval = vval + 1.;
346    denom = 1.;
347  }
348
349  double
350    retval = fabs (vval - fval),
351    ratio = (CoinMax (1., fabs (vval)) / 
352             CoinMax (1., fabs (fval)));
353
354  //printf ("checkinf --> v=%e f=%e den=%e ret=%e ratio=%e\n", vval, fval, denom, retval, ratio);
355
356  if ((ratio < 2)  && 
357      (ratio > .5) &&
358      ((retval /= denom) < CoinMin (COUENNE_EPS, feas_tolerance_)))
359    retval = 0.;
360
361  if (//(retval > 0.) &&
362      (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING))) {
363
364    printf ("  infeas %g: ", retval); 
365    reference_             -> print (); 
366    if (reference_ -> Image ()) {printf (" := "); reference_ -> Image () -> print ();}
367    printf ("\n");
368  }
369
370  return (reference_ -> isInteger ()) ? 
371    CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
372    retval;
373}
374
375
376/// set parameters by reading command line or parameter file
377void CouenneObject::setParameters (Bonmin::BabSetupBase *base) {
378
379  if (!base) return;
380
381  std::string s;
382
383  base -> options () -> GetStringValue ("pseudocost_mult", s, "couenne.");
384
385  if      (s == "interval_lp")     pseudoMultType_ = INTERVAL_LP;
386  else if (s == "interval_lp_rev") pseudoMultType_ = INTERVAL_LP_REV;
387  else if (s == "interval_br")     pseudoMultType_ = INTERVAL_BR;
388  else if (s == "interval_br_rev") pseudoMultType_ = INTERVAL_BR_REV;
389  else if (s == "infeasibility")   pseudoMultType_ = INFEASIBILITY;
390  else if (s == "projectDist")     pseudoMultType_ = PROJECTDIST;
391
392  base -> options() -> GetStringValue ("branch_fbbt",      s, "couenne."); doFBBT_     = (s=="yes");
393  base -> options() -> GetStringValue ("branch_conv_cuts", s, "couenne."); doConvCuts_ = (s=="yes");
394
395  base -> options() -> GetNumericValue ("feas_tolerance", feas_tolerance_, "couenne.");
396
397  std::string brtype;
398  base -> options () -> GetStringValue ("branch_pt_select", brtype, "couenne.");
399
400  if      (brtype == "balanced")    strategy_ = BALANCED;
401  else if (brtype == "lp-clamped")  strategy_ = LP_CLAMPED;
402  else if (brtype == "lp-central")  strategy_ = LP_CENTRAL;
403  else if (brtype == "min-area")    strategy_ = MIN_AREA;
404  else if (brtype == "no-branch")   strategy_ = NO_BRANCH;
405  else if (brtype == "mid-point") {
406    strategy_ = MID_INTERVAL;
407    base -> options () -> GetNumericValue ("branch_midpoint_alpha", alpha_, "couenne.");
408  }
409
410  // accept options for branching rules specific to each operator
411  if (reference_ && reference_ -> Type () == AUX) {
412
413    std::string br_operator = "";
414
415    switch (reference_ -> Image () -> code ()) {
416
417    case COU_EXPRPOW: {
418
419      // begin with default value in case specific exponent are not given
420      base -> options () -> GetStringValue ("branch_pt_select_pow", brtype, "couenne.");
421
422      CouNumber expon = reference_ -> Image () -> ArgList () [1] -> Value ();
423
424      if      (fabs (expon - 2.) < COUENNE_EPS) br_operator = "sqr";
425      else if (fabs (expon - 3.) < COUENNE_EPS) br_operator = "cube";
426      else if (expon             < 0.)          br_operator = "negpow";
427      else                                      br_operator = "pow";
428    } break;
429
430    case COU_EXPRMUL: 
431      br_operator = (reference_ -> Image () -> ArgList () [0] -> Index () !=
432                     reference_ -> Image () -> ArgList () [1] -> Index ()) ?
433        "prod" : "sqr";
434      break;
435    case COU_EXPRINV: br_operator = "negpow"; break;
436    case COU_EXPRDIV: br_operator = "div"; break;
437    case COU_EXPRLOG: br_operator = "log"; break;
438    case COU_EXPREXP: br_operator = "exp"; break;
439    case COU_EXPRSIN:
440    case COU_EXPRCOS: br_operator = "trig"; break;
441    default: break;
442    }
443
444    if (br_operator != "") {
445      // read option
446      char select [40], sel_clamp [40];
447      sprintf (select,    "branch_pt_select_%s", br_operator.c_str ());
448      sprintf (sel_clamp, "branch_lp_clamp_%s",  br_operator.c_str ());
449      base -> options () -> GetStringValue (select, brtype, "couenne.");
450      base -> options () -> GetNumericValue (sel_clamp, lp_clamp_, "couenne.");
451
452      if      (brtype == "balanced")    strategy_ = BALANCED;
453      else if (brtype == "lp-clamped")  strategy_ = LP_CLAMPED;
454      else if (brtype == "lp-central")  strategy_ = LP_CENTRAL;
455      else if (brtype == "min-area")    strategy_ = MIN_AREA;
456      else if (brtype == "no-branch")   strategy_ = NO_BRANCH;
457      else if (brtype == "mid-point") {
458        strategy_ = MID_INTERVAL;
459        base -> options () -> GetNumericValue ("branch_midpoint_alpha", alpha_, "couenne.");
460      }
461    }
462  }
463}
464
465
466/// set up/down estimates based on branching information
467void CouenneObject::setEstimates (const OsiBranchingInformation *info,
468                                  CouNumber *infeasibility,
469                                  CouNumber *brpoint) const {
470
471  int index = reference_ -> Index ();
472
473  CouNumber
474    *up   = &upEstimate_,
475    *down = &downEstimate_,
476     point = 0.;
477
478  ////////////////////////////////////////////////////////////
479  if ((pseudoMultType_ == INTERVAL_LP_REV) ||
480      (pseudoMultType_ == INTERVAL_BR_REV)) {
481
482    up   = &downEstimate_;
483    down = &upEstimate_;
484  }
485
486  ///////////////////////////////////////////////////////////
487  if (info &&
488      ((pseudoMultType_ == INTERVAL_LP) ||
489       (pseudoMultType_ == INTERVAL_LP_REV))) {
490
491    point = info -> solution_ [index];
492
493    CouNumber delta = closeToBounds * (info -> upper_ [index] - info -> lower_ [index]);
494
495    if      (point < info -> lower_ [index] + delta)
496      point        = info -> lower_ [index] + delta;
497    else if (point > info -> upper_ [index] - delta)
498      point        = info -> upper_ [index] - delta;
499  }
500  else if (brpoint &&
501           ((pseudoMultType_ == INTERVAL_BR) ||
502            (pseudoMultType_ == INTERVAL_BR_REV)))
503   
504    point = *brpoint;
505
506  ///////////////////////////////////////////////////////////
507  switch (pseudoMultType_) {
508
509  case INFEASIBILITY: 
510    if (infeasibility) 
511      upEstimate_ = downEstimate_ = *infeasibility;
512    break;
513
514  case INTERVAL_LP:
515  case INTERVAL_LP_REV:
516  case INTERVAL_BR:
517  case INTERVAL_BR_REV:
518    assert (info);
519    *up   = CoinMin (max_pseudocost,         info -> upper_ [index] - point);
520    *down = CoinMin (max_pseudocost, point - info -> lower_ [index]);
521    break;
522
523  case PROJECTDIST: // taken care of in selectBranch procedure
524    break;
525
526  default: 
527    printf ("Couenne: invalid estimate setting procedure\n");
528    exit (-1);
529  }
530
531  if (reference_ -> isInteger ()) {
532    if (downEstimate_ <       point  - floor (point)) downEstimate_ =       point  - floor (point);
533    if (upEstimate_   < ceil (point) -        point)  upEstimate_   = ceil (point) -        point;
534  }
535}
Note: See TracBrowser for help on using the repository browser.