source: trunk/Couenne/src/bound_tightening/obbt.cpp @ 318

Last change on this file since 318 was 318, checked in by pbelotti, 12 years ago

fixed bug that gave 2.982 on minlplib/alan -- wrong OBBT implication on univariate (semi)-aux

  • Property svn:keywords set to Id
File size: 6.6 KB
Line 
1/* $Id: obbt.cpp 318 2010-03-16 03:22:58Z pbelotti $
2 *
3 * Name:    obbt.cpp
4 * Author:  Pietro Belotti
5 * Purpose: Optimality-Based Bound Tightening
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CoinHelperFunctions.hpp"
12#include "CglCutGenerator.hpp"
13
14#include "CouenneExpression.hpp"
15#include "CouenneCutGenerator.hpp"
16#include "CouenneProblem.hpp"
17#include "OsiClpSolverInterface.hpp"
18
19#define OBBT_EPS 1e-3
20#define MAX_OBBT_LP_ITERATION 100
21
22// minimum #bound changed in obbt to generate further cuts
23#define THRES_NBD_CHANGED 1
24
25// maximum number of obbt iterations
26#define MAX_OBBT_ITER 1
27
28// defined in generateCuts.cpp
29void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged);
30
31
32// OBBT for one sense (max/min) and one class of variables (orig/aux)
33int CouenneProblem::call_iter (OsiSolverInterface *csi, 
34                               t_chg_bounds *chg_bds, 
35                               const CoinWarmStart *warmstart, 
36                               Bonmin::BabInfo *babInfo,
37                               double *objcoe,
38                               enum nodeType type,
39                               int sense) const {
40
41  int ncols   = csi -> getNumCols (),
42      nimprov = 0;
43
44  for (int ii=0; ii<ncols; ii++) {
45
46    if (CoinCpuTime () > maxCpuTime_)
47      break;
48
49    int i = evalOrder (ii);
50
51    enum expression::auxSign aSign = Var (i) -> sign ();
52
53    if ((Var (i) -> Type () == type)     &&
54        (Var (i) -> Multiplicity () > 0) &&
55        ((type == VAR)                               || 
56         (aSign  == expression::EQ) ||
57         ((aSign == expression::LEQ) && (sense > 0)) ||
58         ((aSign == expression::GEQ) && (sense < 0)))) {
59
60      int ni = obbt_iter (csi, chg_bds, warmstart, babInfo, objcoe, sense, i);
61
62//       {
63//      // ToDo: Pipe all output through journalist
64//      Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
65//                      "  bounds after obbt step  =====================\n  ");
66//      int j=0;
67//      for (int i=0; i < nVars (); i++)
68//        if (variables_ [i] -> Multiplicity () > 0) {
69//          Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,
70//                          "x_%03d [%+10g %+10g] ", i,
71//                          domain_. lb (i),
72//                          domain_. ub (i));
73//          if (!(++j % 6)) Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n  ");
74//        }
75//      if (j % 6) Jnlst()->Printf(J_VECTOR, J_BOUNDTIGHTENING,"\n");
76//       }
77
78      if (ni < 0) return ni;
79      nimprov += ni;
80    }
81  }
82
83  return nimprov;
84}
85
86
87/// Optimality based bound tightening -- inner loop
88
89int CouenneProblem::obbtInner (OsiSolverInterface *csi,
90                               OsiCuts &cs,
91                               t_chg_bounds *chg_bds,
92                               Bonmin::BabInfo * babInfo) const {
93
94  // set large bounds to infinity (as per suggestion by JJF)
95
96  int ncols = csi -> getNumCols ();
97  const double *lb = csi -> getColLower (),
98               *ub = csi -> getColUpper ();
99
100  double inf = csi -> getInfinity ();
101
102  for (int i=ncols; i--;) {
103    if (lb [i] < - COUENNE_INFINITY) csi -> setColLower (i, -inf);
104    if (ub [i] >   COUENNE_INFINITY) csi -> setColUpper (i,  inf);
105  }
106
107  //  csi -> setHintParam (OsiDoDualInResolve, false);
108
109  // setup cloned interface for later use
110  csi -> setObjSense (1); // minimization
111  csi -> setIntParam (OsiMaxNumIteration, MAX_OBBT_LP_ITERATION);
112  csi -> applyCuts (cs);   // apply all (row+column) cuts to formulation
113  csi -> initialSolve ();
114
115  const CoinWarmStart *warmstart = csi -> getWarmStart ();
116
117  // improve each bound
118
119  double *objcoe = (double *) malloc (ncols * sizeof (double));
120
121  // set obj function coefficients to zero
122  for (int i=ncols; i--;)
123    *objcoe++ = 0.;
124  objcoe -= ncols;
125
126  csi -> setObjective (objcoe);
127  csi -> setObjSense (1);        // minimization
128
129  int nimprov = 0;
130 
131  const int Infeasible = 1;
132
133  try {
134
135    int ni;
136
137    if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, VAR,  1)) < 0) throw Infeasible;
138    nimprov += ni;
139
140    if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, VAR, -1)) < 0) throw Infeasible;
141    nimprov += ni;
142
143    if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, AUX,  1)) < 0) throw Infeasible;
144    nimprov += ni;
145
146    if ((ni = call_iter (csi, chg_bds, warmstart, babInfo, objcoe, AUX, -1)) < 0) throw Infeasible;
147    nimprov += ni;
148  }
149
150  catch (int exception) {
151
152    if (exception == Infeasible)
153      nimprov = -1;
154  }
155
156  free (objcoe);
157  delete warmstart;
158
159  return (nimprov);
160}
161
162
163// Optimality based bound tightening -- main loop
164int CouenneProblem::obbt (const CouenneCutGenerator *cg,
165                          const OsiSolverInterface &si,
166                          OsiCuts &cs,
167                          const CglTreeInfo &info,
168                          Bonmin::BabInfo * babInfo,
169                          t_chg_bounds *chg_bds) {
170
171  // Do OBBT if:
172
173  if (doOBBT_ &&                        // flag is checked, AND
174      ((logObbtLev_ != 0) ||               // (parameter is nonzero OR
175       (info.level == 0)) &&               //  we are at root node), AND
176      (info.pass == 0) &&               // at first round of cuts, AND
177      ((logObbtLev_ < 0) ||               // (logObbtLev = -1, OR
178       (info.level <= logObbtLev_) ||     //  depth is lower than COU_OBBT_CUTOFF_LEVEL, OR
179                                          //  probability inversely proportional to the level)
180       (CoinDrand48 () < pow (2., (double) logObbtLev_ - (info.level + 1))))) {
181
182    jnlst_ -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "----- OBBT\n");
183
184    // TODO: why check info.pass==0? Why not more than one pass? It
185    // should be anyway checked that info.level be >= 0 as <0 means
186    // first call at root node
187
188    OsiSolverInterface *csi = si.clone (true);
189    //dynamic_cast <CouenneSolverInterface<T> *>
190
191    OsiClpSolverInterface *clpcsi = dynamic_cast <OsiClpSolverInterface *> (csi);
192
193    if (clpcsi)
194      clpcsi -> setupForRepeatedUse ();
195    //csi -> doingResolve () = false;
196
197    //csi -> setHintParam (OsiDoDualInResolve, false);
198
199    int nImprov, nIter = 0;
200
201    bool notImproved = false;
202
203    while (!notImproved && 
204           (nIter++ < MAX_OBBT_ITER) &&
205           ((nImprov = obbtInner (csi, cs, chg_bds, babInfo)) > 0) &&
206           (CoinCpuTime () < maxCpuTime_)) {
207
208      int nchanged, *changed = NULL;
209
210      /// OBBT has tightened, add improved bounds
211      sparse2dense (nVars (), chg_bds, changed, nchanged);
212      cg -> genColCuts (*csi, cs, nchanged, changed);
213
214      if ((nIter < MAX_OBBT_ITER) && 
215          (nImprov >= THRES_NBD_CHANGED)) {
216
217        // only generate new row cuts if improvents are enough
218        int nCurCuts = cs.sizeRowCuts ();
219        cg -> genRowCuts (*csi, cs, nchanged, changed, chg_bds);
220
221        if (nCurCuts == cs.sizeRowCuts ())
222          notImproved = true; // repeat only if new cuts available
223
224      } else notImproved = true;
225
226      if (changed) 
227        free (changed);
228    }
229
230    //csi -> doingResolve () = true;
231
232    delete csi;
233
234    if (nImprov < 0) {
235      jnlst_->Printf(J_ITERSUMMARY, J_BOUNDTIGHTENING, "  Couenne: infeasible node after OBBT\n");
236      return -1;
237    }
238  }
239
240  return 0;
241}
Note: See TracBrowser for help on using the repository browser.