source: stable/0.1/Bonmin/experimental/Bcp/BM_lp.cpp @ 52

Last change on this file since 52 was 52, checked in by ladanyi, 13 years ago

Added warmstart

File size: 13.0 KB
Line 
1#include "OsiClpSolverInterface.hpp"
2#include "BM.hpp"
3
4//#############################################################################
5
6BM_lp::BM_lp() :
7    BCP_lp_user(),
8    sos(),
9    babSolver_(3),
10    nlp(),
11    ws(NULL),
12    feasChecker_(NULL),
13    in_strong(0)
14{
15    babSolver_.setSolver(&nlp);
16    setOsiBabSolver(&babSolver_);
17}
18
19/****************************************************************************/
20
21BM_lp::~BM_lp()
22{
23    delete ws;
24    delete feasChecker_;
25    delete[] primal_solution_;
26    delete[] branching_priority_;
27    for (int i = sos.size() - 1; i >= 0; --i) {
28        delete sos[i];
29    }
30}
31
32/****************************************************************************/
33
34OsiSolverInterface *
35BM_lp::initialize_solver_interface()
36{
37    OsiClpSolverInterface * clp = new OsiClpSolverInterface;
38    OsiBabSolver babSolver(3);
39    babSolver.setSolver(clp);
40    clp->setAuxiliaryInfo(&babSolver);
41    clp->messageHandler()->setLogLevel(0);
42    setOsiBabSolver(dynamic_cast<OsiBabSolver *>(clp->getAuxiliaryInfo()));
43    return clp;
44}
45
46/****************************************************************************/
47
48void
49BM_lp::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
50                                       const BCP_vec<BCP_cut*>& cuts,
51                                       const BCP_vec<BCP_obj_status>& vs,
52                                       const BCP_vec<BCP_obj_status>& cs,
53                                       BCP_vec<int>& var_changed_pos,
54                                       BCP_vec<double>& var_new_bd,
55                                       BCP_vec<int>& cut_changed_pos,
56                                       BCP_vec<double>& cut_new_bd)
57{
58    int i;
59    // First copy the bounds into nlp. That way all the branching decisions
60    // will be transferred over.
61    for (i = sos.size() - 1; i >= 0; --i) {
62        sos[i]->first = 0;
63        sos[i]->last = sos[i]->length;
64    }
65    OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
66    const int numCols = osi->getNumCols();
67    const double* clb = osi->getColLower();
68    const double* cub = osi->getColUpper();
69    for (i = 0; i < numCols; ++i) {
70        const BCP_var_core* v =
71            dynamic_cast<const BCP_var_core*>(vars[i]);
72        if (v) {
73            nlp.setColLower(i, clb[i]);
74            nlp.setColUpper(i, cub[i]);
75            continue;
76        }
77        const BM_branching_var* bv =
78            dynamic_cast<const BM_branching_var*>(vars[i]);
79        if (bv) {
80            const int ind = bv->sos_index;
81            const int split = bv->split;
82            const char type = sos[ind]->type;
83            const int *indices = sos[ind]->indices;
84            if (bv->ub() == 0.0) {
85                const int last = sos[ind]->last;
86                for (int j = split + type; j < last; ++j) {
87                    nlp.setColLower(indices[j], 0.0);
88                    nlp.setColUpper(indices[j], 0.0);
89                }
90                sos[ind]->last = split + type;
91            } else {
92                const int first = sos[ind]->first;
93                for (int j = first; j < split; ++j) {
94                    nlp.setColLower(indices[j], 0.0);
95                    nlp.setColUpper(indices[j], 0.0);
96                }
97                sos[ind]->first = split;
98            }
99            continue;
100        }
101        throw BCP_fatal_error("BM_lp: unrecognized variable type\n");
102    }
103
104    in_strong = 0;
105}
106
107/************************************************************************/
108void
109BM_lp::modify_lp_parameters(OsiSolverInterface* lp, bool in_strong_branching)
110
111    // Called each time the node LP is solved
112
113{
114    if (in_strong_branching) {
115        in_strong = 1;
116        //     lp->setIntParam(OsiMaxNumIterationHotStart, 50);
117    }
118}
119
120/****************************************************************************/
121
122BCP_solution*
123BM_lp::test_feasibility(const BCP_lp_result& lp_result,
124                        const BCP_vec<BCP_var*>& vars,
125                        const BCP_vec<BCP_cut*>& cuts)
126{
127    BM_solution* sol = NULL;
128
129    if (par.entry(BM_par::PureBranchAndBound)) {
130        /* PURE_BB TODO:
131           Do whatever must be done to:
132           1) compute a lower bound using NLP and save it in "lower_bound_"
133           2) save the solution of the NLP in "primal_solution_"
134           3) if the optimum (well, local opt anyway...) happens to be
135           feasible (or along the NLP optimization we have blundered into
136           a feas sol) then create "sol"
137        */
138        if (current_index() == 0) {
139            nlp.initialSolve();
140            switch (par.entry(BM_par::WarmStartStrategy)) {
141            case WarmStartNone:
142                break;
143            case WarmStartFromRoot:
144                ws = nlp.getWarmStart();
145                break;
146            case WarmStartFromParent:
147                // FIXME: not yet implemented
148                break;
149            }
150        } else {
151            switch (par.entry(BM_par::WarmStartStrategy)) {
152            case WarmStartNone:
153                nlp.initialSolve();
154                break;
155            case WarmStartFromRoot:
156                nlp.setWarmStart(ws);
157                nlp.resolve();
158                break;
159            case WarmStartFromParent:
160                // FIXME: not yet implemented
161                nlp.initialSolve();
162                break;
163            }
164        }
165        if (nlp.isProvenOptimal()) {
166            const int numCols = nlp.getNumCols();
167            lower_bound_ = nlp.getObjValue();
168            CoinDisjointCopyN(nlp.getColSolution(), numCols, primal_solution_);
169            numNlpFailed_ = 0;
170            Ipopt::SmartPtr<Ipopt::OptionsList> options = nlp.retrieve_options();
171            double intTol;
172            options->GetNumericValue("integer_tolerance",intTol,"bonmin.");
173
174            int i;
175            for (i = 0; i < numCols; ++i) {
176                if (vars[i]->var_type() == BCP_ContinuousVar)
177                    continue;
178                const double psol = CoinMin(CoinMax(vars[i]->lb(),
179                                                    primal_solution_[i]),
180                                            vars[i]->ub());
181                const double frac = psol - floor(psol);
182                const double dist = std::min(frac, 1.0-frac);
183                if (dist >= intTol)
184                    break;
185            }
186            if (i == numCols) {
187                /* yipee! a feasible solution! */
188                sol = new BM_solution;
189                //Just copy the solution stored in solver to sol
190                for (i = 0 ; i < numCols ; i++) {
191                    if (primal_solution_[i] > lp_result.primalTolerance())
192                        sol->add_entry(i, primal_solution_[i]); 
193                }
194                sol->setObjective(lower_bound_);
195            }
196            if (lower_bound_>upper_bound()-get_param(BCP_lp_par::Granularity) &&
197                par.entry(BM_par::PrintBranchingInfo)) {
198                printf("\
199BM_lp: At node %i : will fathom because of high lower bound\n",
200                       current_index());
201            }
202        }
203        else if (nlp.isProvenPrimalInfeasible()) {
204            // prune it!
205            lower_bound_ = 1e200;
206            numNlpFailed_ = 0;
207            if (par.entry(BM_par::PrintBranchingInfo)) {
208                printf("\
209BM_lp: At node %i : will fathom because of infeasibility\n",
210                       current_index());
211            }
212        }
213        else if (nlp.isAbandoned()) {
214            if (nlp.isIterationLimitReached()) {
215                printf("\
216BM_lp: At node %i : WARNING: nlp reached iter limit. Will force branching\n",
217                       current_index());
218            } else {
219                printf("\
220BM_lp: At node %i : WARNING: nlp is abandoned. Will force branching\n",
221                       current_index());
222            }
223            // nlp failed
224            nlp.forceBranchable();
225            lower_bound_ = nlp.getObjValue();
226            CoinDisjointCopyN(nlp.getColSolution(), vars.size(),
227                              primal_solution_);
228            numNlpFailed_ += 1;
229            // will branch, but save in the user data how many times we have
230            // failed, and if we fail too many times in a row then just abandon
231            // the node.
232        }
233        else {
234            // complain loudly about a bug
235            throw BCP_fatal_error("Impossible outcome by nlp.initialSolve()\n");
236        }
237    }
238    else {
239        /* TODO:
240           don't test feasibility in every node
241        */
242        double integerTolerance;
243        double cutOffIncrement;
244        nlp.retrieve_options()->GetNumericValue("integer_tolerance",
245                                                integerTolerance,"bonmin.");
246        /* FIXME: cutoff_decr was called cutoff_incr??? */
247        nlp.retrieve_options()->GetNumericValue("cutoff_decr",
248                                                cutOffIncrement,"bonmin.");
249        OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
250   
251        // Don't test feasibility if the node LP was infeasible
252        if (osi->isProvenPrimalInfeasible() ) {   
253            return sol;
254        }
255
256        if (!feasChecker_) {
257            feasChecker_ = new IpCbcOACutGenerator2(&nlp, NULL, NULL, 
258                                                    cutOffIncrement,
259                                                    integerTolerance, 0, 1);
260            feasChecker_->setLocalSearchNodeLimit(0);
261            feasChecker_->setMaxLocalSearch(0);
262            feasChecker_->setMaxLocalSearchPerNode(0);
263        }
264
265        // The babSolver info used is the one containted in osi
266        OsiBabSolver * babSolver =
267            dynamic_cast<OsiBabSolver *> (osi->getAuxiliaryInfo());
268        //assert(babSolver == &babSolver_);
269        babSolver->setSolver(nlp); 
270        feasChecker_->generateCuts(*osi, cuts_);
271        const int numvar = vars.size();
272        double* solverSol = new double[numvar];
273        double objValue = 1e200;
274   
275        if (babSolver->solution(objValue, solverSol,numvar)) {
276            sol = new BM_solution;
277            //Just copy the solution stored in solver to sol
278            for (int i = 0 ; i < numvar ; i++) {
279                if (solverSol[i] > lp_result.primalTolerance())
280                    sol->add_entry(i, solverSol[i]); 
281            }
282            sol->setObjective(objValue);
283        }
284        delete[] solverSol;
285    }
286    return sol;
287}
288
289/****************************************************************************/
290
291void
292BM_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
293                           const BCP_vec<BCP_var*>& vars,
294                           const BCP_vec<BCP_cut*>& cuts,
295                           BCP_vec<BCP_cut*>& new_cuts,
296                           BCP_vec<BCP_row*>& new_rows)
297{
298    if (! par.entry(BM_par::PureBranchAndBound)) {
299        /* TODO:
300           generate cuts with various other Cgl methods (Gomory, etc)
301        */
302        // fill cuts
303        // eventually fill rows if not in strong branching
304        int numCuts = cuts_.sizeRowCuts();
305        for(int i = 0 ; i < numCuts ; i++) {
306            const OsiRowCut& cut = cuts_.rowCut(i);
307            new_cuts.push_back(new BB_cut(cut));
308            const CoinPackedVector& row = cut.row();
309            new_rows.push_back(new BCP_row(row, cut.lb(), cut.ub()));
310        }
311    }
312    cuts_.dumpCuts();
313}
314
315/****************************************************************************/
316
317void
318BM_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
319                    BCP_vec<BCP_cut*>& cuts,       // what to expand
320                    BCP_vec<BCP_row*>& rows,       // the expanded rows
321                    // things that the user can use for lifting cuts if allowed
322                    const BCP_lp_result& lpres,
323                    BCP_object_origin origin, bool allow_multiple)
324{
325    const int numCuts = cuts.size();
326    for (int i = 0; i < numCuts; ++i) {
327        BB_cut* cut = dynamic_cast<BB_cut*>(cuts[i]);
328        if (!cut) {
329            throw BCP_fatal_error("Non-\"BB_cut\" in cuts_to_rows!\n");
330        }
331        const CoinPackedVector& row = cut->row();
332        rows.push_back(new BCP_row(row, cuts[i]->lb(), cuts[i]->ub()));
333    }
334}
335
336/****************************************************************************/
337
338void
339BM_lp::vars_to_cols(const BCP_vec<BCP_cut*>& cuts,
340                    BCP_vec<BCP_var*>& vars,
341                    BCP_vec<BCP_col*>& cols,
342                    const BCP_lp_result& lpres,
343                    BCP_object_origin origin, bool allow_multiple)
344{
345    /* All vars better be BM_branching_var, and their columns are empty */
346    const int numVars = vars.size();
347    for (int i = 0; i < numVars; ++i) {
348        BM_branching_var* bv = dynamic_cast<BM_branching_var*>(vars[i]);
349        if (!bv) {
350            throw BCP_fatal_error("Non-\"BM_branching_var\" in vars_to_cols!\n");
351        }
352        cols.push_back(new BCP_col());
353    }
354}
355
356/****************************************************************************/
357
358double
359BM_lp::compute_lower_bound(const double old_lower_bound,
360                           const BCP_lp_result& lpres,
361                           const BCP_vec<BCP_var*>& vars,
362                           const BCP_vec<BCP_cut*>& cuts)
363{
364    double bd;
365    if (par.entry(BM_par::PureBranchAndBound)) {
366        bd = lower_bound_;
367    } else {
368        bd = std::max(babSolver_.mipBound(), lpres.objval());
369    }
370    return std::max(bd, old_lower_bound);
371}
372
373/****************************************************************************/
374
375BM_lp::BmSosInfo::BmSosInfo(const TMINLP::SosInfo * sosinfo, int i)
376{
377    priority = sosinfo->priorities ? sosinfo->priorities[i] : 0;
378    length   = sosinfo->starts[i+1] - sosinfo->starts[i];
379    type     = sosinfo->types[i] == 1 ? 0 : 1;
380    indices  = new int[length];
381    weights  = new double[length];
382    first    = 0;
383    last     = length;
384    CoinDisjointCopyN(sosinfo->indices+sosinfo->starts[i], length, indices);
385    CoinDisjointCopyN(sosinfo->weights+sosinfo->starts[i], length, weights);
386    CoinSort_2(weights, weights+length, indices);
387}
388   
389/*---------------------------------------------------------------------------*/
390
391void BM_lp::setSosFrom(const TMINLP::SosInfo * sosinfo)
392{
393    if (!sosinfo || sosinfo->num == 0)
394        return;
395
396    //we have some sos constraints
397    sos.reserve(sosinfo->num);
398    for (int i = 0; i < sosinfo->num; ++i) {
399        sos.push_back(new BmSosInfo(sosinfo, i));
400    }
401
402    if (sosinfo->priorities) {
403        int * priorities = new int[sosinfo->num];
404        CoinDisjointCopyN(sosinfo->priorities, sosinfo->num, priorities);
405        /* FIXME: we may need to go down in order! */
406        if (par.entry(BM_par::SosWithLowPriorityMoreImportant)) {
407            CoinSort_2(priorities, priorities+sosinfo->num, &sos[0],
408                       CoinFirstLess_2<int,BM_lp::BmSosInfo*>());
409        } else {
410            CoinSort_2(priorities, priorities+sosinfo->num, &sos[0],
411                       CoinFirstGreater_2<int,BM_lp::BmSosInfo*>());
412        }
413        delete[] priorities;
414    }
415}
416
417/*---------------------------------------------------------------------------*/
418
419void BM_lp::process_message(BCP_buffer& buf)
420{
421}
422
423/*---------------------------------------------------------------------------*/
424
425#if 0
426void BM_lp::BmSosInfo::shuffle()
427{
428    if (num == 0)
429        return;
430    int pos = 0;
431    while (pos < num) {
432        int cnt = 0;
433        while (priorities[priority_order[pos]] ==
434               priorities[priority_order[pos+cnt]]) cnt++;
435        if (cnt > 1) {
436            std::random_shuffle(priority_order+pos, priority_order+(pos+cnt));
437        }
438        pos += cnt;
439    }
440}
441#endif
Note: See TracBrowser for help on using the repository browser.