source: stable/0.2/Couenne/src/problem/CouenneSolverInterface.cpp @ 159

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

created new stable branch 0.2 from trunk (rev. 157)

File size: 11.6 KB
Line 
1/* $Id: CouenneSolverInterface.cpp 154 2009-06-16 18:52:53Z pbelotti $ */
2/*
3 * Name:    CouenneSolverInterface.cpp
4 * Authors: Pietro Belotti, Carnegie Mellon University
5 *          Andreas Waechter, IBM Corp.
6 * Purpose: Implementation of the OsiSolverInterface::resolve () method
7 *
8 * (C) Carnegie-Mellon University, 2006-09.
9 * This file is licensed under the Common Public License (CPL)
10 */
11
12#include "OsiClpSolverInterface.hpp"
13#include "CouenneProblem.hpp"
14#include "CouenneSolverInterface.hpp"
15#include "CglTreeInfo.hpp"
16
17/// constructor
18CouenneSolverInterface::CouenneSolverInterface (CouenneCutGenerator *cg /*= NULL*/):
19
20  OsiClpSolverInterface(),
21  cutgen_ (cg),
22  knowInfeasible_(false),
23  knowOptimal_(false),
24  knowDualInfeasible_(false),
25  doingResolve_ (true) {
26
27  // prevents from running OsiClpSolverInterface::tightenBounds()
28  specialOptions_ = specialOptions_ | 262144; 
29}
30
31/// copy constructor
32CouenneSolverInterface::CouenneSolverInterface (const CouenneSolverInterface &src):
33
34  OsiSolverInterface    (src),
35  OsiClpSolverInterface (src),
36  cutgen_               (src.cutgen_),
37  knowInfeasible_       (src.knowInfeasible_),
38  knowOptimal_          (src.knowOptimal_),
39  knowDualInfeasible_   (src.knowDualInfeasible_),
40  doingResolve_         (src.doingResolve_) {}
41
42/// Destructor
43CouenneSolverInterface::~CouenneSolverInterface () {
44  //  if (cutgen_)
45  //    delete cutgen_;
46}
47
48
49/// Solve initial LP relaxation
50void CouenneSolverInterface::initialSolve () {
51  /*printf ("------------------------------------- INITIAL SOLVE\n");
52  for (int i=0; i<getNumCols(); i++)
53    printf ("%4d. %20.5g [%20.5g %20.5g]\n",
54            i, getColSolution () [i],
55            getColLower () [i],
56            getColUpper () [i]);
57
58            cutgen_ -> Problem () -> print ();*/
59
60  knowInfeasible_     = 
61  knowOptimal_        = 
62  knowDualInfeasible_ = false;
63
64  OsiClpSolverInterface::initialSolve ();
65  //writeLp ("initialLP");
66
67  if (getObjValue () <= - Couenne_large_bound)
68    knowDualInfeasible_ = true;
69
70  // some originals may be unused due to their zero multiplicity (that
71  // happens when they are duplicates), restore their value
72  if (cutgen_ -> Problem () -> nUnusedOriginals () > 0) {
73    CouNumber *x = new CouNumber [getNumCols ()];
74    CoinCopyN (getColSolution (), getNumCols (), x);
75    cutgen_ -> Problem () -> restoreUnusedOriginals (x);
76    setColSolution (x);
77    delete [] x;
78  }
79
80  /*
81    printf ("------------------------------------- INITSOLV\n");
82  for (int i=0; i<cutgen_ -> Problem () -> nOrigVars (); i++)
83    //if (cutgen_ -> Problem () -> Var (i) -> Multiplicity () <= 0)
84      {
85      printf ("%4d. %20.5g [%20.5g %20.5g]   ",
86              i, getColSolution () [i],
87              getColLower () [i],
88              getColUpper () [i]);
89      cutgen_ -> Problem () -> Var (i) -> print ();
90      printf ("\n");
91    }
92  */
93}
94
95bool CouenneSolverInterface::isProvenPrimalInfeasible() const {
96  return knowInfeasible_ || OsiClpSolverInterface::isProvenPrimalInfeasible();
97}
98
99bool CouenneSolverInterface::isProvenOptimal() const {
100  return knowOptimal_ || OsiClpSolverInterface::isProvenOptimal();
101}
102
103bool CouenneSolverInterface::isProvenDualInfeasible() const {
104  return knowDualInfeasible_ || OsiClpSolverInterface::isProvenDualInfeasible();
105}
106
107/// Defined in Couenne/src/convex/generateCuts.cpp
108void sparse2dense (int, t_chg_bounds *, int *&, int &);
109
110
111/// Resolve an LP relaxation after problem modification
112void CouenneSolverInterface::resolve () {
113  /*printf ("------------------------------------- RESOLVE\n");
114  for (int i=0; i<getNumCols(); i++)
115    printf ("%4d. %20.5g [%20.5g %20.5g]\n",
116            i, getColSolution () [i],
117            getColLower () [i],
118            getColUpper () [i]);*/
119
120  // CUT THIS! Some about-to-be-resolved problems have variables with
121  // lower +inf. That is, between CbcModel::initialSolve() and
122  // CbcModel::resolve(). I couldn't spot where in Couenne this
123  // happens.
124
125  ////////////////////////////////////// Cut {
126  /*const double
127    *lb = getColLower (),
128    *ub = getColUpper ();
129
130  for (int i=getNumCols(); i--;) {
131    if (lb [i] >  COUENNE_INFINITY)
132      setColLower (i, cutgen_ -> Problem () -> Lb (i));
133    //setColLower (i, -COIN_DBL_MAX);//cutgen_ -> Problem () -> Lb (i));
134    if (ub [i] < -COUENNE_INFINITY)
135      setColUpper (i, cutgen_ -> Problem () -> Ub (i));
136    //setColUpper (i,  COIN_DBL_MAX);//cutgen_ -> Problem () -> Ub (i));
137    }*/
138  ////////////////////////////////////// Cut }
139
140  static int count = -1;
141  char filename [30];
142
143  // save problem to be loaded later
144  if (cutgen_ && (cutgen_ -> check_lp ())) {
145    count++;
146    sprintf (filename, "resolve_%d", count);
147    writeMps (filename);
148  }
149
150  knowInfeasible_     = 
151  knowOptimal_        = 
152  knowDualInfeasible_ = false;
153
154  const CoinWarmStart *ws = NULL;
155
156  if (cutgen_ && (cutgen_ -> check_lp ()))
157    ws = getWarmStart ();
158
159  //deleteScaleFactors ();
160
161  // re-solve problem
162  OsiClpSolverInterface::resolve ();
163
164  if (getObjValue () <= - Couenne_large_bound)
165    knowDualInfeasible_ = true;
166
167  CouNumber objval = getObjValue (),
168    curCutoff = cutgen_ -> Problem () -> getCutOff ();
169
170  // check if resolve found new integer solution
171  if (doingResolve () &&                 // this is not called from strong branching
172      isProvenOptimal () &&
173      (objval < curCutoff - COUENNE_EPS) &&
174      (cutgen_ -> Problem () -> checkNLP (getColSolution (), objval, true)) &&
175      (objval < curCutoff - COUENNE_EPS) && // check again as it may have changed
176      (objval > -COUENNE_INFINITY/2)) {    // check if it makes sense
177
178    // also save the solution so that cbcModel::setBestSolution saves it too
179
180    //printf ("new cutoff from CSI: %g\n", objval);
181    cutgen_ -> Problem () -> setCutOff (objval);
182  }
183
184  // some originals may be unused due to their zero multiplicity (that
185  // happens when they are duplicates), restore their value
186  if (cutgen_ -> Problem () -> nUnusedOriginals () > 0) {
187    CouNumber *x = new CouNumber [getNumCols ()];
188    CoinCopyN (getColSolution (), getNumCols (), x);
189    cutgen_ -> Problem () -> restoreUnusedOriginals (x);
190    setColSolution (x);
191    delete [] x;
192  }
193
194  //cutgen_ -> Problem () -> restoreUnusedOriginals (this);
195
196  // check LP independently
197  if (cutgen_ && (cutgen_ -> check_lp ())) {
198
199    OsiSolverInterface
200      *nsi = new OsiClpSolverInterface,
201      *csi = clone ();
202
203    sprintf (filename, "resolve_%d.mps", count);
204    nsi -> readMps (filename);
205
206    nsi -> messageHandler () -> setLogLevel (0);
207    nsi -> setWarmStart (ws);
208
209    nsi -> initialSolve ();
210
211    if ((nsi -> isProvenOptimal () && isProvenOptimal ()) ||
212        !(nsi -> isProvenOptimal ()) && !isProvenOptimal ()) {
213
214      if (nsi -> isProvenOptimal () &&
215          (fabs (nsi -> getObjValue () - getObjValue ()) / 
216           (1. + fabs (nsi -> getObjValue ()) + fabs (getObjValue ())) > 1e-2))
217
218        printf ("Warning: discrepancy between saved %g and current %g [%g], file %s\n", 
219                nsi -> getObjValue (),  getObjValue (),
220                nsi -> getObjValue () - getObjValue (),
221                filename);
222    }
223
224    csi -> messageHandler () -> setLogLevel (0);
225    csi -> setWarmStart (ws);
226
227    csi -> initialSolve ();
228
229    if ((csi -> isProvenOptimal () && isProvenOptimal ()) ||
230        !(csi -> isProvenOptimal ()) && !isProvenOptimal ()) {
231
232      if (csi -> isProvenOptimal () &&
233          (fabs (csi -> getObjValue () - getObjValue ()) / 
234           (1. + fabs (csi -> getObjValue ()) + fabs (getObjValue ())) > 1e-2))
235
236        printf ("Warning: discrepancy between cloned %g and current %g [%g]\n", 
237                csi -> getObjValue (),  getObjValue (),
238                csi -> getObjValue () - getObjValue ());
239    }
240
241    delete nsi;
242    delete csi;
243   
244    delete ws;
245
246    //else printf ("Warning: discrepancy between statuses %s -- %s feasible\n",
247    //filename, isProvenOptimal () ? "current" : "saved");
248  }
249}
250
251
252/// Create a hot start snapshot of the optimization process.
253void CouenneSolverInterface::markHotStart () {
254  //printf(">>>> markHotStart\n");
255  // Using OsiClpSolverInterface doesn't work yet...
256  OsiSolverInterface::markHotStart ();
257}
258
259
260/// Delete the hot start snapshot.
261void CouenneSolverInterface::unmarkHotStart() {
262  //printf("<<<< unmarkHotStart\n");
263  OsiSolverInterface::unmarkHotStart();
264}
265
266
267
268/// Optimize starting from the hot start snapshot.
269void CouenneSolverInterface::solveFromHotStart() {
270
271  //OsiClpSolverInterface::solveFromHotStart ();
272
273  //#if 0
274  knowInfeasible_     = 
275  knowOptimal_        = 
276  knowDualInfeasible_ = false;
277
278  /*
279  const int ncols = cutgen_ -> Problem () -> nVars ();
280
281  cutgen_ -> Problem () -> domain () -> push
282    (cutgen_ -> Problem () -> nVars (),
283     getColSolution (),
284     getColLower    (),
285     getColUpper    ());
286
287  // This vector contains variables whose bounds have changed due to
288  // branching, reduced cost fixing, or bound tightening below. To be
289  // used with malloc/realloc/free
290
291  t_chg_bounds *chg_bds = new t_chg_bounds [ncols];
292
293  OsiCuts cs;
294
295  Bonmin::BabInfo *babInfo = dynamic_cast <Bonmin::BabInfo *> (getAuxiliaryInfo ());
296
297  if (cutgen_ -> Problem () -> doFBBT () &&
298      (!(cutgen_ -> Problem () -> boundTightening (chg_bds, babInfo)))) {
299
300#ifdef DEBUG
301    printf ("#### BT says infeasible before re-solve\n");
302#endif
303    knowInfeasible_ = true;
304    cutgen_ -> Problem () -> domain () -> pop ();
305    return;
306  }
307
308  int *changed = NULL, nchanged;
309  sparse2dense (ncols, chg_bds, changed, nchanged);
310
311  // change tightened bounds through OsiCuts
312  if (nchanged)
313    cutgen_ -> genColCuts (*this, cs, nchanged, changed);
314
315  const int nRowsBeforeRowCuts = getNumRows();
316  //printf("NumRows before getRowCuts = %d\n", getNumRows());
317  cutgen_ -> genRowCuts (*this, cs, nchanged, changed, //CglTreeInfo(),
318                         chg_bds, false);
319
320  // Now go through the list of cuts and apply the column cuts
321  // directly as changes on bounds
322  while(cs.sizeColCuts()) {
323    const OsiColCut& ccut = cs.colCut(0);
324    const CoinPackedVector& lbs = ccut.lbs();
325    int nele = lbs.getNumElements();
326    const int* idxs = lbs.getIndices();
327    const double* eles = lbs.getElements();
328    const double* bnds = getColLower();
329    for (int i=0; i<nele; i++) {
330      if (bnds[*idxs] < *eles) {
331        //printf ("setcolLower %d %g\n", *idxs, *eles);
332        setColLower(*idxs,*eles);
333      }
334      idxs++;
335      eles++;
336    }
337    const CoinPackedVector& ubs = ccut.ubs();
338    nele = ubs.getNumElements();
339    idxs = ubs.getIndices();
340    eles = ubs.getElements();
341    bnds = getColUpper();
342    for (int i=0; i<nele; i++) {
343      if (bnds[*idxs] > *eles) {
344        //printf ("setcolUpper %d %g\n", *idxs, *eles);
345        setColUpper(*idxs,*eles);
346      }
347      idxs++;
348      eles++;
349    }
350    cs.eraseColCut(0);
351  }
352
353  applyCuts (cs);
354
355  const int nRowsAfterRowCuts = getNumRows();
356  //printf("NumRows after applyCuts = %d\n", getNumRows());
357  */
358
359  resolve();
360
361  if (getObjValue () <= - Couenne_large_bound)
362    knowDualInfeasible_ = true;
363
364  // some originals may be unused due to their zero multiplicity (that
365  // happens when they are duplicates), restore their value
366  if (cutgen_ -> Problem () -> nUnusedOriginals () > 0) {
367    CouNumber *x = new CouNumber [getNumCols ()];
368    CoinCopyN (getColSolution (), getNumCols (), x);
369    cutgen_ -> Problem () -> restoreUnusedOriginals (x);
370    setColSolution (x);
371    delete [] x;
372  }
373
374  if (isProvenPrimalInfeasible ()) knowInfeasible_     = true;
375  if (isProvenOptimal          ()) knowOptimal_        = true;
376  if (isProvenDualInfeasible   ()) knowDualInfeasible_ = true;
377
378  //printf("obj value = %e\n",getObjValue());
379
380  // now undo the row cuts
381  /*
382  int nrowsdel = nRowsAfterRowCuts-nRowsBeforeRowCuts;
383  int* rowsdel = new int[nrowsdel];
384  for(int i=0; i<nrowsdel; i++) {
385    rowsdel[i] = nRowsBeforeRowCuts+i;
386  }
387  deleteRows(nrowsdel, rowsdel);
388  delete [] rowsdel;
389  //printf("NumRows after deleting = %d\n", getNumRows());
390
391  cutgen_ -> Problem () -> domain () -> pop ();
392  */
393  //#endif
394}
Note: See TracBrowser for help on using the repository browser.