source: stable/0.4/Couenne/src/disjunctive/generateDisjCuts.cpp @ 946

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

merge r945 from trunk: allow compiling against Cgl 0.58

  • Property svn:keywords set to Author Date Id Revision
File size: 9.8 KB
Line 
1/* $Id: generateDisjCuts.cpp 946 2013-04-15 22:20:38Z 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    bool tightened = false;
266
267    for (int i=0; i<ncols; i++, newLo++, newUp++) {
268
269      if (*newLo > *oldLo++ + COUENNE_EPS) {tighterLower.insert (i, *newLo); tightened = true;}
270      if (*newUp < *oldUp++ - COUENNE_EPS) {tighterUpper.insert (i, *newUp); tightened = true;}
271    }
272
273    if (tightened) {
274      OsiColCut tighter;
275      tighter.setLbs (tighterLower);
276      tighter.setUbs (tighterUpper);
277      if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
278        printf ("tightened bounds in disjunctive cuts:");
279        tighter.print ();
280      }
281      cs.insert (tighter);
282    }
283
284    int deltaNcuts = 
285      cs.sizeRowCuts () - initRowCuts + 
286      cs.sizeColCuts () - initColCuts;
287
288    if (info.level <= 0 && !(info.inTree)) 
289      nrootcuts_ += deltaNcuts;
290    ntotalcuts_ += deltaNcuts;
291
292    if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
293
294      if (cs.sizeRowCuts()>initRowCuts) printf ("added %d row cuts\n", cs.sizeRowCuts () - initRowCuts);
295      if (cs.sizeColCuts()>initColCuts) printf ("added %d col cuts\n", cs.sizeColCuts () - initColCuts);
296    }
297
298    if ((info.level <= 0) && !(info.inTree))
299      jnlst_ -> Printf (J_ERROR, J_COUENNE, 
300                        "%d cuts\n", CoinMax (0, nrootcuts_ - nInitCuts));
301    else if (deltaNcuts)
302      jnlst_ -> Printf (J_WARNING, J_COUENNE, 
303                        "In-BB disjunctive cuts: %d row cuts, %d col cuts\n",
304                        cs.sizeRowCuts () - initRowCuts,
305                        cs.sizeColCuts () - initColCuts);
306  }
307
308  // else {
309  //     jnlst_ -> Printf (J_STRONGWARNING, J_COUENNE,
310  //                    "In-BB disjunctive cuts: %d row cuts, %d col cuts\n",
311  //                    cs.sizeRowCuts () - initRowCuts,
312  //                    cs.sizeColCuts () - initColCuts);
313  //   if ((info.level <= 0) && !(info.inTree))
314  //     jnlst_ -> Printf (J_ERROR, J_COUENNE,
315  //                    "infeasible node\n");
316  // }
317
318  delete csi;
319
320  septime_ += (CoinCpuTime () - time);
321}
Note: See TracBrowser for help on using the repository browser.