source: trunk/Couenne/src/convex/generateCuts.cpp @ 112

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

isolate restoreUnused to when there is at least one unused -- takes care of reformulated stockcycle problem (thanks to C. D'Ambrosio for the bug report). Check linear terms before adding a full exprGroup -- need a constructor pilot. Quicker solution feasibility check (with initial integrality check).

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