source: trunk/Bonmin/experimental/Bcp/BM_lp.cpp @ 592

Last change on this file since 592 was 592, checked in by andreasw, 12 years ago

enabling compilation on MSVC++

File size: 13.2 KB
Line 
1// (C) Copyright International Business Machines Corporation 2006, 2007
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Laszlo Ladanyi, International Business Machines Corporation
7// Pierre Bonami, Carnegie Mellon University
8
9#include "OsiClpSolverInterface.hpp"
10#include "BM.hpp"
11#include "BCP_message_mpi.hpp"
12#include "BCP_lp_node.hpp"
13#include "BCP_lp.hpp"
14
15#include "BonOACutGenerator2.hpp"
16#include "BonEcpCuts.hpp"
17#include "BonOaNlpOptim.hpp"
18
19// The following is included for "min"
20#include "CoinFinite.hpp"
21
22//#############################################################################
23
24BM_lp::BM_lp() :
25    BCP_lp_user(),
26    babSolver_(3),
27    bonmin_(),
28    ws_(NULL),
29    chooseVar_(NULL),
30    freeToBusyRatio_(0.0),
31    in_strong(0)
32{
33      setOsiBabSolver(&babSolver_);
34}
35
36/****************************************************************************/
37
38BM_lp::~BM_lp()
39{
40    delete chooseVar_;
41    delete ws_;
42    delete[] primal_solution_;
43    /* FIXME: CHEATING */
44    for (std::map<int, CoinWarmStart*>::iterator it = warmStart.begin();
45         it != warmStart.end(); ++it) {
46        delete it->second;
47    }
48}
49
50/****************************************************************************/
51
52OsiSolverInterface *
53BM_lp::initialize_solver_interface()
54{
55    OsiClpSolverInterface * clp = new OsiClpSolverInterface;
56    OsiBabSolver babSolver(3);
57    babSolver.setSolver(clp);
58    clp->setAuxiliaryInfo(&babSolver);
59    clp->messageHandler()->setLogLevel(0);
60    setOsiBabSolver(dynamic_cast<OsiBabSolver *>(clp->getAuxiliaryInfo()));
61
62    return clp;
63}
64
65/****************************************************************************/
66
67void
68BM_lp::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
69                                       const BCP_vec<BCP_cut*>& cuts,
70                                       const BCP_vec<BCP_obj_status>& vs,
71                                       const BCP_vec<BCP_obj_status>& cs,
72                                       BCP_vec<int>& var_changed_pos,
73                                       BCP_vec<double>& var_new_bd,
74                                       BCP_vec<int>& cut_changed_pos,
75                                       BCP_vec<double>& cut_new_bd)
76{
77    BM_node* data = dynamic_cast<BM_node*>(get_user_data());
78    if (data) {
79        numNlpFailed_ = data->numNlpFailed_;
80    }
81
82    int i;
83    // First copy the bounds into nlp. That way all the branching decisions
84    // will be transferred over.
85    OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
86    Bonmin::OsiTMINLPInterface& nlp = *bonmin_.nonlinearSolver();
87    nlp.setColLower(osi->getColLower());
88    nlp.setColUpper(osi->getColUpper());
89
90    // Carry the changes over to the object lists in nlp
91    const int numObj = nlp.numberObjects();
92    OsiObject** nlpObj = nlp.objects();
93    for (i = 0; i < numObj; ++i) {
94        OsiSimpleInteger* io = dynamic_cast<OsiSimpleInteger*>(nlpObj[i]);
95        if (io) {
96            io->resetBounds(&nlp);
97        } else {
98            // The rest is OsiSOS where we don't need to do anything
99            break;
100        }
101    }
102
103    // copy over the OsiObjects from the nlp solver if the lp solver is to be
104    // used at all (i.e., not pure B&B)
105    if (bonmin_.getAlgorithm() > 0) {
106        osi->addObjects(nlp.numberObjects(), nlp.objects());
107    }
108
109    in_strong = 0;
110}
111
112/************************************************************************/
113void
114BM_lp::modify_lp_parameters(OsiSolverInterface* lp, bool in_strong_branching)
115
116    // Called each time the node LP is solved
117
118{
119    if (in_strong_branching) {
120        in_strong = 1;
121        //     lp->setIntParam(OsiMaxNumIterationHotStart, 50);
122    }
123}
124
125/****************************************************************************/
126
127BCP_solution*
128BM_lp::test_feasibility(const BCP_lp_result& lp_result,
129                        const BCP_vec<BCP_var*>& vars,
130                        const BCP_vec<BCP_cut*>& cuts)
131{
132    if (bonmin_.getAlgorithm() == 0) {
133        // if pure B&B
134        return test_feasibility_BB(vars);
135    } else {
136        return test_feasibility_hybrid(lp_result, vars, cuts);
137    }
138    return NULL; // fake return to quiet gcc
139}
140
141/****************************************************************************/
142
143BCP_solution*
144BM_lp::test_feasibility_BB(const BCP_vec<BCP_var*>& vars)
145{
146  Bonmin::OsiTMINLPInterface& nlp = *bonmin_.nonlinearSolver();
147    /* PURE_BB TODO:
148       Do whatever must be done to:
149       1) compute a lower bound using NLP and save it in "lower_bound_"
150       2) save the solution of the NLP in "primal_solution_"
151       3) if the optimum (well, local opt anyway...) happens to be
152       feasible (or along the NLP optimization we have blundered into
153       a feas sol) then create "sol"
154    */
155    BM_solution* sol = NULL;
156
157    char prefix[100];
158#ifdef COIN_HAS_MPI
159    const BCP_proc_id* id = getLpProblemPointer()->get_process_id();
160    const BCP_mpi_id* mid = dynamic_cast<const BCP_mpi_id*>(id);
161    sprintf(prefix, "%i", mid->pid());
162#else
163    prefix[0] = 0;
164#endif
165    try {
166    switch (par.entry(BM_par::WarmStartStrategy)) {
167    case WarmStartNone:
168        nlp.initialSolve();
169        break;
170    case WarmStartFromRoot:
171        nlp.setWarmStart(ws_);
172        nlp.resolve();
173        break;
174    case WarmStartFromParent:
175        /* FIXME: CHEAT! this works only in serial mode! */
176        {
177            const int ind = current_index();
178            const int parentind =
179                getLpProblemPointer()->parent->index;
180            std::map<int, CoinWarmStart*>::iterator it =
181                warmStart.find(parentind);
182            nlp.setWarmStart(it->second);
183            nlp.resolve();
184            warmStart[ind] = nlp.getWarmStart();
185            bool sibling_seen =  ((ind & 1) == 0) ?
186                warmStart.find(ind-1) != warmStart.end() :
187                warmStart.find(ind+1) != warmStart.end() ;
188            if (sibling_seen) {
189                delete it->second;
190                warmStart.erase(it);
191            }
192        }
193        break;
194    }
195    }
196    catch(Bonmin::TNLPSolver::UnsolvedError &E) {
197      E.writeDiffFiles(prefix);
198      E.printError(std::cerr);
199      //There has been a failure to solve a problem with Ipopt.  And
200      //we will output file with information on what has been changed
201      //in the problem to make it fail.
202      //Now depending on what algorithm has been called (B-BB or
203      //other) the failed problem may be at different place.
204      //    const OsiSolverInterface &si1 = (algo > 0) ? nlpSolver : *model.solver();
205    }
206    catch(Bonmin::OsiTMINLPInterface::SimpleError &E) {
207      std::cerr<<E.className()<<"::"<<E.methodName()
208               <<std::endl
209               <<E.message()<<std::endl;
210    }
211    catch(CoinError &E) {
212      std::cerr<<E.className()<<"::"<<E.methodName()
213               <<std::endl
214               <<E.message()<<std::endl;
215    }
216    catch (Ipopt::OPTION_INVALID &E) {
217      std::cerr<<"Ipopt exception : "<<E.Message()<<std::endl;
218    }
219    catch(...) {
220      std::cerr<<" unrecognized exception"<<std::endl;
221      throw;
222    }
223
224    const int numCols = nlp.getNumCols();
225    const double* colsol = nlp.getColSolution();
226    if (nlp.isProvenOptimal()) {
227        int i;
228        const double* clb = nlp.getColLower();
229        const double* cub = nlp.getColUpper();
230        // Make sure we are within bounds (get rid of rounding problems)
231        for (i = 0; i < numCols; ++i) {
232            primal_solution_[i] =
233                CoinMin(CoinMax(clb[i], colsol[i]), cub[i]);
234        }
235       
236        lower_bound_ = nlp.getObjValue();
237        numNlpFailed_ = 0;
238
239        for (i = 0; i < numCols; ++i) {
240            if (vars[i]->var_type() == BCP_ContinuousVar)
241                continue;
242            const double psol = primal_solution_[i];
243            const double frac = psol - floor(psol);
244            const double dist = min(frac, 1.0-frac);
245            if (dist >= integerTolerance_)
246                break;
247        }
248        if (i == numCols) {
249            /* yipee! a feasible solution! */
250            sol = new BM_solution;
251            //Just copy the solution stored in solver to sol
252            double ptol;
253            nlp.getDblParam(OsiPrimalTolerance, ptol);
254            for (i = 0 ; i < numCols ; i++) {
255                if (fabs(primal_solution_[i]) > ptol)
256                    sol->add_entry(i, primal_solution_[i]); 
257            }
258            sol->setObjective(lower_bound_);
259        }
260        if (lower_bound_ > upper_bound()-get_param(BCP_lp_par::Granularity) &&
261            par.entry(BM_par::PrintBranchingInfo)) {
262            printf("\
263BM_lp: At node %i : will fathom because of high lower bound\n",
264                   current_index());
265        }
266    }
267    else if (nlp.isProvenPrimalInfeasible()) {
268        // prune it!
269        // FIXME: if nonconvex, restart from a different place...
270        lower_bound_ = 1e200;
271        numNlpFailed_ = 0;
272        if (par.entry(BM_par::PrintBranchingInfo)) {
273            printf("\
274BM_lp: At node %i : will fathom because of infeasibility\n",
275                   current_index());
276        }
277    }
278    else if (nlp.isAbandoned()) {
279        if (nlp.isIterationLimitReached()) {
280            printf("\
281BM_lp: At node %i : WARNING: nlp reached iter limit. Will force branching\n",
282                   current_index());
283        } else {
284            printf("\
285BM_lp: At node %i : WARNING: nlp is abandoned. Will force branching\n",
286                   current_index());
287        }
288        // nlp failed
289        nlp.forceBranchable();
290        lower_bound_ = nlp.getObjValue();
291        CoinDisjointCopyN(colsol, numCols, primal_solution_);
292        numNlpFailed_ += 1;
293        // will branch, but save in the user data how many times we have
294        // failed, and if we fail too many times in a row then just abandon
295        // the node.
296    }
297    else {
298        // complain loudly about a bug
299        throw BCP_fatal_error("Impossible outcome by nlp.initialSolve()\n");
300    }
301    return sol;
302}
303
304/****************************************************************************/
305
306BCP_solution*
307BM_lp::test_feasibility_hybrid(const BCP_lp_result& lp_result,
308                               const BCP_vec<BCP_var*>& vars,
309                               const BCP_vec<BCP_cut*>& cuts)
310{
311    /* First test that the integrality requirements are met. */
312    BCP_solution_generic* gsol = test_full(lp_result, vars, integerTolerance_);
313    if (! gsol) {}
314       
315    /* TODO:
316       don't test feasibility in every node
317    */
318   
319
320    OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
321   
322    // Don't test feasibility if the node LP was infeasible
323    if (osi->isProvenPrimalInfeasible() ) {   
324        return NULL;
325    }
326
327    // The babSolver info used is the one containted in osi
328    OsiBabSolver * babSolver =
329        dynamic_cast<OsiBabSolver *> (osi->getAuxiliaryInfo());
330    //assert(babSolver == &babSolver_);
331    babSolver->setSolver(*bonmin_.nonlinearSolver()); 
332    //Last cut generator is used to check feasibility
333    bonmin_.cutGenerators().back().cgl->generateCuts(*osi, cuts_);
334    const int numvar = vars.size();
335    double* solverSol = new double[numvar];
336    double objValue = 1e200;
337   
338    BM_solution* sol = NULL;
339    if (babSolver->solution(objValue, solverSol,numvar)) {
340        sol = new BM_solution;
341        // Just copy the solution stored in solver to sol
342        for (int i = 0 ; i < numvar ; i++) {
343            if (solverSol[i] > lp_result.primalTolerance())
344                sol->add_entry(i, solverSol[i]); 
345        }
346        sol->setObjective(objValue);
347    }
348    delete[] solverSol;
349    return sol;
350}
351
352/****************************************************************************/
353
354void
355BM_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
356                           const BCP_vec<BCP_var*>& vars,
357                           const BCP_vec<BCP_cut*>& cuts,
358                           BCP_vec<BCP_cut*>& new_cuts,
359                           BCP_vec<BCP_row*>& new_rows)
360{
361    if (bonmin_.getAlgorithm() > 0) { /* if not pure B&B */
362        /* TODO:
363           don't invoke all of them, only the *good* ones. figure out
364           some measurement of how good a generator is.
365        */
366       
367        OsiSolverInterface& si = *getLpProblemPointer()->lp_solver;
368        double rand;
369  Bonmin::BabSetupBase::CuttingMethods::const_iterator end = bonmin_.cutGenerators().end();
370  end--;//have to go back one element because last cut generator checks feasibility
371   
372    /* FIXME: fill out the tree info! */
373    CglTreeInfo info;
374    for(Bonmin::BabSetupBase::CuttingMethods::const_iterator i = bonmin_.cutGenerators().begin() ;
375        i != end ; i++){
376    rand = 1 - CoinDrand48();
377    if(i->frequency > rand){
378      i->cgl->generateCuts(si, cuts_, info);
379    }
380  }
381       
382        // fill cuts
383        // eventually fill rows if not in strong branching
384        int numCuts = cuts_.sizeRowCuts();
385        for(int i = 0 ; i < numCuts ; i++) {
386            const OsiRowCut& cut = cuts_.rowCut(i);
387            new_cuts.push_back(new BB_cut(cut));
388            const CoinPackedVector& row = cut.row();
389            new_rows.push_back(new BCP_row(row, cut.lb(), cut.ub()));
390        }
391    }
392    cuts_.dumpCuts();
393}
394
395/****************************************************************************/
396
397void
398BM_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
399                    BCP_vec<BCP_cut*>& cuts,       // what to expand
400                    BCP_vec<BCP_row*>& rows,       // the expanded rows
401                    // things that the user can use for lifting cuts if allowed
402                    const BCP_lp_result& lpres,
403                    BCP_object_origin origin, bool allow_multiple)
404{
405    const int numCuts = cuts.size();
406    for (int i = 0; i < numCuts; ++i) {
407        BB_cut* cut = dynamic_cast<BB_cut*>(cuts[i]);
408        if (!cut) {
409            throw BCP_fatal_error("Non-\"BB_cut\" in cuts_to_rows!\n");
410        }
411        const CoinPackedVector& row = cut->row();
412        rows.push_back(new BCP_row(row, cuts[i]->lb(), cuts[i]->ub()));
413    }
414    cuts_.dumpCuts();
415}
416
417/****************************************************************************/
418
419double
420BM_lp::compute_lower_bound(const double old_lower_bound,
421                           const BCP_lp_result& lpres,
422                           const BCP_vec<BCP_var*>& vars,
423                           const BCP_vec<BCP_cut*>& cuts)
424{
425    double bd;
426    if (bonmin_.getAlgorithm() == 0) {
427        /* if pure B&B */
428        bd = lower_bound_;
429    } else {
430        bd = CoinMax(babSolver_.mipBound(), lpres.objval());
431    }
432    return CoinMax(bd, old_lower_bound);
433}
434
435/*---------------------------------------------------------------------------*/
Note: See TracBrowser for help on using the repository browser.