source: trunk/Couenne/src/bound_tightening/aggressiveBT.cpp @ 589

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

restored scanf of cutoff. FBBT on real bounds within non-CouenneObject? branching (thanks to F. Margot for noticing).

  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
Line 
1/* $Id: aggressiveBT.cpp 589 2011-05-31 08:55:27Z pbelotti $
2 *
3 * Name:    aggressiveBT.cpp
4 * Author:  Pietro Belotti
5 * Purpose: aggressive bound tightening -- fake bounds in variables to
6 *          exclude parts of the solution space through fathoming on
7 *          bounds/infeasibility
8 *
9 * (C) Carnegie-Mellon University, 2007-10.
10 * This file is licensed under the Eclipse Public License (EPL)
11 */
12
13#include "CouenneCutGenerator.hpp"
14
15#include "BonTNLPSolver.hpp"
16#include "BonNlpHeuristic.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "BonCouenneInfo.hpp"
19
20#include "CouenneProblem.hpp"
21#include "CouenneExprVar.hpp"
22#include "CouenneProblemElem.hpp"
23
24using namespace Couenne;
25
26namespace Bonmin {
27  class OsiTMINLPInterface;
28  class BabInfo;
29  class TNLPSolver;
30}
31
32#define MAX_ABT_ITER           1  // max # aggressive BT iterations
33#define THRES_ABT_IMPROVED     0  // only continue ABT if at least these bounds have improved
34#define THRES_ABT_ORIG       100  // only do ABT on auxiliaries if they are less originals than this
35
36static double distanceToBound (int n, const double* xOrig,
37                               const double* lower, const double* upper) {
38
39  double Xdist = 0.;
40
41  for (int i=0, j=n; j--; i++) {
42
43    CouNumber
44      diff, 
45      xO = xOrig [i];
46
47    if      ((diff = lower [i] - xO) > 0) Xdist += diff;
48    else if ((diff = xO - upper [i]) > 0) Xdist += diff;
49  }
50
51  return Xdist;
52}
53
54
55// Probing: for each variable, fake new bounds [l,b] or [b,u], with
56// given b, and apply bound tightening. If the interval is fathomed on
57// bounds or on infeasibility, the complementary bound interval is a
58// valid tightening.
59
60bool CouenneProblem::aggressiveBT (Bonmin::OsiTMINLPInterface *nlp,
61                                   t_chg_bounds *chg_bds, 
62                                   const CglTreeInfo &info,
63                                   Bonmin::BabInfo * babInfo) const {
64
65  Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "Aggressive FBBT\n");
66
67  if (info.level <= 0 && !(info.inTree))  {
68    jnlst_ -> Printf (J_ERROR, J_COUENNE, "Probing: ");
69    fflush (stdout);
70  }
71
72  CouenneInfo* couInfo =
73    dynamic_cast <CouenneInfo *> (babInfo);
74
75  int  ncols  = nVars ();
76  bool retval = false;
77
78  CouNumber
79    *olb = new CouNumber [ncols],
80    *oub = new CouNumber [ncols];
81
82  // save current bounds
83  CoinCopyN (Lb (), ncols, olb);
84  CoinCopyN (Ub (), ncols, oub);
85
86  // Find the solution that is closest to the current bounds
87  // TODO: Also check obj value
88  SmartPtr<const CouenneInfo::NlpSolution> closestSol;
89  double dist = 1e50;
90
91  const std::list<SmartPtr<const CouenneInfo::NlpSolution> >& solList =
92    couInfo->NlpSolutions();
93
94  for (std::list<SmartPtr<const CouenneInfo::NlpSolution> >::const_iterator
95         i = solList.begin();
96       i != solList.end(); i++) {
97    assert(nOrigVars_ == (*i)->nVars());
98    const double thisDist = distanceToBound(nOrigVars_, (*i)->solution(), olb, oub);
99    if (thisDist < dist) {
100      closestSol = *i;
101      dist = thisDist;
102    }
103  }
104
105  jnlst_ -> Printf (J_VECTOR, J_BOUNDTIGHTENING, "best dist = %e\n", dist);
106
107  // If this solution is not sufficiently inside the bounds, we solve the NLP now
108  if (dist > 0.1) { // TODO: Find tolerance
109
110    // find integer initial point /////////////////////////
111
112    int nvars = nVars ();
113
114    double *lower = new double [nvars];
115    double *upper = new double [nvars];
116
117    CoinFillN (lower, nvars, -COUENNE_INFINITY);
118    CoinFillN (upper, nvars,  COUENNE_INFINITY);
119
120    CoinCopyN (nlp -> getColLower (), nOrigVars_, lower);
121    CoinCopyN (nlp -> getColUpper (), nOrigVars_, upper);
122
123    double *Y = new double [nvars];
124
125    CoinFillN (Y,    nvars,      0.);
126    CoinCopyN (X (), nOrigVars_, Y);
127
128    if (getIntegerCandidate (nlp -> getColSolution (), Y, lower, upper) < 0) {
129
130      jnlst_ -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: find NLP point in ABT failed\n");
131      retval = true;
132
133    } else {
134
135      //////////////////////////////////////////////////////
136
137      nlp -> setColLower    (lower);
138      nlp -> setColUpper    (upper);
139      nlp -> setColSolution (Y);
140
141      try {
142        nlp -> options () -> SetNumericValue ("max_cpu_time", CoinMax (0., getMaxCpuTime () - CoinCpuTime ()));
143        nlp -> initialSolve ();
144      }
145      catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
146   
147      delete [] Y;
148      delete [] lower;
149      delete [] upper;
150
151      if (nlp->isProvenOptimal()) {
152        closestSol = new CouenneInfo::NlpSolution
153          (nOrigVars_, nlp->getColSolution(), nlp->getObjValue());
154        couInfo->addSolution(closestSol);
155        dist = 0.;
156      }
157      else {
158        jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: NLP solve in ABT failed\n");
159        retval = true;
160      }
161    }
162  }
163
164  int nTotImproved = 0;
165
166  if (!retval && (dist < 1e10)) {
167
168    retval = true;
169
170    // X is now the NLP solution, but in a low-dimensional space. We
171    // have to get the corresponding point in higher dimensional space
172    // through getAuxs()
173
174    double *X = new double [ncols];
175    CoinCopyN (closestSol->solution(), nOrigVars_, X);
176    getAuxs (X);
177
178    // create a new, fictitious, bound bookkeeping structure
179    t_chg_bounds *f_chg = new t_chg_bounds [ncols];
180
181    if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
182      //    CouNumber cutoff = getCutOff ();
183      int       objind = Obj (0) -> Body  () -> Index ();
184      for (int i=0; i<nOrigVars_; i++)
185        Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
186                        "   %2d %+20g [%+20g %+20g]\n",
187                        i, X [i], Lb (i), Ub (i));
188      Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
189                      "-------------\nAggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
190                      Lb (objind), getCutOff (), ncols);
191    }
192
193    int improved, second, iter = 0;
194
195    // Repeatedly fake tightening bounds on both sides of every variable
196    // to concentrate around current NLP point.
197    //
198    // MAX_ABT_ITER is the maximum # of outer cycles. Each call to
199    // fake_tighten in turn has an iterative algorithm for a
200    // derivative-free, uni-dimensional optimization problem on a
201    // monotone function.
202
203    do {
204
205      bool maxTimeReached = false;
206
207      improved = 0;
208
209      // scan all variables
210      for (int i=0; i<ncols; i++) {
211
212        if (CoinCpuTime () > maxCpuTime_) {
213          maxTimeReached = true; // avoids getrusage...
214          break;
215        }
216
217        int index = evalOrder (i);
218
219        if (Var (index) -> Multiplicity () <= 0) 
220          continue;
221
222        // AW: We only want to do the loop that temporarily changes
223        // bounds around the NLP solution only for those points from the
224        // NLP solution (no auxiliary vars)?
225
226        // PBe: if we do want that, index should be initialized as i, as
227        // evalOrder gives a variable index out of an array index.
228
229        // PBe: That makes a lot of sense when problems are really
230        // big. Instances arki000[24].nl spend a lot of time here
231
232        if ((nOrigVars_ < THRES_ABT_ORIG) || (index < nOrigVars_)) {
233
234          // if (index == objind) continue; // don't do it on objective function
235
236          improved = 0;
237
238          if ((variables_ [index] -> sign () != expression::AUX_GEQ) &&
239              (X [index] >= Lb (index) + COUENNE_EPS)) {
240
241            Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
242                            "------------- tighten left x%d\n", index);
243
244            // tighten on left
245            if ((improved = fake_tighten (0, index, X, olb, oub, chg_bds, f_chg)) < 0) {
246              retval = false;
247              break;
248            }
249          }
250
251          second = 0;
252
253          if (retval && (variables_ [index] -> sign () != expression::AUX_LEQ) &&
254              (X [index] <= Ub (index) - COUENNE_EPS)) {
255            Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
256                            "------------- tighten right x%d\n", index);
257
258            // tighten on right
259            if ((second = fake_tighten (1, index, X, olb, oub, chg_bds, f_chg) < 0)) {
260              retval = false;
261              break;
262            }
263          }
264
265          improved += second;
266          nTotImproved += improved;
267        }
268      }
269
270      if (maxTimeReached)
271        break;
272
273    } while (retval && (improved > THRES_ABT_IMPROVED) && (iter++ < MAX_ABT_ITER));
274
275    // store new valid bounds, or restore old ones if none changed
276    CoinCopyN (olb, ncols, Lb ());
277    CoinCopyN (oub, ncols, Ub ());
278
279    if (Jnlst()->ProduceOutput(J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
280
281      Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,"------------------\n");
282
283      if (!retval) Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
284                                   "Couenne infeasible node from aggressive BT\n");
285
286      int objind = Obj (0) -> Body  () -> Index ();
287
288      Jnlst()->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING,
289                      "-------------\ndone Aggressive BT. Current bound = %g, cutoff = %g, %d vars\n", 
290                      Lb (objind), getCutOff (), ncols);
291
292      if (Jnlst()->ProduceOutput(J_DETAILED, J_BOUNDTIGHTENING))
293        for (int i=0; i<nOrigVars_; i++)
294          printf("   x%02d [%+20g %+20g]  | %+20g\n",
295                 i, Lb (i), Ub (i), X [i]);
296
297      if (Jnlst()->ProduceOutput(J_MOREDETAILED, J_BOUNDTIGHTENING))
298        for (int i=nOrigVars_; i<ncols; i++)
299          printf ("   w%02d [%+20g %+20g]  | %+20g\n", i, Lb (i), Ub (i), X [i]);
300    }
301
302    delete [] X;
303    delete [] f_chg;
304   
305  } else
306
307    if ((dist > 1e10) && !retval)
308      jnlst_ -> Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "TODO: Don't have point for ABT\n");
309
310  delete [] olb;
311  delete [] oub;
312
313  if (info.level <= 0 && !(info.inTree))  {
314    if (!retval) jnlst_ -> Printf (J_ERROR, J_COUENNE, "done (infeasible)\n");
315    else         jnlst_ -> Printf (J_ERROR, J_COUENNE, "done (%d improved bounds)\n", nTotImproved);
316  }
317
318  return retval;// && btCore (psi, cs, chg_bds, babInfo, true); // !!!
319  //return retval && btCore (psi, cs, chg_bds, babInfo, true);
320}
Note: See TracBrowser for help on using the repository browser.