source: branches/0.2a/Bonmin/experimental/Bcp/BM_lp.cpp @ 528

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

Test a 0.2a branch

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("BM_lp: At node %i : will fathom because of high lower bound\n",
254                   current_index());
255        }
256    }
257    else if (nlp.isProvenPrimalInfeasible()) {
258        // prune it!
259        // FIXME: if nonconvex, restart from a different place...
260        lower_bound_ = 1e200;
261        numNlpFailed_ = 0;
262        if (par.entry(BM_par::PrintBranchingInfo)) {
263            printf("BM_lp: At node %i : will fathom because of infeasibility\n",
264                   current_index());
265        }
266    }
267    else if (nlp.isAbandoned()) {
268        if (nlp.isIterationLimitReached()) {
269            printf("\
270BM_lp: At node %i : WARNING: nlp reached iter limit. Will force branching\n",
271                   current_index());
272        } else {
273            printf("\
274BM_lp: At node %i : WARNING: nlp is abandoned. Will force branching\n",
275                   current_index());
276        }
277        // nlp failed
278        nlp.forceBranchable();
279        lower_bound_ = nlp.getObjValue();
280        CoinDisjointCopyN(colsol, numCols, primal_solution_);
281        numNlpFailed_ += 1;
282        // will branch, but save in the user data how many times we have
283        // failed, and if we fail too many times in a row then just abandon
284        // the node.
285    }
286    else {
287        // complain loudly about a bug
288        throw BCP_fatal_error("Impossible outcome by nlp.initialSolve()\n");
289    }
290    return sol;
291}
292
293/****************************************************************************/
294
295BCP_solution*
296BM_lp::test_feasibility_hybrid(const BCP_lp_result& lp_result,
297                               const BCP_vec<BCP_var*>& vars,
298                               const BCP_vec<BCP_cut*>& cuts)
299{
300    /* First test that the integrality requirements are met. */
301    BCP_solution_generic* gsol = test_full(lp_result, vars, integerTolerance_);
302    if (! gsol) {}
303       
304    /* TODO:
305       don't test feasibility in every node
306    */
307   
308
309    OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
310   
311    // Don't test feasibility if the node LP was infeasible
312    if (osi->isProvenPrimalInfeasible() ) {   
313        return NULL;
314    }
315
316    // The babSolver info used is the one containted in osi
317    OsiBabSolver * babSolver =
318        dynamic_cast<OsiBabSolver *> (osi->getAuxiliaryInfo());
319    //assert(babSolver == &babSolver_);
320    babSolver->setSolver(*bonmin_.nonlinearSolver()); 
321    //Last cut generator is used to check feasibility
322    bonmin_.cutGenerators().back().cgl->generateCuts(*osi, cuts_);
323    const int numvar = vars.size();
324    double* solverSol = new double[numvar];
325    double objValue = 1e200;
326   
327    BM_solution* sol = NULL;
328    if (babSolver->solution(objValue, solverSol,numvar)) {
329        sol = new BM_solution;
330        // Just copy the solution stored in solver to sol
331        for (int i = 0 ; i < numvar ; i++) {
332            if (solverSol[i] > lp_result.primalTolerance())
333                sol->add_entry(i, solverSol[i]); 
334        }
335        sol->setObjective(objValue);
336    }
337    delete[] solverSol;
338    return sol;
339}
340
341/****************************************************************************/
342
343void
344BM_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
345                           const BCP_vec<BCP_var*>& vars,
346                           const BCP_vec<BCP_cut*>& cuts,
347                           BCP_vec<BCP_cut*>& new_cuts,
348                           BCP_vec<BCP_row*>& new_rows)
349{
350    if (bonmin_.getAlgorithm() > 0) { /* if not pure B&B */
351        /* TODO:
352           don't invoke all of them, only the *good* ones. figure out
353           some measurement of how good a generator is.
354        */
355       
356        OsiSolverInterface& si = *getLpProblemPointer()->lp_solver;
357        double rand;
358  Bonmin::BabSetupBase::CuttingMethods::const_iterator end = bonmin_.cutGenerators().end();
359  end--;//have to go back one element because last cut generator checks feasibility
360   
361    /* FIXME: fill out the tree info! */
362    CglTreeInfo info;
363    for(Bonmin::BabSetupBase::CuttingMethods::const_iterator i = bonmin_.cutGenerators().begin() ;
364        i != end ; i++){
365    rand = 1 - CoinDrand48();
366    if(i->frequency > rand){
367      i->cgl->generateCuts(si, cuts_, info);
368    }
369  }
370       
371        // fill cuts
372        // eventually fill rows if not in strong branching
373        int numCuts = cuts_.sizeRowCuts();
374        for(int i = 0 ; i < numCuts ; i++) {
375            const OsiRowCut& cut = cuts_.rowCut(i);
376            new_cuts.push_back(new BB_cut(cut));
377            const CoinPackedVector& row = cut.row();
378            new_rows.push_back(new BCP_row(row, cut.lb(), cut.ub()));
379        }
380    }
381    cuts_.dumpCuts();
382}
383
384/****************************************************************************/
385
386void
387BM_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
388                    BCP_vec<BCP_cut*>& cuts,       // what to expand
389                    BCP_vec<BCP_row*>& rows,       // the expanded rows
390                    // things that the user can use for lifting cuts if allowed
391                    const BCP_lp_result& lpres,
392                    BCP_object_origin origin, bool allow_multiple)
393{
394    const int numCuts = cuts.size();
395    for (int i = 0; i < numCuts; ++i) {
396        BB_cut* cut = dynamic_cast<BB_cut*>(cuts[i]);
397        if (!cut) {
398            throw BCP_fatal_error("Non-\"BB_cut\" in cuts_to_rows!\n");
399        }
400        const CoinPackedVector& row = cut->row();
401        rows.push_back(new BCP_row(row, cuts[i]->lb(), cuts[i]->ub()));
402    }
403    cuts_.dumpCuts();
404}
405
406/****************************************************************************/
407
408double
409BM_lp::compute_lower_bound(const double old_lower_bound,
410                           const BCP_lp_result& lpres,
411                           const BCP_vec<BCP_var*>& vars,
412                           const BCP_vec<BCP_cut*>& cuts)
413{
414    double bd;
415    if (bonmin_.getAlgorithm() == 0) {
416        /* if pure B&B */
417        bd = lower_bound_;
418    } else {
419        bd = CoinMax(babSolver_.mipBound(), lpres.objval());
420    }
421    return CoinMax(bd, old_lower_bound);
422}
423
424/*---------------------------------------------------------------------------*/
Note: See TracBrowser for help on using the repository browser.