CouenneFPscipSolve.cpp
Go to the documentation of this file.
1 /* $Id: CouenneFPscipSolve.cpp 1075 2014-04-11 23:38:43Z pbelotti $
2  *
3  * Name: CouenneFPscipSolve.cpp
4  * Authors: Timo Berthold, ZIB Berlin
5  * Purpose: call SCIP
6  *
7  * This file is licensed under the Eclipse Public License (EPL)
8  */
9 
10 #include "CouenneProblem.hpp"
11 #include "CouenneFeasPump.hpp"
12 #include "CouenneExprVar.hpp"
13 
14 #include "cons_rowcuts.h"
15 
16 using namespace Couenne;
17 
18 #ifdef COIN_HAS_SCIP
19 
20 /* general SCIP includes */
21 #include "scip/scip.h"
22 #include "scip/cons_linear.h"
23 #include "scip/scipdefplugins.h"
24 #include "scip/cons_bounddisjunction.h"
25 
26 SCIP_RETCODE CouenneFeasPump::ScipSolve (const double *nSol, double* &sol, int niter, int* nsuciter, CouNumber &obj) {
27 
28  static int currentmilpmethod = 0;
29 
30  int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
31 
32  SCIP* scip;
33 
34  SCIP_VAR** vars;
35  SCIP_Real* lbs;
36  SCIP_Real* ubs;
37  const SCIP_Real* objs;
38  const char* vartypes;
39  const CoinPackedMatrix * matrix;
40  const CoinBigIndex* rowstarts;
41  const int* rowlengths;
42  const SCIP_Real* coeffs;
43  const SCIP_Real* lhss;
44  const SCIP_Real* rhss;
45  const int* indices;
46 
47  SCIP_Real timelimit;
48 
49  SCIP_VAR **tabuvars;
50  SCIP_Real *tabubounds;
51  SCIP_BOUNDTYPE *tabuboundtypes;
52 
53  double infinity;
54  int nvars;
55  int nconss;
56  int nscipsols;
57 
58  bool solveagain;
59 
60  // COUENNE_INFINITY , getInfinity()
61 
62  // get problem data
63  nvars = milp_ -> getNumCols ();
64  nconss = milp_ -> getNumRows ();
65  infinity = milp_ -> getInfinity ();
66 
67  // get variable data
68  lbs = CoinCopyOfArray (milp_ -> getColLower (), milp_ -> getNumCols ());
69  ubs = CoinCopyOfArray (milp_ -> getColUpper (), milp_ -> getNumCols ());
70 
71  objs = milp_ -> getObjCoefficients ();
72  vartypes = milp_ -> getColType ();
73 
74  // tighten bounds with data from problem_ (which may have run btCore())
75 
76  for (int i = 0; i < problem_->nVars(); ++i)
77  if (problem_ -> Var (i) -> Multiplicity () > 0) {
78  if (problem_ -> Lb (i) > lbs [i]) lbs [i] = problem_ -> Lb (i);
79  if (problem_ -> Ub (i) < ubs [i]) ubs [i] = problem_ -> Ub (i);
80  }
81 
82  // get row data
83  lhss = milp_ -> getRowLower ();
84  rhss = milp_ -> getRowUpper ();
85 
86  // get matrix data
87  matrix = milp_ -> getMatrixByRow();
88  rowstarts = matrix -> getVectorStarts();
89  rowlengths = matrix -> getVectorLengths();
90  coeffs = matrix -> getElements();
91  indices = matrix -> getIndices();
92 
93  // access tabuPool_ in this class
94 
95  if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
96  SCIPdebugMessage("create SCIP problem instance with %d variables and %d constraints.\n", nvars, nconss);
97  }
98 
99  // initialize SCIP
100  SCIP_CALL( SCIPcreate (&scip) );
101  assert(scip != NULL);
102 
103  // include default SCIP plugins
104  SCIP_CALL( SCIPincludeDefaultPlugins(scip) );
105 
106  // include row cut constraint hanlder
107  if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
108  {
109  SCIP_CALL( SCIPincludeConshdlrRowcuts(scip, couenneCG_, milp_) );
110  }
111 
112  // create problem instance in SCIP
113  SCIP_CALL( SCIPcreateProb(scip, "auxiliary FeasPump MILP", NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
114 
115  // allocate local memory for SCIP variable array
116  SCIP_CALL( SCIPallocMemoryArray(scip, &vars, nvars) );
117 
118  // Allocating space for tabu stuff
119 
120  SCIP_CALL( SCIPallocMemoryArray(scip, &tabuvars , 2*nvars) ); // fix to nvars + nIntvars
121  SCIP_CALL( SCIPallocMemoryArray(scip, &tabubounds , 2*nvars) );
122  SCIP_CALL( SCIPallocMemoryArray(scip, &tabuboundtypes, 2*nvars) );
123 
124  // one variable for objective !!!!!!!!!
125 
126  // create variables
127  for (int i=0; i<nvars; i++) {
128 
129  char varname[SCIP_MAXSTRLEN];
130 
131  bool neglect =
132  (i < problem_ -> nVars ()) &&
133  (problem_ -> Var (i) -> Multiplicity () <= 0);
134 
135  // check that all data is in valid ranges
136  assert( 0 <= vartypes[i] && vartypes[i] <= 2);
137 
138  if (!neglect) {
139 
140  checkInfinity(scip, lbs[i], infinity);
141  checkInfinity(scip, ubs[i], infinity);
142  }
143 
144  // all variables are named x_i
145  (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "x_%d", i);
146  SCIP_CALL( SCIPcreateVar(scip, &vars[i], varname,
147  CoinMin (lbs [i], ubs [i]),
148  CoinMax (lbs [i], ubs [i]),
149  objs [i],
150  vartypes[i] == 0 ? SCIP_VARTYPE_CONTINUOUS : (vartypes[i] == 1 ? SCIP_VARTYPE_BINARY : SCIP_VARTYPE_INTEGER),
151  TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
152 
153  // add the variable to SCIP
154  SCIP_CALL( SCIPaddVar(scip, vars[i]) );
155  }
156 
157  for (std::set <CouenneFPsolution, compareSol>:: iterator i = tabuPool_ . begin ();
158  i != tabuPool_ . end (); ++i) {
159 
160  const double *x = i -> x ();
161 
162  int nEntries = 0;
163 
164  SCIP_CONS *tabucons = NULL;
165 
166  for (int j = 0; j < i -> n (); ++j) {
167 
168  if (problem_ -> Var (j) -> isInteger () &&
169  (problem_ -> Var (j) -> Multiplicity () > 0) &&
170  (fabs (ubs [j] - lbs [j]) > .5)) {
171 
172  // problem_ -> Ub (j) - problem_ -> Lb (j) > .5) {
173 
174  // if (fabs (x [j] - floor (x [j] + .5)) >= SCIPfeastol (scip)) {
175  // printf ("integer var x%d not really integer: %e\n", j, x [j]);
176  // }
177 
178  //assert (fabs (lbs [j] - problem_ -> Lb (j)) < SCIPfeastol (scip));
179  //assert (fabs (ubs [j] - problem_ -> Ub (j)) < SCIPfeastol (scip));
180  assert (fabs (x [j] - floor (x [j] + .5)) < SCIPfeastol (scip) * 1.e3);
181 
182  assert (nEntries <= 2*nvars - 2);
183 
184  double x_rounded = floor (x [j] + .5);
185 
186  if (x [j] >= ubs [j] - COUENNE_EPS) {
187 
188  tabuvars [nEntries] = vars [j];
189  tabubounds [nEntries] = x_rounded - 1.;
190  tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_UPPER;
191 
192  } else if (x [j] <= lbs [j] + COUENNE_EPS) {
193 
194  tabuvars [nEntries] = vars [j];
195  tabubounds [nEntries] = x_rounded + 1.;
196  tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_LOWER;
197 
198  } else {
199 
200  tabuvars [nEntries] = vars [j];
201  tabubounds [nEntries] = x_rounded - 1.;
202  tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_UPPER;
203 
204  tabuvars [nEntries] = vars [j];
205  tabubounds [nEntries] = x_rounded + 1.;
206  tabuboundtypes [nEntries++] = SCIP_BOUNDTYPE_LOWER;
207  }
208  }
209  }
210 
211  if (nEntries != 0) {
212 
213  SCIP_CALL (SCIPcreateConsBounddisjunction (scip, &tabucons, "Tabu Solution", nEntries,
214  tabuvars, tabuboundtypes, tabubounds,
215  TRUE, TRUE, TRUE, TRUE, TRUE,
216  FALSE, FALSE, FALSE, FALSE, FALSE));
217 
218  SCIP_CALL( SCIPaddCons(scip, tabucons) );
219  //SCIP_CALL (SCIPprintCons (scip, tabucons, NULL));
220  SCIP_CALL (SCIPreleaseCons (scip, &tabucons));
221  }
222  }
223 
224  // create constraints
225  for (int i=0; i<nconss; i++) {
226 
227  SCIP_CONS* cons;
228 
229  char consname[SCIP_MAXSTRLEN];
230  (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN, "row_%d", i);
231 
232  // check that all data is in valid ranges
233  checkInfinity(scip, lhss[i], infinity);
234  checkInfinity(scip, rhss[i], infinity);
235 
236  // create an empty linear constraint
237  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, consname, 0, NULL, NULL, lhss[i], rhss[i],
238  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
239 
240  // add variables to constraint
241  for(int j=rowstarts[i]; j<rowstarts[i]+rowlengths[i]; j++)
242  {
243  checkInfinity(scip, coeffs[j], infinity);
244  SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[indices[j]], coeffs[j]) );
245  }
246 
247  // add constraint to SCIP
248  SCIP_CALL( SCIPaddCons(scip, cons) );
249  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
250  }
251 
252  // determine the method to solve the MILP
253  if (milpMethod_ == 0 && niter == 0)
254  {
255  // initialize currentmilpmethod: at the root node we solve the MILP with SCIP default
256  // deeper in the tree, we will solve the MILP by a FP heuristic
257  if( depth == 0 )
258  currentmilpmethod = 2;
259  else
260  currentmilpmethod = 4;
261  }
262  else if (milpMethod_ != 0)
263  currentmilpmethod = milpMethod_; // use a fixed method to solve the MILP
264 
265  // MILP solving loop. If the MILP terminates without a solution, it might get resolved with a more expensive atrategy
266  do {
267  solveagain = false;
268 
269  // reset parameters if MILP is solved agian
270  SCIP_CALL( SCIPresetParams(scip) );
271 
272  // deactivate SCIP output
273  if (!(problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_WARNING, J_NLPHEURISTIC))) {
274  SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 0) );
275  }
276 
277  // do not abort subproblem on CTRL-C
278  SCIP_CALL( SCIPsetBoolParam(scip, "misc/catchctrlc", FALSE) );
279 
280  // set time limit
281  timelimit = problem_ -> getMaxCpuTime () - CoinCpuTime ();
282 
283  if (timelimit < 0)
284  break;
285 
286  SCIP_CALL( SCIPsetRealParam(scip, "limits/time", timelimit) );
287 
288  if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_WARNING, J_NLPHEURISTIC)) {
289  SCIPinfoMessage(scip, NULL, "using MILP method: %d\n",currentmilpmethod);
290  }
291 
292  // tune SCIP differently, depending on the chosen method to solve the MILP
299  switch(currentmilpmethod)
300  {
301  case -1: // solve the MILP completely. SCIP's default setting should be best for this
302  if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
303  {
304  SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
305  }
306 
307  // As this is expensive, stop as soon as a solution is found
308 
309  SCIP_CALL( SCIPsetIntParam(scip, "limits/bestsol", 1) );
310  break;
311 
312  case 1: // Be aggressive in finding feasible solutions, but lazy about the dual bound.
313  // Set limits on overall nodes and stall nodes (nodes without incumbent improvement).
314  SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 1000) );
315  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 10000) );
316  SCIP_CALL( SCIPsetRealParam (scip, "limits/gap", .001) );
317 
318  // disable expensive cutting plane separation
319  //SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_FAST, TRUE) );
320 
321  // disable expensive presolving
322  //SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
323 
324  // use aggressive primal heuristics
325  SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) );
326 
327  // use best estimate node selection
328  /*
329  if( SCIPfindNodesel(scip, "estimate") != NULL )
330  {
331  SCIP_CALL( SCIPsetIntParam(scip, "nodeselection/estimate/stdpriority", INT_MAX/4) );
332  }
333  */
334 
335  // use inference branching
336  /*
337  if( SCIPfindBranchrule(scip, "inference") != NULL )
338  {
339  SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
340  }
341  */
342 
343  // disable conflict analysis
344  /*
345  SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useprop", FALSE) );
346  SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useinflp", FALSE) );
347  SCIP_CALL( SCIPsetBoolParam(scip, "conflict/useboundlp", FALSE) );
348  SCIP_CALL( SCIPsetBoolParam(scip, "conflict/usesb", FALSE) );
349  SCIP_CALL( SCIPsetBoolParam(scip, "conflict/usepseudo", FALSE) );
350  */
351 
352  break;
353 
354  case 2: // default SCIP with node limits
355  SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
356  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
357  SCIP_CALL( SCIPsetRealParam (scip, "limits/gap", .005) );
358 
359  // disable expensive dual techniques
360 
361  SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_FAST, TRUE) );
362  SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
363  SCIP_CALL( SCIPsetIntParam(scip, "branching/pscost/priority", INT_MAX/4) );
364 
365  break;
366  case 3: // solve the MILP with RENS. Disable most other features, enable RENS
367  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 1) );
368 
369  // disable cutting plane separation
370  SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
371 
372  // disable expensive presolving
373  SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
374 
375  // besides RENS, only use cheap heuristics
376  SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
377 
378  // use inference branching
379  if( SCIPfindBranchrule(scip, "inference") != NULL )
380  {
381  SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
382  }
383 
384  // ensure that RENS is called
385  if( SCIPfindHeur(scip, "rens") != NULL )
386  {
387  SCIP_CALL( SCIPsetIntParam(scip, "heuristics/rens/freq", 0) );
388  SCIP_CALL( SCIPsetRealParam(scip, "heuristics/rens/minfixingrate", 0.0) );
389  }
390  if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
391  {
392  SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
393  }
394  break;
395 
396  case 4: // solve the MILP with Feasibility Pump. Disable most other features, enable stage 3 for feaspump
397  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 1) );
398 
399  // disable cutting plane separation
400  SCIP_CALL( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
401 
402  // disable expensive presolving
403  SCIP_CALL( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
404 
405  // besides feaspump, only use cheap heuristics
406  SCIP_CALL( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
407 
408  // use inference branching
409  if( SCIPfindBranchrule(scip, "inference") != NULL )
410  {
411  SCIP_CALL( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
412  }
413 
414  // ensure that feasibility pump is called
415  if( SCIPfindHeur(scip, "feaspump") != NULL )
416  {
417  SCIP_CALL( SCIPsetIntParam(scip, "heuristics/feaspump/freq", 0) );
418  SCIP_CALL( SCIPsetIntParam(scip, "heuristics/feaspump/maxsols", -1) );
419 
420  // Comment next 8 lines as Workaround for mac/conv/batch
421  // to avoid assertion in stage 3 on nsolutions==0
422 
423  if( SCIPversion() <= 2.01 )
424  {
425  SCIP_CALL( SCIPsetBoolParam(scip, "heuristics/feaspump2/stage3", TRUE) );
426  }
427  else
428  {
429  SCIP_CALL( SCIPsetBoolParam(scip, "heuristics/feaspump/stage3", TRUE) );
430  }
431  }
432  if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
433  {
434  SCIP_CALL( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
435  }
436  break;
437 
438  case 5: // use rounding heuristic in Couenne
439 
440  SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
441  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
442  SCIP_CALL( SCIPsetRealParam (scip, "limits/gap", .05) );
443 
444  {
445  if (!sol)
446  sol = new CouNumber [problem_ -> nVars ()];
447 
448  if (!nSol)
449  nSol = milp_ -> getColSolution ();
450 
451  // whatever nSol, pick solution through getIntCand
452 
453  problem_ -> getIntegerCandidate (nSol, sol, problem_ -> Lb (), problem_ -> Ub ());
454 
455  SCIP_SOL* scipSol;
456  SCIP_CALL( SCIPcreateSol(scip, &scipSol, NULL ));
457  SCIP_CALL(SCIPsetSolVals(scip, scipSol, nvars, vars,(double*) sol));
458  }
459 
460  break;
461 
462  case 6: // round; TODO: but perturb first if we are cycling
463 
464  SCIP_CALL( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
465  SCIP_CALL( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
466  SCIP_CALL( SCIPsetRealParam (scip, "limits/gap", .05) );
467 
468  if (!sol)
469  sol = new CouNumber [problem_ -> nVars ()];
470 
471  if (!nSol)
472  nSol = milp_ -> getColSolution ();
473 
474  for (int j = 0; j < problem_ -> nVars (); ++j)
475  sol [j] =
476  (problem_ -> Var (j) -> Multiplicity () <= 0) ? 0. :
477  (problem_ -> Var (j) -> isInteger ()) ? COUENNE_round (nSol [j]) : nSol [j];
478 
479  SCIP_SOL* scipSol;
480  SCIP_CALL( SCIPcreateSol(scip, &scipSol, NULL ));
481  SCIP_CALL(SCIPsetSolVals(scip, scipSol, nvars, vars, (double*) sol));
482 
483  break;
484 
485 
486  default:
487  printf("Invalid MILP method in feasibility pump: %d\n", currentmilpmethod);
488  assert(false);
489  break;
490  }
491 
493 
494  SCIP_CALL(SCIPtransformProb(scip));
495 
496  {
497  SCIP_CONSHDLR* conshdlr_bounddisjunction = SCIPfindConshdlr(scip, "bounddisjunction");
498  int nbdconss= SCIPconshdlrGetNConss(conshdlr_bounddisjunction);
499  SCIP_CONS** bdconss = SCIPconshdlrGetConss(conshdlr_bounddisjunction);
500 
501  for(int i=0; i < nbdconss; i++)
502  {
503  SCIP_CONS* checkcons = bdconss[i];
504  SCIP_VAR** checkvars = SCIPgetVarsBounddisjunction(scip, checkcons);
505  SCIP_Real* checkbounds = SCIPgetBoundsBounddisjunction(scip,checkcons);
506  SCIP_BOUNDTYPE* checkboundtypes =SCIPgetBoundtypesBounddisjunction(scip, checkcons);
507  int checknvars = SCIPgetNVarsBounddisjunction(scip, checkcons);
508 
509  for(int j=i+1; j < nbdconss; j++)
510  {
511  SCIP_CONS* tmpcons =SCIPconshdlrGetConss(conshdlr_bounddisjunction)[j];
512  SCIP_VAR **tmpvars = SCIPgetVarsBounddisjunction(scip,tmpcons);
513  SCIP_Real *tmpbounds =SCIPgetBoundsBounddisjunction(scip, tmpcons);
514  SCIP_BOUNDTYPE *tmpboundtypes =SCIPgetBoundtypesBounddisjunction(scip, tmpcons);
515  int tmpnvars = SCIPgetNVarsBounddisjunction(scip, tmpcons);
516  int k;
517 
518  if( checknvars != tmpnvars )
519  continue;
520 
521  for(k=0; k < tmpnvars; k++)
522  if(tmpvars[k] != checkvars[k] || tmpbounds[k] !=checkbounds[k] || tmpboundtypes[k] != checkboundtypes[k])
523  break;
524 
525  if (k == tmpnvars)
526  problem_ -> Jnlst () -> Printf (Ipopt::J_WARNING, J_NLPHEURISTIC, "ZZZ identical bounddisjunction constraints\n");
527  }
528  }
529  }
530 
532 
533 #if 0
534  // writes MILP problem and SCIP settings into a file
535  SCIP_CALL( SCIPwriteOrigProblem(scip, "debug.lp", NULL, FALSE) );
536  SCIP_CALL( SCIPwriteParams(scip, "debug.set", FALSE,TRUE) );
537 #endif
538 
539  // solve the MILP
540  // /|
541  // / +---------+
542  SCIP_RETCODE retcode = SCIPsolve (scip); // < MILP |
543  // \ +---------+
544  // \|
545 #if 0
546  if (SCIPgetStatus (scip) == SCIP_STATUS_INFEASIBLE) {
547 
548  // writes MILP problem and SCIP settings into a file
549  SCIP_CALL( SCIPwriteOrigProblem(scip, "debug.lp", NULL, FALSE) );
550  SCIP_CALL( SCIPwriteParams(scip, "debug.set", FALSE, TRUE));
551 
552  printf ("SCIP found that the problem is infeasible, exiting\n");
553 
554  exit (-1);
555  }
556 #endif
557 
558  if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_WARNING, J_NLPHEURISTIC))
559  SCIP_CALL( SCIPprintStatistics(scip, NULL) );
560 
561  // Greppable line with condensed info
562 
563  problem_ -> Jnlst () -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC, "[FeasPump-SCIP] %5d %5d %7.2f\n", SCIPgetNVars(scip), SCIPgetNConss(scip), SCIPgetSolvingTime(scip));
564 
565  if (retcode != SCIP_OKAY) {
566  problem_ -> Jnlst () -> Printf (Ipopt::J_WARNING, J_NLPHEURISTIC, "Couenne FP: SCIPsolve did not succeed\n");
567  goto TERMINATION;
568  }
569 
570  nscipsols = SCIPgetNSols(scip);
571 
572  // copy the solution
573  if( nscipsols)
574  {
575  SCIP_SOL** scipsols;
576  SCIP_SOL* bestsol;
577  SCIP_Real cutoffbound;
578 
579  int nstoredsols;
580 
581  double objval = 0;
582 
583  /* get incumbent solution */
584  bestsol = SCIPgetBestSol(scip);
585  assert(bestsol != NULL);
586 
587  /* get SCIP solution pool */
588  scipsols = SCIPgetSols(scip);
589  assert(scipsols != NULL);
590 
591  if (!sol)
592  sol = new CouNumber [problem_ -> nVars ()];
593 
594  // get solution values and objective of incumbent
595  SCIP_CALL( SCIPgetSolVals(scip, bestsol, problem_ -> nVars (), vars, sol) );
596  obj = SCIPgetSolOrigObj(scip, bestsol);
597 
598  nstoredsols = 0;
599 
600  // if we use an objective feaspump, the objective function might be negative
601  cutoffbound = obj > 0 ? 2*obj : obj / 2.0;
602 
603  // insert other SCIP solutions into solution pool
604  // do not store too many or too poor solutions
605  for(int i=1; i<nscipsols && nstoredsols < 10 &&
606  (objval = SCIPgetSolOrigObj(scip,scipsols[i])) <= cutoffbound; i++){
607 
608  double* tmpsol = new CouNumber [nvars];
609 
610  // get solution values
611  SCIP_CALL( SCIPgetSolVals(scip, scipsols[i], problem_ -> nVars (), vars, tmpsol) );
612  CouenneFPsolution couennesol = CouenneFPsolution (problem_, tmpsol);
613 
614  delete [] tmpsol;
615 
616  // add solutions to the pool if they are not in the tabu list
617  if ( tabuPool_ . find (couennesol) == tabuPool_ . end ()
618  && pool_ -> Set (). find (couennesol) == pool_ -> Set() . end ()
619  ){
620 
621  pool_ -> Set (). insert (couennesol);
622 
623  ++nstoredsols;
624  }
625  }
626 
627  ++(*nsuciter);
628 
629  // if we succeeded three times in a row, try a cheaper MILP_ solving method next time
630  // TODO: if we want to use time limits, hitting the time limit would be another good reason to switch
631  if( milpMethod_ == 0 && *nsuciter >= 3 && currentmilpmethod < 4 )
632  {
633  ++currentmilpmethod;
634  *nsuciter = 0;
635  }
636  }
637  //try to use a more aggressive, more expensive way to solve the MILP
638  else if( milpMethod_ == 0 && currentmilpmethod > 1 )
639  {
640  --currentmilpmethod;
641  solveagain = true;
642  *nsuciter = 0;
643 
644  // throw away the current solution process
645  SCIP_CALL( SCIPfreeTransform(scip) );
646  }
647  else {
648 
649  obj = COIN_DBL_MAX;
650  }
651 
652  } while (solveagain);
653 
655 
656  TERMINATION:
657 
658  // release variables before freeing them
659  for (int i=0; i<nvars; i++) {
660  SCIP_CALL( SCIPreleaseVar(scip, &vars[i]) );
661  }
662 
663  // free memory
664  SCIPfreeMemoryArray(scip, &vars);
665 
666  // free tabu stuff
667  SCIPfreeMemoryArray(scip, &tabuvars );
668  SCIPfreeMemoryArray(scip, &tabubounds );
669  SCIPfreeMemoryArray(scip, &tabuboundtypes);
670 
671  delete [] lbs;
672  delete [] ubs;
673 
674  SCIP_CALL( SCIPfree(&scip) );
675 
676  BMScheckEmptyMemory();
677 
678  return SCIP_OKAY;
679 }
680 
681 #endif
COUENNE_round
#define COUENNE_round(x)
Definition: CouennePrecisions.hpp:42
Couenne::CouenneProblem::nVars
int nVars() const
Total number of variables.
Definition: CouenneProblem.hpp:380
Couenne::CouenneFeasPump::couenneCG_
CouenneCutGenerator * couenneCG_
CouenneCutGenerator for linearization cuts.
Definition: CouenneFeasPump.hpp:182
Couenne::CouenneFeasPump::milpCuttingPlane_
enum fpCutPlane milpCuttingPlane_
Separate convexification cuts during or after MILP.
Definition: CouenneFeasPump.hpp:238
Couenne::CouenneFeasPump::pool_
CouenneFPpool * pool_
Pool of solutions.
Definition: CouenneFeasPump.hpp:208
j
static char * j
Definition: OSdtoa.cpp:3622
Couenne::CouenneFeasPump::milpMethod_
int milpMethod_
Which SCIP MILP method to use.
Definition: CouenneFeasPump.hpp:250
Couenne::CouenneFPsolution
Class containing a solution with infeasibility evaluation.
Definition: CouenneFPpool.hpp:32
CouenneFeasPump.hpp
cons_rowcuts.h
constraint handler for rowcuts constraints enables separation of convexification cuts during SCIP sol...
end
fint end
Definition: BonFilterSolver.cpp:135
k
void fint fint * k
Definition: BonBqpdSolver.cpp:27
COUENNE_EPS
#define COUENNE_EPS
Definition: CouennePrecisions.hpp:19
Couenne::CouenneFeasPump::tabuPool_
std::set< CouenneFPsolution, compareSol > tabuPool_
Solutions to avoid.
Definition: CouenneFeasPump.hpp:211
Couenne::CouNumber
double CouNumber
main number type in Couenne
Definition: CouenneTypes.hpp:100
Couenne::J_NLPHEURISTIC
const Ipopt::EJournalCategory J_NLPHEURISTIC(Ipopt::J_USER5)
void
static void
Definition: OSdtoa.cpp:1741
CouenneExprVar.hpp
Couenne::CouenneFeasPump::problem_
CouenneProblem * problem_
Couenne representation of the problem.
Definition: CouenneFeasPump.hpp:179
CouenneProblem.hpp
n
void fint * n
Definition: BonFilterSolver.cpp:96
Couenne::CouenneFeasPump::milp_
OsiSolverInterface * milp_
MILP relaxation of the MINLP (used to find integer, non-NLP-feasible solutions)
Definition: CouenneFeasPump.hpp:200
x
void fint fint fint real fint real * x
Definition: BonBqpdSolver.cpp:27
Couenne::isInteger
bool isInteger(CouNumber x)
is this number integer?
Definition: CouenneExpression.hpp:309
Generated by   doxygen 1.8.5