source: trunk/Bonmin/src/Algorithms/OaGenerators/BonOACutGenerator2.cpp @ 1536

Last change on this file since 1536 was 1536, checked in by pbonami, 10 years ago

Port Claudia's code to trunk

File size: 9.5 KB
Line 
1// (C) Copyright Carnegie Mellon University 2005, 2006
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// P. Bonami, Carnegie Mellon University
7//
8// Date : 05/26/2005
9
10#include "BonOACutGenerator2.hpp"
11#include "BonminConfig.h"
12
13#include "OsiClpSolverInterface.hpp"
14
15#include "CbcModel.hpp"
16#include "BonCbcLpStrategy.hpp"
17#ifdef COIN_HAS_CPX
18#include "OsiCpxSolverInterface.hpp"
19#endif
20#include "OsiAuxInfo.hpp"
21#include "BonSolverHelp.hpp"
22
23#include <climits>
24
25namespace Bonmin
26{
27   static const char * txt_id = "OA decomposition";
28
29
30/// Constructor with basic setup
31  OACutGenerator2::OACutGenerator2(BabSetupBase & b):
32      OaDecompositionBase(b, true, false)
33  {
34    int ivalue;
35    std::string bonmin="bonmin.";
36    std::string prefix = (b.prefix() == bonmin) ? "" : b.prefix();
37    prefix += "oa_decomposition.";
38    b.options()->GetEnumValue("milp_solver",ivalue,prefix);
39    if (ivalue <= 0) {//uses cbc
40      CbcStrategyDefault strategy;
41      setStrategy(strategy);
42    }
43    else if (ivalue == 1) {
44      CbcStrategyChooseCuts strategy(b, prefix);
45      setStrategy(strategy);
46    }
47    else if (ivalue == 2) {
48#ifdef COIN_HAS_CPX
49      OsiCpxSolverInterface * cpxSolver = new OsiCpxSolverInterface;
50      b.nonlinearSolver()->extractLinearRelaxation(*cpxSolver);
51      assignLpInterface(cpxSolver);
52#else
53      std::cerr << "You have set an option to use CPLEX as the milp\n"
54      << "subsolver in oa decomposition. However, apparently\n"
55      << "CPLEX is not configured to be used in bonmin.\n"
56      << "See the manual for configuring CPLEX\n";
57      throw -1;
58#endif
59    }
60
61    double oaTime;
62    b.options()->GetNumericValue("time_limit",oaTime,prefix);
63    parameter().maxLocalSearchTime_ =
64    std::min(b.getDoubleParameter(BabSetupBase::MaxTime), oaTime);
65    parameter().maxLocalSearch_ = INT_MAX;
66    b.options()->GetIntegerValue("solution_limit", parameter().maxSols_,prefix);
67  }
68  OACutGenerator2::~OACutGenerator2()
69  {}
70
71  /// virutal method to decide if local search is performed
72  bool
73  OACutGenerator2::doLocalSearch(BabInfo * babInfo) const
74  {
75    return (nLocalSearch_<parameters_.maxLocalSearch_ &&
76            CoinCpuTime() - timeBegin_ < parameters_.maxLocalSearchTime_);
77  }
78  /// virtual method which performs the OA algorithm by modifying lp and nlp.
79  double
80  OACutGenerator2::performOa(OsiCuts &cs,
81      solverManip &lpManip,
82      SubMipSolver * &subMip,
83      BabInfo * babInfo,
84      double & cutoff, const CglTreeInfo & info) const
85  {
86
87    double lastPeriodicLog = CoinCpuTime();
88
89    //const int numcols = nlp_->getNumCols();
90
91
92    bool isInteger = false;
93
94    OsiSolverInterface * lp = lpManip.si();
95    OsiBranchingInformation branch_info(lp, false);
96    bool milpOptimal = 1;
97
98
99    double milpBound = -COIN_DBL_MAX;
100    bool milpFeasible = 1;
101    bool feasible = 1;
102
103    if (subMip)//Perform a local search
104    {
105      subMip->find_good_sol(cutoff, parameters_.subMilpLogLevel_,
106          (parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime()));
107      milpBound = subMip->lowBound();
108      milpOptimal = subMip->optimal();
109
110      feasible = milpBound < cutoff;
111      milpFeasible = feasible;
112      isInteger = (subMip->getLastSolution() != NULL);
113      nLocalSearch_++;
114
115      if (milpOptimal)
116        handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
117      else
118      {
119        handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
120      }
121    }
122    int numberPasses = 0;
123
124#ifdef OA_DEBUG
125    bool foundSolution = 0;
126#endif
127    double * nlpSol = NULL;
128
129    while (isInteger && feasible ) {
130      numberPasses++;
131      //after a prescribed elapsed time give some information to user
132      double time = CoinCpuTime();
133      if (time - lastPeriodicLog > parameters_.logFrequency_) {
134        double lb = (subMip == NULL) ?lp->getObjValue() : subMip->getLowerBound();
135        handler_->message(PERIODIC_MSG,messages_)
136        <<time - timeBegin_<<cutoff
137        <<lb
138        <<CoinMessageEol;
139        lastPeriodicLog = CoinCpuTime();
140      }
141
142
143      //setup the nlp
144      int numberCutsBefore = cs.sizeRowCuts();
145
146      //Fix the variable which have to be fixed, after having saved the bounds
147      const double * colsol = subMip == NULL ? lp->getColSolution():
148          subMip->getLastSolution();
149      branch_info.solution_ = colsol;
150
151      fixIntegers(*nlp_,branch_info, parameters_.cbcIntegerTolerance_,objects_, nObjects_);
152
153      nlp_->resolve(txt_id);
154      if (post_nlp_solve(babInfo, cutoff)) {
155        //nlp solved and feasible
156        // Update the cutoff
157        cutoff = nlp_->getObjValue() *(1 - parameters_.cbcCutoffIncrement_);
158        // Update the lp solver cutoff
159        lp->setDblParam(OsiDualObjectiveLimit, cutoff);
160        numSols_++;
161      }
162
163      nlpSol = const_cast<double *>(nlp_->getColSolution());
164
165      // Get the cuts outer approximation at the current point
166      const double * toCut = (parameter().addOnlyViolated_)?
167          colsol:NULL;
168      nlp_->getOuterApproximation(cs, nlpSol, 1, toCut,
169                                  parameter().global_);
170
171      int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
172      assert(numberCuts);
173      installCuts(*lp, cs, numberCuts);
174
175      lp->resolve();
176
177      double objvalue = lp->getObjValue();
178      //milpBound = max(milpBound, lp->getObjValue());
179      feasible = (lp->isProvenOptimal() &&
180          !lp->isDualObjectiveLimitReached() && (objvalue<cutoff)) ;
181      //if value of integers are unchanged then we have to get out
182      bool changed = !feasible;//if lp is infeasible we don't have to check anything
183      branch_info.solution_ = lp->getColSolution();
184      if (!changed)
185        changed = isDifferentOnIntegers(*nlp_, objects_, nObjects_,
186                                        0.1,
187                                        nlp_->getColSolution(), lp->getColSolution());
188      if (changed) {
189
190        isInteger = integerFeasible(*lp, branch_info, parameters_.cbcIntegerTolerance_,
191                                     objects_, nObjects_);
192      }
193      else {
194        isInteger = 0;
195        //        if(!fixed)//fathom on bounds
196        milpBound = 1e200;
197      }
198#ifdef OA_DEBUG
199      printf("Obj value after cuts %g %d rows\n",lp->getObjValue(),
200          numberCuts) ;
201#endif
202      //check time
203      if (CoinCpuTime() - timeBegin_ > parameters_.maxLocalSearchTime_)
204        break;
205      //do we perform a new local search ?
206      if (feasible && !isInteger &&
207          nLocalSearch_ < parameters_.maxLocalSearch_ &&
208          numSols_ < parameters_.maxSols_) {
209
210        /** do we have a subMip? if not create a new one. */
211        if (subMip == NULL) subMip = new SubMipSolver(lp, parameters_.strategy());
212
213        nLocalSearch_++;
214
215        subMip->find_good_sol(cutoff, parameters_.subMilpLogLevel_,
216            parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime()
217            );
218
219        milpBound = subMip->lowBound();
220
221        if (subMip->optimal())
222          handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
223        else
224          handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
225
226
227        colsol = const_cast<double *> (subMip->getLastSolution());
228        isInteger = (colsol != 0);
229
230        feasible =  (milpBound < cutoff);
231
232        if (feasible && isInteger) {
233          bool changed = false;
234          changed = isDifferentOnIntegers(*nlp_, objects_, nObjects_,
235                                          0.1,
236                                          nlp_->getColSolution(), colsol);
237          //solution problem is solved
238          if (!changed) {
239            feasible = 0;
240            milpBound = 1e50;
241          }
242          milpFeasible = feasible;
243        }
244        if (subMip->optimal())
245          milpOptimal = 1;
246        else {
247          milpOptimal = 0;
248        }
249
250        if (milpBound < cutoff)
251          handler_->message(UPDATE_LB, messages_)<<milpBound<<CoinCpuTime() - timeBegin_<<CoinMessageEol;
252        else {
253          milpBound = 1e50;
254          feasible = 0;
255          handler_->message(OASUCCESS, messages_)<<CoinCpuTime() - timeBegin_ <<CoinMessageEol;
256        }
257      }/** endif localSearch*/
258      else if (subMip!=NULL) {
259        delete subMip;
260        subMip = NULL;
261      }
262    }
263
264#ifdef OA_DEBUG
265    debug_.printEndOfProcedureDebugMessage(cs, foundSolution, cutoff, milpBound, isInteger, feasible, std::cout);
266#endif
267    return milpBound;
268  }
269
270  /** Register OA options.*/
271  void
272  OACutGenerator2::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
273  {
274    roptions->SetRegisteringCategory("Options for OA decomposition", RegisteredOptions::BonminCategory);
275    roptions->AddStringOption2("oa_decomposition", "If yes do initial OA decomposition",
276                               "no",
277                               "no","",
278                               "yes","",
279                               "");
280    roptions->AddBoundedIntegerOption("oa_log_level",
281        "specify OA iterations log level.",
282        0,2,1,
283        "Set the level of output of OA decomposition solver : "
284        "0 - none, 1 - normal, 2 - verbose"
285                                     );
286    roptions->setOptionExtraInfo("oa_decomposition", 7);
287
288    roptions->AddLowerBoundedNumberOption("oa_log_frequency",
289        "display an update on lower and upper bounds in OA every n seconds",
290        0.,1.,100.,
291        "");
292    roptions->setOptionExtraInfo("oa_log_frequency", 15);
293  }
294}/* End namespace Bonmin. */
295
Note: See TracBrowser for help on using the repository browser.