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

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

finer output control

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