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

Last change on this file since 508 was 508, checked in by pbonami, 12 years ago

Make bcp-bonmin work with new setup

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