source: trunk/Couenne/src/disjunctive/generateDisjCuts.cpp @ 945

Last change on this file since 945 was 945, checked in by stefan, 8 years ago

generateCuts in Cgl >= 0.58 will not be const anymore

  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 KB
Line 
1/* $Id: generateDisjCuts.cpp 945 2013-04-06 20:25:21Z stefan $
2 *
3 * Name:    generateDisjCuts.cpp
4 * Author:  Pietro Belotti
5 * Purpose: separation method for disjunctive cuts
6 *
7 * (C) Carnegie-Mellon University, 2008-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include <vector>
12#include "CoinTime.hpp"
13
14#include "CouenneCutGenerator.hpp"
15#include "CouenneDisjCuts.hpp"
16#include "CouenneProblem.hpp"
17#include "CouenneInfeasCut.hpp"
18
19using namespace Ipopt;
20using namespace Couenne;
21
22/// generate disjunctive cuts
23void CouenneDisjCuts::generateCuts (const OsiSolverInterface &si, 
24                                    OsiCuts &cs, 
25                                    const CglTreeInfo info)
26#if CGL_VERSION_MAJOR == 0 && CGL_VERSION_MINOR <= 57
27  const
28#endif
29{
30
31  // Check if cs contains only one cut and if it is of the form 1 <=
32  // x0 <= -1. That means a previous cut generator has determined that
33  // this node is infeasible and we shouldn't take the pain of running
34  // this CGL.
35
36  if (isWiped (cs))
37    return;
38 
39  if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) 
40    printf ("--- generateDisjCuts: level = %d, pass = %d, intree = %d [%d]\n",
41            info.level, info.pass, info.inTree, depthStopSeparate_);
42
43  if ((depthStopSeparate_ >= 0) &&        // if -1 no limit on depth
44      (info.level > depthStopSeparate_))  // check if too deep for adding these cuts
45    return;
46
47  int nInitCuts = nrootcuts_;
48
49  if ((info.level <= 0) && !(info.inTree)) {
50
51    jnlst_ -> Printf (J_ERROR, J_COUENNE, 
52                      "Disjunctive cuts (root node): ");
53    fflush (stdout);
54  }
55
56  double time = CoinCpuTime ();
57
58  // use clone of solver interface
59  OsiSolverInterface *csi = si.clone ();
60
61  int
62    initRowCuts = cs.sizeRowCuts (),
63    initColCuts = cs.sizeColCuts ();
64
65  // get disjunctions /////////////////////////////
66  //
67  // start:
68  //
69  // consider problem Ax <= a_0, x in [l,u]
70  //
71  // // A and a_0 may be, depending on depth of BB tree and size of the
72  // // problem, from the linearization at the root node (in this case
73  // // a clone() is needed at the first iteration and branching rules
74  // // should be applied explicitly) or the current LP as passed by
75  // // si.
76  //
77  // Get set of disjunctions (in the form of OsiBranchingObjects) in
78  // number limited by depth, problem size, and sorted according to
79  // branching method (HotInfo if strong branching?)
80  //
81  // preprocessing /////////////////////////////////
82  //
83  // for each disjunction (x_i <= or >= x_i^d) {
84  //
85  //   1) apply left  disj., bound tightening returns x in [l_1,u_1]
86  //   2) apply right disj., bound tightening returns x in [l_2,u_2]
87  //
88  //   3) if both infeasible, bail out
89  //   4) if one infeasible only, save column cut, apply tightening to
90  //      si and goto start
91  // }
92  //
93  // CGLP //////////////////////////////////////////
94  //
95  // for each disjunction (x_i <= or >= x_i^d) {
96  //
97  //   5) get cuts Bx <= b_0 for left  disj.
98  //   6) get cuts Cx <= c_0 for right disj.
99  //
100  //   7) if both infeasible, bail out
101  //   8) if one infeasible only, save column cut, apply tightening to
102  //      si and continue
103  //
104  //   9) generate CGLP:
105  //         consider subset of Ax <= a_0:
106  //           a) active only, or
107  //           b) {root linearization cuts} union {currently active}, or
108  //           c) all cuts (!)
109  //
110  //   10) solve CGLP
111  //
112  //   11) add corresponding cut, possibly to CGLP as well?
113  // }
114
115  int maxDisj = (initDisjNumber_ >= 0) ? 
116    CoinMin ((int) (csi -> numberObjects () * initDisjPercentage_), initDisjNumber_) :
117    (int) (csi -> numberObjects () * initDisjPercentage_);
118
119  // number of disjunctions to consider (branching objects)
120  numDisjunctions_ = (depthLevelling_ < 0 || info.level < depthLevelling_) ? 
121    (int) (maxDisj) :
122    (int) (maxDisj / (2 + info.level - depthLevelling_));
123
124  if (numDisjunctions_ < 1) numDisjunctions_ = 1;
125
126  const int 
127    exc_infeasible = 1, 
128    exc_normal     = 2, 
129    max_iterations = 1;
130
131  couenneCG_ -> Problem () -> domain () -> push (&si, &cs, false);
132
133  std::vector <std::pair <OsiCuts *, OsiCuts *> > disjunctions;
134
135  bool infeasNode = false;
136
137  try {
138
139    // get disjunctions (rows and cols) ////////////////////////////////////////////////////
140
141    bool start_over;
142    int iterations = 0;
143
144    do { // repeat as long as preprocessing or separation of convCuts
145         // shrink the whole problem
146
147      ++iterations;
148
149      start_over = false;
150
151      // preprocess, get column cuts (disjoint bounding boxes)
152      int result = getDisjunctions (disjunctions, *csi, cs, info);
153
154      if      (result == COUENNE_INFEASIBLE) throw exc_infeasible; // fathom node
155      else if (result == COUENNE_TIGHTENED && 
156               iterations < max_iterations) start_over = true;     // tightened node
157
158      if (disjunctions.empty ())
159        throw exc_normal;
160
161      // generate convexification cuts for each disjunction
162      for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
163           disjI != disjunctions.end (); ++disjI) {
164
165        if (CoinCpuTime () > couenneCG_ -> Problem () -> getMaxCpuTime ()) {
166          start_over = false;
167          break;
168        }
169
170        // separate on single disjunction
171
172        // left
173        result = separateWithDisjunction (disjI -> first, *csi, cs, info);
174        if      (result == COUENNE_INFEASIBLE) throw exc_infeasible;           // fathom node
175        else if (result == COUENNE_TIGHTENED && iterations < max_iterations) { // tightened node
176          start_over = true;
177          break;
178        }
179
180        // right
181        result = separateWithDisjunction (disjI -> second, *csi, cs, info);
182        if      (result == COUENNE_INFEASIBLE) throw exc_infeasible;           // fathom node
183        else if (result == COUENNE_TIGHTENED && iterations < max_iterations) { // tightened node
184          start_over = true;
185          break;
186        }
187      }
188
189      if (start_over) {
190
191        for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
192             disjI != disjunctions.end (); ++disjI) {
193          delete disjI -> first;
194          delete disjI -> second;
195        }
196
197        disjunctions.erase (disjunctions.begin (), disjunctions.end ());
198      }
199
200      if (!start_over && jnlst_ -> ProduceOutput (J_VECTOR, J_DISJCUTS))
201
202        // generate convexification cuts for each disjunction
203        for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
204             disjI != disjunctions.end (); ++disjI) {
205
206          printf ("=========================== CUTS for the LEFT part\n");
207          for (int i=0; i<disjI->first->sizeColCuts (); i++) disjI->first->colCutPtr(i)->print();
208          for (int i=0; i<disjI->first->sizeRowCuts (); i++) disjI->first->rowCutPtr(i)->print();
209          printf ("=========================== CUTS for the RIGHT part\n");
210          for (int i=0; i<disjI->second->sizeColCuts (); i++) disjI->second->colCutPtr(i)->print();
211          for (int i=0; i<disjI->second->sizeRowCuts (); i++) disjI->second->rowCutPtr(i)->print();
212          printf ("===========================\n");
213        }
214
215    } while (start_over && (iterations < max_iterations));
216
217    // si contains the tightest bounding box. Use it to update
218    // CouenneProblem's bounds AND add to cs
219
220    // already done above
221
222    // CGLP //////////////////////////////////////////////////////////////////////////
223
224    // maybe one last FBBT before big CGLP?
225
226    // generate all cuts
227    if (generateDisjCuts (disjunctions, *csi, cs, info) == COUENNE_INFEASIBLE) // node can be fathomed
228      throw exc_infeasible;
229  }
230
231  catch (int exception) {
232
233    if (exception == exc_infeasible) { // add infeasible column cut 1 <= x_0 <= -1
234
235      jnlst_ -> Printf (J_DETAILED, J_DISJCUTS, "--- Disjunctive Cut separator: infeasible node\n");
236      WipeMakeInfeas (cs);
237      infeasNode = true;
238    }
239  }
240
241  // cleanup
242  for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
243       disjI != disjunctions.end (); ++disjI) {
244
245    delete disjI -> first;
246    delete disjI -> second;
247  }
248
249  couenneCG_ -> Problem () -> domain () -> pop ();
250
251  if (!infeasNode) {
252
253    // tighten bounds of si based on those tightened of csi
254
255    CoinPackedVector
256      tighterLower, 
257      tighterUpper;
258
259    const double 
260      *oldLo = si. getColLower (),   *newLo = csi -> getColLower (), 
261      *oldUp = si. getColUpper (),   *newUp = csi -> getColUpper ();
262
263    int ncols = si.getNumCols ();
264
265    for (int i=0; i<ncols; i++, newLo++, newUp++) {
266
267      if (*newLo > *oldLo++ + COUENNE_EPS) tighterLower.insert (i, *newLo);
268      if (*newUp < *oldUp++ - COUENNE_EPS) tighterUpper.insert (i, *newUp);
269    }
270
271    if ((tighterLower.getNumElements () > 0) || 
272        (tighterUpper.getNumElements () > 0)) {
273      OsiColCut tighter;
274      tighter.setLbs (tighterLower);
275      tighter.setUbs (tighterUpper);
276      if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
277        printf ("tightened bounds in disjunctive cuts:");
278        tighter.print ();
279      }
280      cs.insert (tighter);
281    }
282
283    int deltaNcuts = 
284      cs.sizeRowCuts () - initRowCuts + 
285      cs.sizeColCuts () - initColCuts;
286
287    if (info.level <= 0 && !(info.inTree)) 
288      nrootcuts_ += deltaNcuts;
289    ntotalcuts_ += deltaNcuts;
290
291    if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
292
293      if (cs.sizeRowCuts()>initRowCuts) printf ("added %d row cuts\n", cs.sizeRowCuts () - initRowCuts);
294      if (cs.sizeColCuts()>initColCuts) printf ("added %d col cuts\n", cs.sizeColCuts () - initColCuts);
295    }
296
297    if ((info.level <= 0) && !(info.inTree))
298      jnlst_ -> Printf (J_ERROR, J_COUENNE, 
299                        "%d cuts\n", CoinMax (0, nrootcuts_ - nInitCuts));
300    else if (deltaNcuts)
301      jnlst_ -> Printf (J_WARNING, J_COUENNE, 
302                        "In-BB disjunctive cuts: %d row cuts, %d col cuts\n",
303                        cs.sizeRowCuts () - initRowCuts,
304                        cs.sizeColCuts () - initColCuts);
305  }
306
307  if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
308
309    printf ("Disjunctive cuts (%d row + %d col):\n", cs.sizeRowCuts () - initRowCuts, cs.sizeColCuts () - initColCuts);
310    for (int i=initRowCuts; i<cs.sizeRowCuts (); ++i) cs.rowCutPtr (i) -> print ();
311    for (int i=initColCuts; i<cs.sizeColCuts (); ++i) cs.colCutPtr (i) -> print ();
312  }
313
314  // else {
315  //     jnlst_ -> Printf (J_STRONGWARNING, J_COUENNE,
316  //                    "In-BB disjunctive cuts: %d row cuts, %d col cuts\n",
317  //                    cs.sizeRowCuts () - initRowCuts,
318  //                    cs.sizeColCuts () - initColCuts);
319  //   if ((info.level <= 0) && !(info.inTree))
320  //     jnlst_ -> Printf (J_ERROR, J_COUENNE,
321  //                    "infeasible node\n");
322  // }
323
324  delete csi;
325
326  septime_ += (CoinCpuTime () - time);
327}
Note: See TracBrowser for help on using the repository browser.