source: stable/0.3/Couenne/src/branch/CouenneObject.cpp @ 524

Last change on this file since 524 was 524, checked in by pbelotti, 9 years ago

added handler for lp_clamp

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