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

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

intInfeasibility() checks lower/upper bounds to avoid non-null infeasibility upon out-of-bounds values, which triggers infeasibility in Couenne but not in Cbc, and this in turn makes the BB loop

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