source: stable/0.2/Couenne/src/convex/generateCuts.cpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

File size: 17.3 KB
Line 
1/* $Id: generateCuts.cpp 141 2009-06-03 04:19:19Z pbelotti $
2 *
3 * Name:    generateCuts.cpp
4 * Author:  Pietro Belotti
5 * Purpose: the generateCuts() method of the convexification class CouenneCutGenerator
6 *
7 * (C) Carnegie-Mellon University, 2006-08.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "BonAuxInfos.hpp"
12#include "CglCutGenerator.hpp"
13
14#include "CouenneCutGenerator.hpp"
15#include "CouenneProblem.hpp"
16#include "CouenneSolverInterface.hpp"
17
18#define Couenne_large_bound2 9.99e30
19
20// checks bad cuts against known optimum
21bool isOptimumCut (const CouNumber *opt, OsiCuts &cs, CouenneProblem *p);
22
23// set and lift bound for auxiliary variable associated with objective
24// function
25void fictitiousBound (OsiCuts &cs,
26                      CouenneProblem *p, 
27                      bool action) {     // true before convexifying, false afterwards
28
29  // fictitious bound for initial unbounded lp relaxations
30  const CouNumber large_tol = (Couenne_large_bound2 / 1e6);
31
32  // set trivial dual bound to objective function, if there is none
33
34  int ind_obj = p -> Obj (0) -> Body () -> Index ();
35
36  if (ind_obj < 0) return;
37
38  // we have a single variable objective function
39
40  //int sense = -1; //(p -> Obj (0) -> Sense () == MINIMIZE) ? -1 : 1;
41
42  if (action)
43    //if (sense<0)
44      {if (p -> Lb (ind_obj) < - Couenne_large_bound2) p -> Lb (ind_obj) = - Couenne_large_bound2;}
45  //else         {if (p -> Ub (ind_obj) >   large_bound2) p -> Ub (ind_obj) =   large_bound2;}
46  else
47    //if (sense>0) {if (fabs (p->Ub(ind_obj)-large_bound2)<large_tol) p->Ub(ind_obj)=COUENNE_INFINITY;}
48    //else         
49      {if (fabs (p->Lb(ind_obj)+Couenne_large_bound2)<large_tol) p->Lb(ind_obj) =-COUENNE_INFINITY;}
50}
51
52
53// translate changed bound sparse array into a dense one
54void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged) {
55
56  // convert sparse chg_bds in something handier
57
58  changed  = (int *) realloc (changed, ncols * sizeof (int));
59  nchanged = 0;
60
61  for (register int i=ncols, j=0; i--; j++, chg_bds++)
62    if (chg_bds -> lower() != t_chg_bounds::UNCHANGED ||
63        chg_bds -> upper() != t_chg_bounds::UNCHANGED ) {
64      *changed++ = j;
65      nchanged++;
66    }
67
68  changed -= nchanged;
69  //changed = (int *) realloc (changed, nchanged * sizeof (int));
70}
71
72
73/// get new bounds from parents' bounds + branching rules
74void updateBranchInfo (const OsiSolverInterface &si, CouenneProblem *p, 
75                       t_chg_bounds *chg, const CglTreeInfo &info);
76
77/// a convexifier cut generator
78
79void CouenneCutGenerator::generateCuts (const OsiSolverInterface &si,
80                                        OsiCuts &cs, 
81                                        const CglTreeInfo info) const {
82  const int infeasible = 1;
83
84  int nInitCuts = cs.sizeRowCuts ();
85
86  CouNumber
87    *&realOpt = problem_ -> bestSol (),
88    *saveOptimum = realOpt;
89
90  if (!firstcall_ && realOpt) { 
91
92    // have a debug optimal solution. Check if current bounds
93    // contain it, otherwise pretend it does not exist
94
95    CouNumber *opt = realOpt;
96
97    const CouNumber
98      *sol = si.getColSolution (),
99      *lb  = si.getColLower (),
100      *ub  = si.getColUpper ();
101
102    int objind = problem_ -> Obj (0) -> Body () -> Index ();
103
104    for (int j=0, i=problem_ -> nVars (); i--; j++, opt++, lb++, ub++)
105      if ((j != objind) && 
106          ((*opt < *lb - COUENNE_EPS * (1 + CoinMin (fabs (*opt), fabs (*lb)))) || 
107           (*opt > *ub + COUENNE_EPS * (1 + CoinMin (fabs (*opt), fabs (*ub)))))) {
108       
109        jnlst_ -> Printf (J_VECTOR, J_CONVEXIFYING, 
110                          "out of bounds, ignore x%d = %g [%g,%g] opt = %g\n", 
111                          problem_ -> nVars () - i - 1, *sol, *lb, *ub, *opt);
112
113        // optimal point is not in current bounding box,
114        // pretend realOpt is NULL until we return from this procedure
115        realOpt = NULL;
116        break;
117      }
118  }
119
120  /*static int count = 0;
121  char fname [20];
122  sprintf (fname, "relax_%d", count++);
123  si.writeLp (fname);
124  printf ("writing %s\n", fname);*/
125
126  jnlst_ -> Printf (J_DETAILED, J_CONVEXIFYING,
127                    "generateCuts: level = %d, pass = %d, intree = %d\n",
128                    info.level, info.pass, info.inTree);
129
130  Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
131
132  if (babInfo)
133    babInfo -> setFeasibleNode ();
134
135  double now   = CoinCpuTime ();
136  int    ncols = problem_ -> nVars ();
137
138  // This vector contains variables whose bounds have changed due to
139  // branching, reduced cost fixing, or bound tightening below. To be
140  // used with malloc/realloc/free
141
142  t_chg_bounds *chg_bds = new t_chg_bounds [ncols];
143
144  /*for (int i=0; i < ncols; i++)
145    if (problem_ -> Var (i) -> Multiplicity () <= 0) {
146      chg_bds [i].setLower (t_chg_bounds::UNCHANGED);
147      chg_bds [i].setUpper (t_chg_bounds::UNCHANGED);
148      }*/
149
150  problem_ -> installCutOff (); // install upper bound
151
152  if (firstcall_) {
153
154    // First convexification //////////////////////////////////////
155
156    // OsiSolverInterface is empty yet, no information can be obtained
157    // on variables or bounds -- and none is needed since our
158    // constructor populated *problem_ with variables and bounds. We
159    // only need to update the auxiliary variables and bounds with
160    // their current value.
161
162    for (int i=0; i < ncols; i++) 
163      if (problem_ -> Var (i) -> Multiplicity () > 0) {
164        chg_bds [i].setLower (t_chg_bounds::CHANGED);
165        chg_bds [i].setUpper (t_chg_bounds::CHANGED);
166      }
167
168    // start with FBBT, should take advantage of cutoff found by NLP
169    // run AFTER initial FBBT...
170    if (problem_ -> doFBBT () &&
171        (! (problem_ -> boundTightening (chg_bds, babInfo))))
172      printf ("Couenne: WARNING, first convexification is infeasible\n");
173
174    // For each auxiliary variable replacing the original (nonlinear)
175    // constraints, check if corresponding bounds are violated, and
176    // add cut to cs
177
178    int nnlc = problem_ -> nCons ();
179
180    for (int i=0; i<nnlc; i++) {
181
182      // for each constraint
183      CouenneConstraint *con = problem_ -> Con (i);
184
185      // (which has an aux as its body)
186      int index = con -> Body () -> Index ();
187
188      if ((index >= 0) && 
189          ((con -> Body () -> Type () == AUX) ||
190           (con -> Body () -> Type () == VAR))) {
191
192        // get the auxiliary that is at the lhs
193        exprVar *conaux = problem_ -> Var (index);
194
195        if (conaux &&
196            (conaux -> Type () == AUX) &&
197            (conaux -> Image ()) && 
198            (conaux -> Image () -> Linearity () <= LINEAR)) {
199
200          // reduce density of problem by adding w >= l rather than
201          // ax + b >= l for any linear auxiliary defined as w := ax+b
202
203          double 
204            lb = (*(con -> Lb ())) (), 
205            ub = (*(con -> Ub ())) ();
206
207          OsiColCut newBound;
208          if (lb > -COUENNE_INFINITY) newBound.setLbs (1, &index, &lb);
209          if (ub <  COUENNE_INFINITY) newBound.setUbs (1, &index, &ub);
210
211          cs.insert (newBound);
212
213          // the auxiliary w of constraint w <= b is associated with a
214          // linear expression w = ax: add constraint ax <= b
215          /*conaux -> Image () -> generateCuts (conaux, si, cs, this, chg_bds,
216                                              conaux -> Index (),
217                                              (*(con -> Lb ())) (),
218                                              (*(con -> Ub ())) ());*/
219
220          // take it from the list of the variables to be linearized
221          //
222          // DO NOT decrease multiplicity. Even if it is a linear
223          // term, its bounds can still be used in implied bounds
224          //
225          // Are we sure? That will happen only if its multiplicity is
226          // nonzero, for otherwise this aux is only used here, and is
227          // useless elsewhere
228          //
229          //conaux -> decreaseMult (); // !!!
230        }
231
232        // also, add constraint w <= b
233
234        // not now, do it later
235
236//      // if there exists violation, add constraint
237//      CouNumber l = con -> Lb () -> Value (),
238//                u = con -> Ub () -> Value ();
239
240//      // tighten bounds in Couenne's problem representation
241//      problem_ -> Lb (index) = CoinMax (l, problem_ -> Lb (index));
242//      problem_ -> Ub (index) = CoinMin (u, problem_ -> Ub (index));
243
244      } else { // body is more than just a variable, but it should be
245               // linear. If so, generate equivalent linear cut
246
247        assert (false); // TODO
248      }
249    }
250
251    if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
252      if (cs.sizeRowCuts ()) {
253        jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint row cuts\n",
254                          cs.sizeRowCuts ());
255        for (int i=0; i<cs.sizeRowCuts (); i++) 
256          cs.rowCutPtr (i) -> print ();
257      }
258      if (cs.sizeColCuts ()) {
259        jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne: %d constraint col cuts\n",
260                          cs.sizeColCuts ());
261        for (int i=0; i<cs.sizeColCuts (); i++) 
262          cs.colCutPtr (i) -> print ();
263      }
264    }
265  } else {
266
267    // use new optimum as lower bound for variable associated w/objective
268    int indobj = problem_ -> Obj (0) -> Body () -> Index ();
269
270    assert (indobj >= 0);
271
272    // transmit solution from OsiSolverInterface to problem
273    problem_ -> domain () -> push
274      (problem_ -> nVars (),
275       si. getColSolution (), 
276       si. getColLower    (),
277       si. getColUpper    ());
278
279    if (indobj >= 0) {
280
281      // Use current value of objvalue's x as a lower bound for bound
282      // tightening
283      double lp_bound = problem_ -> domain () -> x (indobj);
284
285      //if (problem_ -> Obj (0) -> Sense () == MINIMIZE)
286      {if (lp_bound > problem_ -> Lb (indobj)) problem_ -> Lb (indobj) = lp_bound;}
287           //else {if (lp_bound < problem_ -> Ub (indobj)) problem_ -> Ub (indobj) = lp_bound;}
288    }
289
290    updateBranchInfo (si, problem_, chg_bds, info); // info.depth >= 0 || info.pass >= 0
291  }
292
293  // restore constraint bounds before tightening and cut generation
294  for (int i = problem_ -> nCons (); i--;) {
295
296    // for each constraint
297    CouenneConstraint *con = problem_ -> Con (i);
298
299    // (which has an aux as its body)
300    int index = con -> Body () -> Index ();
301
302    if ((index >= 0) && 
303        ((con -> Body () -> Type () == AUX) ||
304         (con -> Body () -> Type () == VAR))) {
305
306      // if there exists violation, add constraint
307      CouNumber l = con -> Lb () -> Value (),   
308        u = con -> Ub () -> Value ();
309
310      // tighten bounds in Couenne's problem representation
311      problem_ -> Lb (index) = CoinMax (l, problem_ -> Lb (index));
312      problem_ -> Ub (index) = CoinMin (u, problem_ -> Ub (index));
313    }
314  }
315
316  problem_ -> installCutOff (); // install upper bound
317
318  fictitiousBound (cs, problem_, false); // install finite lower bound, if currently -inf
319
320  int *changed = NULL, nchanged;
321
322  // Bound tightening ///////////////////////////////////////////
323
324  // do bound tightening only at first pass of cutting plane in a node
325  // of BB tree (info.pass == 0) or if first call (creation of RLT,
326  // info.pass == -1)
327
328  try {
329
330    // Bound tightening ////////////////////////////////////
331
332    /*printf ("== BT ================\n");
333    for (int i = 0; i < problem_ -> nVars (); i++)
334      if (problem_ -> Var (i) -> Multiplicity () > 0)
335        printf ("%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
336                problem_ -> X  (i), problem_ -> Lb (i), problem_ -> Ub (i));
337                printf("=============================\n");*/
338
339    // Reduced Cost BT -- to be done first to use rcost correctly
340    if (!firstcall_  &&                         // have a linearization already
341        problem_ -> doRCBT () &&                // authorized to do reduced cost tightening
342        problem_ -> redCostBT (&si, chg_bds) && // some variables were tightened with reduced cost
343        !(problem_ -> btCore (chg_bds)))        // in this case, do another round of FBBT
344      throw infeasible;
345
346    // FBBT
347    if (problem_ -> doFBBT () && 
348        //(info.pass <= 0) && // do it in subsequent rounds too
349        (! (problem_ -> boundTightening (chg_bds, babInfo))))
350      throw infeasible;
351
352    // OBBT
353    if (!firstcall_ && // no obbt if first call (there is no LP to work with)
354        problem_ -> obbt (this, si, cs, info, babInfo, chg_bds) < 0)
355      throw infeasible;
356
357    // Bound tightening done /////////////////////////////
358
359    if ((problem_ -> doFBBT () ||
360         problem_ -> doOBBT () ||
361         problem_ -> doABT  ()) &&
362        (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING))) {
363
364      jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== after bt =============\n");
365      for (int i = 0; i < problem_ -> nVars (); i++)
366        if (problem_ -> Var (i) -> Multiplicity () > 0)
367          jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
368                         problem_ -> X  (i), problem_ -> Lb (i), problem_ -> Ub (i));
369      jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
370    }
371
372    // Generate convexification cuts //////////////////////////////
373
374    sparse2dense (ncols, chg_bds, changed, nchanged);
375
376    double *nlpSol;
377
378    //--------------------------------------------
379
380    if (babInfo && ((nlpSol = const_cast <double *> (babInfo -> nlpSolution ())))) {
381
382      // Aggressive Bound Tightening ////////////////////////////////
383
384      int logAbtLev = problem_ -> logAbtLev ();
385
386      if (problem_ -> doABT () &&           // flag is checked, AND
387          ((logAbtLev != 0) ||                // (parameter is nonzero OR
388           (info.level == 0)) &&              //  we are at root node), AND
389          (info.pass == 0) &&               // at first round of cuts, AND
390          ((logAbtLev < 0) ||                 // (logAbtLev = -1, OR
391           (info.level <= logAbtLev) ||       //  depth is lower than COU_OBBT_CUTOFF_LEVEL, OR
392           (CoinDrand48 () <                  //  probability inversely proportional to the level)
393            pow (2., (double) logAbtLev - (info.level + 1))))) {
394
395        jnlst_ -> Printf(J_ITERSUMMARY, J_CONVEXIFYING,"  performing ABT\n");
396        if (! (problem_ -> aggressiveBT (nlp_, chg_bds, babInfo)))
397          throw infeasible;
398
399        sparse2dense (ncols, chg_bds, changed, nchanged);
400      }
401
402      // obtain solution just found by nlp solver
403
404      // Auxiliaries should be correct. solution should be the one found
405      // at the node even if not as good as the best known.
406
407      // save violation flag and disregard it while adding cut at NLP
408      // point (which are not violated by the current, NLP, solution)
409      bool save_av = addviolated_;
410      addviolated_ = false;
411
412      // save values
413      problem_ -> domain () -> push
414        (problem_ -> nVars (), 
415         problem_ -> domain () -> x  (), 
416         problem_ -> domain () -> lb (), 
417         problem_ -> domain () -> ub (), false);
418
419      // fill originals with nlp values
420      CoinCopyN (nlpSol, problem_ -> nOrigVars (), problem_ -> domain () -> x ());
421      //problem_ -> initAuxs ();
422
423      problem_ -> getAuxs (problem_ -> domain () -> x ());
424
425      if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
426        jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on NLP =============\n");
427        for (int i = 0; i < problem_ -> nVars (); i++)
428          if (problem_ -> Var (i) -> Multiplicity () > 0)
429            jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
430                           problem_ -> X  (i),
431                           problem_ -> Lb (i),
432                           problem_ -> Ub (i));
433        jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
434      }
435
436      problem_ -> domain () -> current () -> isNlp () = true;
437      genRowCuts (si, cs, nchanged, changed, chg_bds);  // add cuts
438
439      problem_ -> domain () -> pop (); // restore point
440
441      addviolated_ = save_av;     // restore previous value
442
443      //    if (!firstcall_) // keep solution if called from extractLinearRelaxation()
444      babInfo -> setHasNlpSolution (false); // reset it after use //AW HERE
445    } else {
446
447      if (jnlst_ -> ProduceOutput (J_VECTOR, J_CONVEXIFYING)) {
448        jnlst_ -> Printf(J_VECTOR, J_CONVEXIFYING,"== genrowcuts on LP =============\n");
449        for (int i = 0; i < problem_ -> nVars (); i++)
450          if (problem_ -> Var (i) -> Multiplicity () > 0)
451            jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"%4d %+20.8g [%+20.8g,%+20.8g]\n", i,
452                           problem_ -> X  (i),
453                           problem_ -> Lb (i),
454                           problem_ -> Ub (i));
455        jnlst_->Printf(J_VECTOR, J_CONVEXIFYING,"=============================\n");
456      }
457
458      genRowCuts (si, cs, nchanged, changed, chg_bds);
459    }
460
461    // change tightened bounds through OsiCuts
462    if (nchanged)
463      genColCuts (si, cs, nchanged, changed);
464
465    if (firstcall_ && (cs.sizeRowCuts () >= 1))
466      jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
467                     "Couenne: %d initial row cuts\n", cs.sizeRowCuts ());
468
469    if (realOpt && // this is a good time to check if we have messed up with the optimal solution
470        isOptimumCut (realOpt, cs, problem_))
471      jnlst_->Printf(J_ITERSUMMARY, J_CONVEXIFYING,
472                     "\n\n CUT OPTIMUM\n\n");
473  }
474
475  catch (int exception) {
476
477    if ((exception == infeasible) && (!firstcall_)) {
478
479      jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,
480                        "Couenne: Infeasible node\n");
481
482      OsiColCut *infeascut = new OsiColCut;
483
484      if (infeascut) {
485        int i=0;
486        double upper = -1., lower = +1.;
487        infeascut -> setLbs (1, &i, &lower);
488        infeascut -> setUbs (1, &i, &upper);
489        cs.insert (infeascut);
490        delete infeascut;
491      }
492    }
493
494    if (babInfo) // set infeasibility to true in order to skip NLP heuristic
495      babInfo -> setInfeasibleNode ();
496  }
497
498  delete [] chg_bds;
499
500  if (changed) 
501    free (changed);
502
503  if (firstcall_) {
504
505    jnlst_ -> Printf (J_SUMMARY, J_CONVEXIFYING, 
506                      "Couenne: %d cuts (%d row, %d col) for linearization\n", 
507                      cs.sizeRowCuts () + cs.sizeColCuts (),
508                      cs.sizeRowCuts (),  cs.sizeColCuts ());
509
510    fictitiousBound (cs, problem_, true);
511    firstcall_  = false;
512    ntotalcuts_ = nrootcuts_ = cs.sizeRowCuts ();
513  } else { 
514    problem_ -> domain () -> pop ();
515
516    ntotalcuts_ += (cs.sizeRowCuts () - nInitCuts);
517
518    if (saveOptimum)
519      realOpt = saveOptimum; // restore debug optimum
520  }
521
522  septime_ += CoinCpuTime () - now;
523
524  if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_CONVEXIFYING)) {
525
526    if (cs.sizeColCuts ()) {
527      jnlst_ -> Printf (J_ITERSUMMARY, J_CONVEXIFYING,"Couenne col cuts:\n");
528      for (int i=0; i<cs.sizeColCuts (); i++) 
529        cs.colCutPtr (i) -> print ();
530    }
531  }
532
533  if (!(info.inTree)) 
534    rootTime_ = CoinCpuTime ();
535}
Note: See TracBrowser for help on using the repository browser.