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

Last change on this file since 530 was 530, checked in by ladanyi, 12 years ago

prepare Bcp for changing some LPs into TM storage places. Serial version works, currently debugging parallel version.

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