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

Last change on this file since 908 was 908, checked in by pbelotti, 8 years ago

clean-up of sdp cuts

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