source: trunk/Couenne/src/problem/CouenneSolverInterface.cpp @ 112

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

isolate restoreUnused to when there is at least one unused -- takes care of reformulated stockcycle problem (thanks to C. D'Ambrosio for the bug report). Check linear terms before adding a full exprGroup -- need a constructor pilot. Quicker solution feasibility check (with initial integrality check).

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