source: stable/0.1/Couenne/src/main/BonCouenneSetup.cpp @ 45

Last change on this file since 45 was 45, checked in by pbelotti, 13 years ago

cut references to disjunctive cuts

File size: 17.9 KB
Line 
1// (C) Copyright International Business Machines Corporation 2007
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Pierre Bonami, International Business Machines Corporation
7//
8// Date : 04/18/2007
9
10#include "BonCouenneSetup.hpp"
11#include "BonInitHeuristic.hpp"
12#include "BonNlpHeuristic.hpp"
13#include "BonCouenneInterface.hpp"
14
15#include "BonGuessHeuristic.hpp"
16#include "CbcCompareActual.hpp"
17
18#include "CouenneObject.hpp"
19#include "CouenneVarObject.hpp"
20#include "CouenneVTObject.hpp"
21#include "CouenneChooseVariable.hpp"
22#include "CouenneChooseStrong.hpp"
23#include "CouenneSolverInterface.hpp"
24#include "CouenneCutGenerator.hpp"
25#include "BonCouenneInfo.hpp"
26#include "BonCbcNode.hpp"
27
28// MILP cuts
29#include "CglGomory.hpp"
30#include "CglProbing.hpp"
31#include "CglKnapsackCover.hpp"
32#include "CglOddHole.hpp"
33#include "CglClique.hpp"
34#include "CglFlowCover.hpp"
35#include "CglMixedIntegerRounding2.hpp"
36#include "CglTwomir.hpp"
37#include "CglPreProcess.hpp"
38#include "CglLandP.hpp"
39#include "CglRedSplit.hpp"
40
41// Ampl includes
42#include "asl.h"
43#include "getstub.h"
44
45
46namespace Bonmin{
47 
48  SmartAsl::~SmartAsl(){
49    //Code from Ipopt::AmplTNLP to free asl
50    if(asl != NULL){
51        if (X0) {
52          delete [] X0;
53          X0 = NULL;
54        }
55        if (havex0) {
56          delete [] havex0;
57          havex0 = NULL;
58        }
59        if (pi0) {
60          delete [] pi0;
61          pi0 = NULL;
62        }
63        if (havepi0) {
64          delete [] havepi0;
65          havepi0 = NULL;
66        }
67        ASL* asl_to_free = (ASL*)asl;
68        ASL_free(&asl_to_free);
69        asl = NULL;
70    }
71    ASL_free(&asl);
72  }
73 
74  CouenneSetup::~CouenneSetup(){
75    //if (CouennePtr_)
76    //delete CouennePtr_;
77  }
78
79  bool CouenneSetup::InitializeCouenne (char **& argv) {
80    /* Get the basic options. */
81    readOptionsFile();
82 
83    /** Change default value for failure behavior so that code doesn't crash
84        when Ipopt does not solve a sub-problem.*/
85
86    options_->SetStringValue ("nlp_failure_behavior", "fathom", "bonmin.");
87
88    gatherParametersValues(options_);
89
90    continuousSolver_ = new CouenneSolverInterface;
91    CouenneInterface * ci = new CouenneInterface;
92    nonlinearSolver_ = ci;
93    /* Read the model in various places. */
94    ci->readAmplNlFile(argv,roptions(),options(),journalist());
95    aslfg_ = new SmartAsl;
96    aslfg_->asl = readASLfg (argv);
97
98    /** Set the output level for the journalist for all Couenne
99     categories.  We probably want to make that a bit more flexible
100     later. */
101    int i;
102
103    options()->GetIntegerValue("boundtightening_print_level", i, "bonmin.");
104    journalist()->GetJournal("console")-> SetPrintLevel(J_BOUNDTIGHTENING, (EJournalLevel) i);
105
106    options()->GetIntegerValue("branching_print_level", i, "bonmin.");
107    journalist()->GetJournal("console")-> SetPrintLevel(J_BRANCHING, (EJournalLevel) i);
108
109    options()->GetIntegerValue("convexifying_print_level", i, "bonmin.");
110    journalist()->GetJournal("console")-> SetPrintLevel(J_CONVEXIFYING, (EJournalLevel) i);
111
112    options()->GetIntegerValue("problem_print_level", i, "bonmin.");
113    journalist()->GetJournal("console")-> SetPrintLevel(J_PROBLEM, (EJournalLevel) i);
114
115    options()->GetIntegerValue("nlpheur_print_level", i, "bonmin.");
116    journalist()->GetJournal("console")-> SetPrintLevel(J_NLPHEURISTIC, (EJournalLevel) i);
117
118    /* Initialize Couenne cut generator.*/
119    //int ivalue, num_points;
120    //options()->GetEnumValue("convexification_type", ivalue,"bonmin.");
121    //options()->GetIntegerValue("convexification_points",num_points,"bonmin.");
122   
123    CouenneCutGenerator * couenneCg = 
124      new CouenneCutGenerator (ci, this, aslfg_->asl, journalist());
125
126    CouenneProblem * couenneProb = couenneCg -> Problem();
127
128    Bonmin::BabInfo * extraStuff = new Bonmin::CouenneInfo(0);
129
130    // as per instructions by John Forrest, to get changed bounds
131    extraStuff -> setExtraCharacteristics (extraStuff -> extraCharacteristics () | 2);
132
133    continuousSolver_ -> setAuxiliaryInfo (extraStuff);
134    delete extraStuff;
135   
136    extraStuff = dynamic_cast <Bonmin::BabInfo *> (continuousSolver_ -> getAuxiliaryInfo ());
137   
138    /* Setup log level*/
139    int lpLogLevel;
140    options()->GetIntegerValue("lp_log_level",lpLogLevel,"bonmin.");
141    continuousSolver_->messageHandler()->setLogLevel(lpLogLevel);
142
143    //////////////////////////////////////////////////////////////
144
145    couenneCg -> Problem () -> setMaxCpuTime (getDoubleParameter (BabSetupBase::MaxTime));
146
147    ci -> extractLinearRelaxation (*continuousSolver_, *couenneCg);
148
149    // In case there are no discrete variables, we have already a
150    // heuristic solution for which create a initialization heuristic
151    if (!(extraStuff -> infeasibleNode ()) &&
152        ci -> isProvenOptimal () && 
153        ci -> haveNlpSolution ()) {
154
155      /// setup initial heuristic (in principle it should only run once...)
156      InitHeuristic* initHeuristic = new InitHeuristic
157        (ci -> getObjValue (), ci -> getColSolution (), *couenneProb);
158      HeuristicMethod h;
159      h.id = "Init Rounding NLP";
160      h.heuristic = initHeuristic;
161      heuristics_.push_back(h);
162    }
163
164    if (extraStuff -> infeasibleNode ()){
165      std::cout<<"Initial linear relaxation constructed by Couenne is infeasible, exiting..."
166               <<std::endl;
167      return false;
168    }
169
170    //continuousSolver_ -> findIntegersAndSOS (false);
171    //addSos (); // only adds embedded SOS objects
172
173    // Add Couenne SOS ///////////////////////////////////////////////////////////////
174
175    std::string s;
176    int nSOS = 0;
177
178    // allocate sufficient space for both nonlinear variables and SOS's
179    OsiObject ** objects;
180
181    options () -> GetStringValue ("enable_sos", s, "couenne.");
182    if (s == "yes") {
183
184      objects = new OsiObject* [couenneProb -> nVars ()];
185
186      nSOS = couenneProb -> findSOS (nonlinearSolver (), objects);
187      //printf ("==================== found %d SOS\n", nSOS);
188      //nonlinearSolver () -> addObjects (nSOS, objects);
189      continuousSolver () -> addObjects (nSOS, objects);
190
191      for (int i=0; i<nSOS; i++)
192        delete objects [i];
193      delete [] objects;
194    }
195
196    //model -> assignSolver (continuousSolver_, true);
197    //continuousSolver_ = model -> solver();
198
199    // Add Couenne objects for branching /////////////////////////////////////////////
200
201    options () -> GetStringValue ("display_stats", s, "couenne.");
202    displayStats_ = (s == "yes");
203
204    options () -> GetStringValue ("branching_object", s, "couenne.");
205
206    enum CouenneObject::branch_obj objType = CouenneObject::VAR_OBJ;
207
208    if      (s == "vt_obj")   objType = CouenneObject::VT_OBJ;
209    else if (s == "var_obj")  objType = CouenneObject::VAR_OBJ;
210    else if (s == "expr_obj") objType = CouenneObject::EXPR_OBJ;
211    else {
212      printf ("CouenneSetup: Unknown branching object type\n");
213      exit (-1);
214    }
215
216    int 
217      nobj  = 0, // if no SOS then objects is empty
218      nVars = couenneProb -> nVars ();
219
220    objects = new OsiObject* [couenneProb -> nVars ()];
221
222    int contObjPriority = 2000; // default object priority -- it is 1000 for integers and 10 for SOS
223
224    options () -> GetIntegerValue ("cont_var_priority", contObjPriority, "bonmin.");
225
226    for (int i = 0; i < nVars; i++) { // for each variable
227
228      exprVar *var = couenneProb -> Var (i);
229
230      // we only want enabled variables
231      if (var -> Multiplicity () <= 0) 
232        continue;
233
234      switch (objType) {
235
236      case CouenneObject::EXPR_OBJ:
237
238        // if this variable is associated with a nonlinear function
239        if (var -> isInteger () || 
240            (var -> Type  () == AUX) && 
241            (var -> Image () -> Linearity () > LINEAR)) {
242
243          objects [nobj] = new CouenneObject (couenneProb, var, this, journalist ());
244          objects [nobj++] -> setPriority (contObjPriority);
245          //objects [nobj++] -> setPriority (contObjPriority + var -> rank ());
246        }
247
248        break;
249
250      case CouenneObject::VAR_OBJ:
251
252        // branching objects on variables
253        if // comment three lines below for linear variables too
254          (var -> isInteger () || 
255           (couenneProb -> Dependence () [var -> Index ()] . size () > 0)) {  // has indep
256           //|| ((var -> Type () == AUX) &&                                  // or, aux
257           //    (var -> Image () -> Linearity () > LINEAR))) {              // of nonlinear
258
259          objects [nobj] = new CouenneVarObject (couenneProb, var, this, journalist ());
260          objects [nobj++] -> setPriority (contObjPriority);
261          //objects [nobj++] -> setPriority (contObjPriority + var -> rank ());
262        }
263
264        break;
265
266      default:
267      case CouenneObject::VT_OBJ:
268
269        // branching objects on variables
270        if // comment three lines below for linear variables too
271          (var -> isInteger () || 
272           (couenneProb -> Dependence () [var -> Index ()] . size () > 0)) { // has indep
273          //|| ((var -> Type () == AUX) &&                      // or, aux
274          //(var -> Image () -> Linearity () > LINEAR))) { // of nonlinear
275
276          objects [nobj] = new CouenneVTObject (couenneProb, var, this, journalist ());
277          objects [nobj++] -> setPriority (contObjPriority);
278          //objects [nobj++] -> setPriority (contObjPriority + var -> rank ());
279        }
280
281        break;
282      }
283    }
284
285    // Add objects /////////////////////////////////
286
287    continuousSolver_ -> addObjects (nobj, objects);
288
289    for (int i = 0 ; i < nobj ; i++)
290      delete objects [i];
291
292    delete [] objects;
293
294    // Setup Convexifier generators ////////////////////////////////////////////////
295
296    int freq;
297
298    options()->GetIntegerValue("convexification_cuts",freq,"couenne.");
299
300    if (freq != 0) {
301
302      CuttingMethod cg;
303      cg.frequency = freq;
304      cg.cgl = couenneCg;
305      cg.id = "Couenne convexifier cuts";
306      cutGenerators().push_back(cg);
307
308      // set cut gen pointer
309      dynamic_cast <CouenneSolverInterface *> 
310        (continuousSolver_) -> setCutGenPtr (couenneCg);
311    }
312
313    // add other cut generators -- test for integer variables first
314    if (couenneCg -> Problem () -> nIntVars () > 0)
315      addMilpCutGenerators ();
316
317    CouennePtr_ = couenneCg;
318
319    // Setup heuristic to solve nlp problems. /////////////////////////////////
320
321    int doNlpHeurisitic = 0;
322    options()->GetEnumValue("local_optimization_heuristic", doNlpHeurisitic, "couenne.");
323    if(doNlpHeurisitic)
324    {
325      int numSolve;
326      options()->GetIntegerValue("log_num_local_optimization_per_level",numSolve,"couenne.");
327      NlpSolveHeuristic * nlpHeuristic = new NlpSolveHeuristic;
328      nlpHeuristic->setNlp(*ci,false);
329      nlpHeuristic->setCouenneProblem(couenneProb);
330      //nlpHeuristic->setMaxNlpInf(1e-4);
331      nlpHeuristic->setMaxNlpInf(maxNlpInf_0);
332      nlpHeuristic->setNumberSolvePerLevel(numSolve);
333      HeuristicMethod h;
334      h.id = "Couenne Rounding NLP";
335      h.heuristic = nlpHeuristic;
336      heuristics_.push_back(h);
337    }
338
339#if 0
340    {
341      CbcCompareEstimate compare;
342      model -> setNodeComparison(compare);
343      GuessHeuristic * guessHeu = new GuessHeuristic (*model);
344      HeuristicMethod h;
345      h.id = "Bonmin Guessing";
346      h.heuristic = guessHeu;
347      heuristics_.push_back (h);
348      //model_.addHeuristic(guessHeu);
349      //delete guessHeu;
350    }
351#endif
352
353    // Add Branching rules ///////////////////////////////////////////////////////
354
355    int varSelection;
356    if (!options_->GetEnumValue("variable_selection",varSelection,"bonmin.")) {
357      // change the default for Couenne
358      varSelection = OSI_SIMPLE;
359    }
360
361    switch (varSelection) {
362
363    case OSI_STRONG: { // strong branching
364      CouenneChooseStrong * chooseVariable = new CouenneChooseStrong
365        (*this, couenneProb, journalist ());
366      chooseVariable->setTrustStrongForSolution(false);
367      chooseVariable->setTrustStrongForBound(false);
368      chooseVariable->setOnlyPseudoWhenTrusted(true);
369      branchingMethod_ = chooseVariable;
370      break;
371    }
372
373    case OSI_SIMPLE: // default choice
374      branchingMethod_ = new CouenneChooseVariable
375        (continuousSolver_, couenneProb, journalist ());
376      break;
377
378    default:
379      std::cerr << "Unknown variable_selection for Couenne\n" << std::endl;
380      throw;
381      break;
382    }
383
384    int ival;
385    if (!options_->GetEnumValue("node_comparison",ival,"bonmin.")) {
386      // change default for Couenne
387      nodeComparisonMethod_ = bestBound;
388    }
389    else {
390      nodeComparisonMethod_ = NodeComparison(ival);
391    }
392
393    if(intParam_[NumCutPasses] < 2)
394    intParam_[NumCutPasses] = 2;
395
396    // Tell Cbc not to check again if a solution returned from
397    // heuristic is indeed feasible
398    intParam_[BabSetupBase::SpecialOption] = 16 | 4;
399
400    return true;
401  }
402 
403  void CouenneSetup::registerOptions(){
404    registerAllOptions(roptions());
405  }
406
407
408  void
409  CouenneSetup::registerAllOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions){
410    BabSetupBase::registerAllOptions(roptions);
411    BonCbcFullNodeInfo::registerOptions(roptions);
412    CouenneCutGenerator::registerOptions (roptions);
413
414    roptions -> AddNumberOption
415      ("couenne_check",
416       "known value of a global optimum",
417       COIN_DBL_MAX,
418       "Default value is +infinity.");
419
420    roptions -> AddStringOption2 (
421      "display_stats",
422      "display statistics at the end of the run",
423      "no",
424      "yes", "",
425      "no", "");
426
427    roptions->AddBoundedIntegerOption(
428      "branching_print_level",
429      "Output level for braching code in Couenne",
430      -2, J_LAST_LEVEL-1, J_NONE,
431      "");
432    roptions->AddBoundedIntegerOption(
433      "boundtightening_print_level",
434      "Output level for bound tightening code in Couenne",
435      -2, J_LAST_LEVEL-1, J_NONE,
436      "");
437    roptions->AddBoundedIntegerOption(
438      "convexifying_print_level",
439      "Output level for convexifying code in Couenne",
440      -2, J_LAST_LEVEL-1, J_NONE,
441      "");
442    roptions->AddBoundedIntegerOption(
443      "problem_print_level",
444      "Output level for problem manipulation code in Couenne",
445      -2, J_LAST_LEVEL-1, J_WARNING,
446      "");
447    roptions->AddBoundedIntegerOption(
448      "nlpheur_print_level",
449      "Output level for NLP heuristic in Couenne",
450      -2, J_LAST_LEVEL-1, J_WARNING,
451      "");
452
453
454    // copied from BonminSetup::registerMilpCutGenerators(), in
455    // BonBonminSetup.cpp
456
457    struct cutOption_ {
458
459      const char *cgname;
460      int         defaultFreq;
461
462    } cutOption [] = {
463      {(const char *) "Gomory_cuts",             0},
464      {(const char *) "probing_cuts",            0},
465      {(const char *) "cover_cuts",              0},
466      {(const char *) "mir_cuts",                0},
467      {(const char *) "2mir_cuts",               0},
468      {(const char *) "flow_covers_cuts",        0},
469      {(const char *) "lift_and_project_cuts",   0},
470      {(const char *) "reduce_split_cuts",       0},
471      {(const char *) "clique_cuts",             0},
472      {NULL, 0}};
473
474    for (int i=0; cutOption [i].cgname; i++) {
475
476      char descr [150];
477
478      sprintf (descr, "Frequency k (in terms of nodes) for generating %s cuts in branch-and-cut.",
479               cutOption [i].cgname);
480
481      roptions -> AddLowerBoundedIntegerOption
482        (cutOption [i].cgname,
483         descr,
484         -100, cutOption [i].defaultFreq,
485         "If k > 0, cuts are generated every k nodes, "
486         "if -99 < k < 0 cuts are generated every -k nodes but "
487         "Cbc may decide to stop generating cuts, if not enough are generated at the root node, "
488         "if k=-99 generate cuts only at the root node, if k=0 or 100 do not generate cuts.");
489
490      roptions->setOptionExtraInfo (cutOption [i].cgname, 5);
491    }
492  }
493
494
495
496  /** Add milp cut generators according to options.*/
497  void CouenneSetup::addMilpCutGenerators () {
498
499    enum extraInfo_ {CUTINFO_NONE, CUTINFO_MIG, CUTINFO_PROBING, CUTINFO_CLIQUE};
500
501    struct cutInfo {
502
503      const char      *optname;
504      CglCutGenerator *cglptr;
505      const char      *cglId;
506      enum extraInfo_  extraInfo;
507
508    } cutList [] = {
509      {(const char*)"Gomory_cuts",new CglGomory,      (const char*)"Mixed Integer Gomory",CUTINFO_MIG},
510      {(const char*)"probing_cuts",new CglProbing,        (const char*) "Probing", CUTINFO_PROBING},
511      {(const char*)"mir_cuts",new CglMixedIntegerRounding2, (const char*) "Mixed Integer Rounding", 
512                                                                                        CUTINFO_NONE},
513      {(const char*)"2mir_cuts",    new CglTwomir,         (const char*) "2-MIR",    CUTINFO_NONE},
514      {(const char*)"cover_cuts",   new CglKnapsackCover,  (const char*) "Cover",    CUTINFO_NONE},
515      {(const char*)"clique_cuts",  new CglClique,         (const char*) "Clique",   CUTINFO_CLIQUE},
516      {(const char*)"lift_and_project_cuts",new CglLandP,(const char*)"Lift and Project",CUTINFO_NONE},
517      {(const char*)"reduce_split_cuts",new CglRedSplit,(const char*) "Reduce and Split",CUTINFO_NONE},
518      {(const char*)"flow_covers_cuts",new CglFlowCover,(const char*) "Flow cover cuts", CUTINFO_NONE},
519      {NULL, NULL, NULL, CUTINFO_NONE}};
520
521    int freq;
522
523    for (int i=0; cutList [i]. optname; i++) {
524
525      options_ -> GetIntegerValue (std::string (cutList [i]. optname), freq, "bonmin.");
526
527      if (!freq) {
528        delete cutList [i].cglptr;
529        continue;
530      }
531
532      CuttingMethod cg;
533      cg.frequency = freq;
534      cg.cgl       = cutList [i].cglptr;
535      cg.id        = std::string (cutList [i]. cglId);
536      cutGenerators_.push_back (cg);
537
538      // options for particular cases
539      switch (cutList [i].extraInfo) {
540
541      case CUTINFO_MIG: {
542        CglGomory *gc = dynamic_cast <CglGomory *> (cutList [i].cglptr);
543
544        if (!gc) break;
545
546        gc -> setLimitAtRoot(512);
547        gc -> setLimit(50);
548      }
549        break;
550
551      case CUTINFO_PROBING: {
552        CglProbing *pc = dynamic_cast <CglProbing *> (cutList [i].cglptr);
553
554        if (!pc) break;
555
556        pc->setUsingObjective(1);
557        pc->setMaxPass(3);
558        pc->setMaxPassRoot(3);
559        // Number of unsatisfied variables to look at
560        pc->setMaxProbe(10);
561        pc->setMaxProbeRoot(50);
562        // How far to follow the consequences
563        pc->setMaxLook(10);
564        pc->setMaxLookRoot(50);
565        pc->setMaxLookRoot(10);
566        // Only look at rows with fewer than this number of elements
567        pc->setMaxElements(200);
568        pc->setRowCuts(3);
569      }
570        break;
571
572      case CUTINFO_CLIQUE: {
573        CglClique *clique = dynamic_cast <CglClique *> (cutList [i].cglptr);
574
575        if (!clique) break;
576
577        clique -> setStarCliqueReport(false);
578        clique -> setRowCliqueReport(false);
579        clique -> setMinViolation(0.1);
580      }
581        break;
582
583        //case CUTINFO_NONE:
584      default:
585        break;
586      }
587    }
588  }
589}
Note: See TracBrowser for help on using the repository browser.