/home/coin/SVN-release/OS-2.4.1/Bonmin/experimental/Bcp/BM_lp.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation 2006, 2007
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Laszlo Ladanyi, International Business Machines Corporation
00007 // Pierre Bonami, Carnegie Mellon University
00008 
00009 #include "OsiClpSolverInterface.hpp"
00010 #include "BM.hpp"
00011 #include "BCP_message_mpi.hpp"
00012 #include "BCP_lp_node.hpp"
00013 #include "BCP_lp.hpp"
00014 
00015 #include "BonOACutGenerator2.hpp"
00016 #include "BonEcpCuts.hpp"
00017 #include "BonOaNlpOptim.hpp"
00018 
00019 // The following is included for "min"
00020 #include "CoinFinite.hpp"
00021 
00022 #ifndef BM_DEBUG_PRINT
00023 #define BM_DEBUG_PRINT 0
00024 #endif
00025 
00026 static char prefix[100];
00027 
00028 //#############################################################################
00029 
00030 BM_lp::BM_lp() :
00031     BCP_lp_user(),
00032     in_strong(0),
00033     bonmin_(),
00034     objInd_(NULL),
00035     objNum_(0),
00036     infInd_(NULL),
00037     infUseful_(NULL),
00038     infNum_(0),
00039     feasInd_(NULL),
00040     feasUseful_(NULL),
00041     feasNum_(0),
00042     sbResult_(NULL),
00043     bestSbResult_(NULL)
00044 {
00045 }
00046 
00047 /****************************************************************************/
00048 
00049 BM_lp::~BM_lp()
00050 {
00051   delete[] objInd_;
00052   delete[] infInd_;
00053   delete[] infUseful_;
00054   delete[] feasInd_;
00055   delete[] feasUseful_;
00056   delete[] sbResult_;
00057 }
00058 
00059 /****************************************************************************/
00060 
00061 OsiSolverInterface *
00062 BM_lp::initialize_solver_interface()
00063 {
00064 #ifdef COIN_HAS_MPI
00065   sprintf(prefix, "%i", getLpProblemPointer()->get_process_id());
00066 #else
00067   prefix[0] = 0;
00068 #endif
00069 
00070   OsiSolverInterface* solver = NULL;
00071   if (bonmin_.getAlgorithm() == 0) {
00072     // Pure B&B
00073     solver = bonmin_.nonlinearSolver()->clone();
00074   } else {
00075     OsiClpSolverInterface * clp = new OsiClpSolverInterface;
00076     OsiBabSolver babSolver(3);
00077     babSolver.setSolver(clp);
00078     clp->setAuxiliaryInfo(&babSolver);
00079     clp->messageHandler()->setLogLevel(0);
00080     setOsiBabSolver(dynamic_cast<OsiBabSolver *>(clp->getAuxiliaryInfo()));
00081     solver = clp;
00082   }
00083   return solver;
00084 }
00085 
00086 /****************************************************************************/
00087 
00088 void
00089 BM_lp::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
00090                                        const BCP_vec<BCP_cut*>& cuts,
00091                                        const BCP_vec<BCP_obj_status>& vs,
00092                                        const BCP_vec<BCP_obj_status>& cs,
00093                                        BCP_vec<int>& var_changed_pos,
00094                                        BCP_vec<double>& var_new_bd,
00095                                        BCP_vec<int>& cut_changed_pos,
00096                                        BCP_vec<double>& cut_new_bd)
00097 {
00098     node_start_time = CoinWallclockTime();
00099     BM_node* data = dynamic_cast<BM_node*>(get_user_data());
00100     numNlpFailed_ = data->numNlpFailed_;
00101 
00102     if (bonmin_.getAlgorithm() != 0) {
00103       // Not pure BB, so an LP solver will be used. Now we have to...
00104 
00105       // First copy the bounds into nlp. That way all the branching decisions
00106       // will be transferred over.
00107       int i;
00108 
00109       OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
00110       Bonmin::OsiTMINLPInterface& nlp = *bonmin_.nonlinearSolver();
00111       nlp.setColLower(osi->getColLower());
00112       nlp.setColUpper(osi->getColUpper());
00113 
00114       // Carry the changes over to the object lists in nlp
00115       const int numObj = nlp.numberObjects();
00116       OsiObject** nlpObj = nlp.objects();
00117       for (i = 0; i < numObj; ++i) {
00118         OsiSimpleInteger* io = dynamic_cast<OsiSimpleInteger*>(nlpObj[i]);
00119         if (io) {
00120           io->resetBounds(&nlp);
00121         } else {
00122           // The rest is OsiSOS where we don't need to do anything
00123           break;
00124         }
00125       }
00126 
00127       // copy over the OsiObjects from the nlp solver
00128       osi->addObjects(nlp.numberObjects(), nlp.objects());
00129     }
00130 
00131     in_strong = 0;
00132 }
00133 
00134 /************************************************************************/
00135 void
00136 BM_lp::load_problem(OsiSolverInterface& osi, BCP_problem_core* core,
00137                     BCP_var_set& vars, BCP_cut_set& cuts)
00138 {
00139   if (bonmin_.getAlgorithm() != 0) {
00140     // We are doing hybrid, so osi is an LP solver. Call the default.
00141     BCP_lp_user::load_problem(osi, core, vars, cuts);
00142     return;
00143   }
00144   // Otherwise we do B&B and osi is an NLP solver.
00145   // There is no need to fill it with the data from bonmin_.nonlinearSolver()
00146   // since osi is a clone() of the master_lp in the BCP_lp object, and the
00147   // master_lp is a clone() of bonmin_.nonlinearSolver(), and that clone()
00148   // copies the data as well.
00149 }
00150 
00151 /************************************************************************/
00152 void
00153 BM_lp::modify_lp_parameters(OsiSolverInterface* lp, bool in_strong_branching)
00154     // Called each time the node LP is solved
00155 {
00156   if (in_strong_branching) {
00157     in_strong = 1;
00158     //     lp->setIntParam(OsiMaxNumIterationHotStart, 50);
00159   }
00160 }
00161 
00162 /****************************************************************************/
00163 
00164 BCP_solution*
00165 BM_lp::test_feasibility(const BCP_lp_result& lp_result,
00166                         const BCP_vec<BCP_var*>& vars,
00167                         const BCP_vec<BCP_cut*>& cuts)
00168 {
00169   if (bonmin_.getAlgorithm() == 0) {
00170     // if pure B&B
00171     return test_feasibility_BB(lp_result, vars);
00172   } else {
00173     return test_feasibility_hybrid(lp_result, vars, cuts);
00174   }
00175   return NULL; // fake return to quiet gcc
00176 }
00177 
00178 /****************************************************************************/
00179 
00180 BCP_solution*
00181 BM_lp::test_feasibility_BB(const BCP_lp_result& lpres,
00182                            const BCP_vec<BCP_var*>& vars)
00183 {
00184   // We can just take the primal solution and test whether it satisfies
00185   // integrality requirements
00186   BCP_solution_generic* sol = test_full(lpres, vars, integerTolerance_);
00187   if (sol) {
00188     sol->set_objective_value(lpres.objval());
00189 #if (BM_DEBUG_PRINT != 0)
00190     printf("LP %.3f: Solution found. node: %i  depth: %i  value: %f\n",
00191            CoinWallclockTime() - start_time(),
00192            current_index(), current_level(), lpres.objval());
00193 #endif
00194   }
00195   return sol;
00196 }
00197 
00198 /****************************************************************************/
00199 
00200 BCP_solution*
00201 BM_lp::test_feasibility_hybrid(const BCP_lp_result& lp_result,
00202                                const BCP_vec<BCP_var*>& vars,
00203                                const BCP_vec<BCP_cut*>& cuts)
00204 {
00205     /* First test that the integrality requirements are met. */
00206     BCP_solution_generic* gsol = test_full(lp_result, vars, integerTolerance_);
00207     if (! gsol) {}
00208         
00209     /* TODO:
00210        don't test feasibility in every node
00211     */
00212     
00213 
00214     OsiSolverInterface * osi = getLpProblemPointer()->lp_solver;
00215     
00216     // Don't test feasibility if the node LP was infeasible
00217     if (osi->isProvenPrimalInfeasible() ) {   
00218         return NULL;
00219     }
00220 
00221     // The babSolver info used is the one containted in osi
00222     OsiBabSolver * babSolver =
00223         dynamic_cast<OsiBabSolver *> (osi->getAuxiliaryInfo());
00224     babSolver->setSolver(*bonmin_.nonlinearSolver()); 
00225     //Last cut generator is used to check feasibility
00226     bonmin_.cutGenerators().back().cgl->generateCuts(*osi, cuts_);
00227     const int numvar = vars.size();
00228     double* solverSol = new double[numvar];
00229     double objValue = 1e200;
00230     
00231     BCP_solution_generic* sol = NULL;
00232     if (babSolver->solution(objValue, solverSol, numvar)) {
00233         sol = new BCP_solution_generic(false);
00234         // Just copy the solution stored in solver to sol
00235         for (int i = 0 ; i < numvar ; i++) {
00236             if (solverSol[i] > lp_result.primalTolerance())
00237                 sol->add_entry(vars[i], solverSol[i]); 
00238         }
00239         sol->set_objective_value(objValue);
00240     }
00241     delete[] solverSol;
00242     return sol;
00243 }
00244 
00245 /****************************************************************************/
00246 
00247 void
00248 BM_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
00249                            const BCP_vec<BCP_var*>& vars,
00250                            const BCP_vec<BCP_cut*>& cuts,
00251                            BCP_vec<BCP_cut*>& new_cuts,
00252                            BCP_vec<BCP_row*>& new_rows)
00253 {
00254     if (bonmin_.getAlgorithm() > 0) { /* if not pure B&B */
00255         /* TODO:
00256            don't invoke all of them, only the *good* ones. figure out
00257            some measurement of how good a generator is.
00258         */
00259         
00260         OsiSolverInterface& si = *getLpProblemPointer()->lp_solver;
00261         double rand;
00262   Bonmin::BabSetupBase::CuttingMethods::const_iterator end = bonmin_.cutGenerators().end();
00263   end--;//have to go back one element because last cut generator checks feasibility
00264     
00265     /* FIXME: fill out the tree info! */
00266     CglTreeInfo info;
00267     for(Bonmin::BabSetupBase::CuttingMethods::const_iterator i = bonmin_.cutGenerators().begin() ;
00268         i != end ; i++){
00269     rand = 1 - CoinDrand48();
00270     if(i->frequency > rand){
00271       i->cgl->generateCuts(si, cuts_, info);
00272     }
00273   }
00274         
00275         // fill cuts
00276         // eventually fill rows if not in strong branching
00277         int numCuts = cuts_.sizeRowCuts();
00278         for(int i = 0 ; i < numCuts ; i++) {
00279             const OsiRowCut& cut = cuts_.rowCut(i);
00280             new_cuts.push_back(new BB_cut(cut));
00281             const CoinPackedVector& row = cut.row();
00282             new_rows.push_back(new BCP_row(row, cut.lb(), cut.ub()));
00283         }
00284     }
00285     cuts_.dumpCuts();
00286 }
00287 
00288 /****************************************************************************/
00289 
00290 void
00291 BM_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
00292                     BCP_vec<BCP_cut*>& cuts,       // what to expand
00293                     BCP_vec<BCP_row*>& rows,       // the expanded rows
00294                     // things that the user can use for lifting cuts if allowed
00295                     const BCP_lp_result& lpres,
00296                     BCP_object_origin origin, bool allow_multiple)
00297 {
00298     const int numCuts = cuts.size();
00299     for (int i = 0; i < numCuts; ++i) {
00300         BB_cut* cut = dynamic_cast<BB_cut*>(cuts[i]);
00301         if (!cut) {
00302             throw BCP_fatal_error("Non-\"BB_cut\" in cuts_to_rows!\n");
00303         }
00304         const CoinPackedVector& row = cut->row();
00305         rows.push_back(new BCP_row(row, cuts[i]->lb(), cuts[i]->ub()));
00306     }
00307     cuts_.dumpCuts();
00308 }
00309 
00310 /****************************************************************************/
00311 
00312 double
00313 BM_lp::compute_lower_bound(const double old_lower_bound,
00314                            const BCP_lp_result& lpres,
00315                            const BCP_vec<BCP_var*>& vars,
00316                            const BCP_vec<BCP_cut*>& cuts)
00317 {
00318     double bd;
00319     if (bonmin_.getAlgorithm() == 0) {
00320         /* if pure B&B */
00321         bd = lpres.objval();
00322     } else {
00323       /* FIXME: what the heck was this...
00324         bd = CoinMax(babSolver_.mipBound(), lpres.objval());
00325       */
00326     }
00327     return CoinMax(bd, old_lower_bound);
00328 }
00329 
00330 /*---------------------------------------------------------------------------*/

Generated on Thu Nov 10 03:05:39 2011 by  doxygen 1.4.7