/home/coin/SVN-release/OS-2.4.2/Bcp/examples/BAC/LP/BB_lp.cpp

Go to the documentation of this file.
00001 // Last edit: 1/21/09
00002 //
00003 // Name:     BB_lp.cpp
00004 // Author:   Francois Margot
00005 //           Tepper School of Business
00006 //           Carnegie Mellon University, Pittsburgh, PA 15213
00007 //           email: fmargot@andrew.cmu.edu
00008 // Date:     12/28/03
00009 //-----------------------------------------------------------------------------
00010 // Copyright (C) 2003, Francois Margot, International Business Machines
00011 // Corporation and others.  All Rights Reserved.
00012 
00013 #include <vector>
00014 
00015 #include "BB_lp.hpp"
00016 #include "BB_cut.hpp"
00017 #include "OsiClpSolverInterface.hpp"
00018 
00019 #include "BCP_math.hpp"
00020 #include "BCP_lp.hpp"
00021 #include "BCP_problem_core.hpp"
00022 
00023 using namespace std;
00024 
00025 /************************************************************************/
00026 void
00027 BB_lp::unpack_module_data(BCP_buffer& buf)
00028 {
00029   buf.unpack(p_desc);
00030   EPS = p_desc->EPSILON;
00031 }
00032 
00033 /************************************************************************/
00034 OsiSolverInterface *
00035 BB_lp::initialize_solver_interface()
00036 
00037   // Called once at the beginning, from the root node
00038 
00039 {
00040   OsiClpSolverInterface * clp = new OsiClpSolverInterface;
00041   clp->messageHandler()->setLogLevel(0);
00042   clp->getModelPtr()->messageHandler()->setLogLevel(0);
00043 
00044   return clp;
00045 }
00046 
00047 /************************************************************************/
00048 void
00049 BB_lp::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
00050                                        const BCP_vec<BCP_cut*>& cuts,
00051                                        const BCP_vec<BCP_obj_status>& vstatus,
00052                                        const BCP_vec<BCP_obj_status>& cstatus,
00053                                        BCP_vec<int>& var_changed_pos,
00054                                        BCP_vec<double>& var_new_bd,
00055                                        BCP_vec<int>& cut_changed_pos,
00056                                        BCP_vec<double>& cut_new_bd)
00057 
00058   // Called at the beginning of the optimization of a node. 
00059   // The node LP is not yet solved.
00060 
00061 {
00062   in_strong = 0;
00063 
00064 #ifdef USER_DATA
00065   MY_user_data *curr_ud = dynamic_cast<MY_user_data*> (get_user_data());
00066   curr_ud->is_processed = 1;
00067 #endif
00068 
00069 }
00070 
00071 /************************************************************************/
00072 void
00073 BB_lp::modify_lp_parameters(OsiSolverInterface* lp, const int changeType,
00074                             bool in_strong_branching)
00075 
00076   // Called each time the node LP is solved
00077 
00078 {
00079    if (in_strong_branching) {
00080      in_strong = 1;
00081      lp->setIntParam(OsiMaxNumIterationHotStart, 50);
00082    }
00083 
00084    // write the current LP in file lpnode.lp
00085    lp->writeLp("lpnode", "lp");
00086    cout << "LP node written in file lpnode.lp" << endl;
00087 
00088    // to write the current LP file in file lpnode.mps use the following:
00089    // lp->writeMps("lpnode", "mps", 0.0);
00090    // cout << "LP node written in file lpnode.mps" << endl;
00091 
00092 }
00093 
00094 /************************************************************************/
00095 BCP_solution*
00096 BB_lp::test_feasibility(const BCP_lp_result& lp_result,
00097                         const BCP_vec<BCP_var*>& vars,
00098                         const BCP_vec<BCP_cut*>& cuts)
00099 
00100   // Called after each node LP has been solved.
00101   // Called even if the node LP was infeasible
00102   // Called also during strong branching
00103 
00104 {
00105   // Don't test feasibility during strong branching
00106 
00107   if (in_strong) {
00108     return(0);
00109   }
00110   
00111   // Don't test feasibility if the node LP was infeasible or
00112   // termination was not clean
00113 
00114   if(getLpProblemPointer()->lp_solver->isAbandoned() ||
00115      getLpProblemPointer()->lp_solver->isProvenPrimalInfeasible() ||
00116      getLpProblemPointer()->lp_solver->isDualObjectiveLimitReached() ||
00117      getLpProblemPointer()->lp_solver->isIterationLimitReached()) {
00118     return(0);
00119   }
00120 
00121   int i, j, k, cn = p_desc->colnum;;
00122   const double *x = lp_result.x();
00123   double *check_lhs = new double[p_desc->indexed->getNumRows()];
00124 
00125   // check indexed constraints
00126 
00127   violated_cuts.clear();
00128   p_desc->indexed->times(lp_result.x(), check_lhs);
00129   for (i=0; i<p_desc->indexed->getNumRows(); ++i) {
00130     if ((check_lhs[i] < p_desc->rlb_indexed[i] - EPS) || 
00131         (check_lhs[i] > p_desc->rub_indexed[i] + EPS))
00132       violated_cuts.push_back(i);
00133   }
00134 
00135   delete[] check_lhs;
00136 
00137   OsiRowCut* rcut = new OsiRowCut();
00138   int *cutind = new int[cn], cut_nz;
00139   double* cutcoef = new double[cn], cutrhs = 1;
00140 
00141   // check algorithmic cuts
00142 
00143   for(i=0; i<cn; i++) {
00144     j = (i+1) % cn;
00145     k = (i+2) % cn;
00146 
00147     cutcoef[0] = 1;
00148     cutcoef[1] = 1;
00149     cutcoef[2] = 1;
00150     cut_nz = 3;
00151 
00152     if(x[i] + x[j] + x[k] > 1 + EPS) {
00153 
00154       // cut is violated 
00155 
00156       cutind[0] = i;
00157       cutind[1] = j;
00158       cutind[2] = k;
00159       
00160       rcut->setLb(-BCP_DBL_MAX);
00161       rcut->setUb(cutrhs);
00162       rcut->setRow(cut_nz, cutind, cutcoef);
00163       
00164       BB_cut* cut = new BB_cut(*rcut);
00165       algo_cuts.push_back(cut);
00166       
00167     }
00168   }
00169   
00170   delete rcut;
00171   delete[] cutind;
00172   delete[] cutcoef;
00173 
00174   // if all constraints are satisfied, check integrality of the vars
00175 
00176   return (violated_cuts.empty() +  algo_cuts.empty() == 2 ?
00177           BCP_lp_user::test_feasibility(lp_result, vars, cuts) : 0);
00178 }
00179 
00180 /********************************************************************/
00181 void
00182 BB_lp::logical_fixing(const BCP_lp_result& lpres,
00183                const BCP_vec<BCP_var*>& vars,
00184                const BCP_vec<BCP_cut*>& cuts,
00185                const BCP_vec<BCP_obj_status>& var_status,
00186                const BCP_vec<BCP_obj_status>& cut_status,
00187                const int var_bound_changes_since_logical_fixing,
00188                BCP_vec<int>& changed_pos, BCP_vec<double>& new_bd) 
00189 
00190   // Called at each iteration, after test_feasibility, if the node is not
00191   // fathomable
00192 {
00193 }
00194 
00195 /************************************************************************/
00196 void
00197 BB_lp::generate_cuts_in_lp(const BCP_lp_result& lpres,
00198                            const BCP_vec<BCP_var*>& vars,
00199                            const BCP_vec<BCP_cut*>& cuts,
00200                            BCP_vec<BCP_cut*>& new_cuts,
00201                            BCP_vec<BCP_row*>& new_rows)
00202 
00203   // Send to BCP the cuts generated in test_feasibility.
00204   // Use this function to generate standard cuts (Knapsack covers,
00205   // Lift-and-Project, odd holes, ...).
00206 
00207 {
00208    int i;
00209    
00210    for (i=violated_cuts.size()-1; i>=0; --i) {
00211       const int ind = violated_cuts[i];
00212       new_cuts.push_back(new BB_indexed_cut(ind, p_desc->rlb_indexed[ind], 
00213                                             p_desc->rub_indexed[ind]));
00214    }
00215    cout << "generate_cuts_in_lp(): found " << new_cuts.size() 
00216         << " indexed cuts" << endl;
00217 
00218    for(i=algo_cuts.size()-1; i>=0; --i) {
00219      new_cuts.push_back(algo_cuts[i]);
00220    }
00221    cout << "generate_cuts_in_lp(): found " << algo_cuts.size() 
00222         << " algorithmic cuts" << endl;
00223 
00224    algo_cuts.clear();
00225 }
00226 
00227 /************************************************************************/
00228 BCP_solution*
00229 BB_lp::generate_heuristic_solution(const BCP_lp_result& lpres,
00230                             const BCP_vec<BCP_var*>& vars,
00231                             const BCP_vec<BCP_cut*>& cuts) 
00232 {
00233 #ifdef HEUR_SOL
00234 
00235   // Simple rounding heuristic 
00236 
00237   int i, j, k, cn = p_desc->colnum;
00238   const double *x = lpres.x();
00239   double *rounded_x = new double[cn];
00240 
00241   for(i=0; i<cn; i++) {
00242     if(x[i] > 0.5 + EPS) {
00243         rounded_x[i] = 1;
00244     }
00245     else {
00246       rounded_x[i] = 0;
00247     }
00248   }
00249 
00250   // check if rounded_x is feasible or not
00251 
00252   int nb_rows_core = p_desc->core->getNumRows();
00253   int nb_rows_indexed = p_desc->indexed->getNumRows();
00254   int max_core_indexed;
00255 
00256   if(nb_rows_indexed > nb_rows_core) {
00257     max_core_indexed = nb_rows_indexed;
00258   }
00259   else {
00260     max_core_indexed = nb_rows_core;
00261   }
00262   
00263   double *check_lhs = new double[max_core_indexed];
00264 
00265   // check core constraints
00266 
00267   p_desc->core->times(rounded_x, check_lhs);
00268 
00269   for (i=0; i<nb_rows_core; ++i) {
00270     if ((check_lhs[i] < p_desc->rlb_core[i] - EPS) || 
00271         (check_lhs[i] > p_desc->rub_core[i] + EPS)) {
00272 
00273       delete[] check_lhs;
00274       delete[] rounded_x;
00275       cout << "generate_heuristic_solution() returns nothing.\n" 
00276            << endl;
00277       return(0);
00278     }
00279   }
00280 
00281   // check indexed constraints
00282 
00283   p_desc->indexed->times(rounded_x, check_lhs);
00284   for (i=0; i<nb_rows_indexed; ++i) {
00285     if ((check_lhs[i] < p_desc->rlb_indexed[i] - EPS) || 
00286         (check_lhs[i] > p_desc->rub_indexed[i] + EPS)) {
00287 
00288       delete[] check_lhs;
00289       delete[] rounded_x;
00290       cout << "generate_heuristic_solution() returns nothing.\n" 
00291            << endl;
00292       return(0);
00293     }
00294   }
00295 
00296   delete[] check_lhs;
00297 
00298   // check algorithmic constraints
00299 
00300   for(i=0; i<cn; i++) {
00301     j = (i+1) % cn;
00302     k = (i+2) % cn;
00303 
00304     if(x[i] + x[j] + x[k] > 1 + EPS) {
00305 
00306       delete[] rounded_x;
00307       cout << "generate_heuristic_solution() returns nothing.\n" 
00308            << endl;
00309       return(0);
00310     }
00311   }
00312 
00313   // rounded_x is feasible
00314     
00315   BCP_solution_generic *heur_sol = new BCP_solution_generic();
00316   heur_sol->_delete_vars = false; // otherwise BCP will delete the variables
00317                                   // appearing in the solution 
00318 
00319   BCP_vec<BCP_var_core *> core_vars = getLpProblemPointer()->core->vars;
00320   int ind;
00321 
00322   for(i=0; i<cn; i++) {
00323     ind = core_vars[i]->bcpind();
00324     if(rounded_x[ind] > EPS) {
00325       heur_sol->add_entry(core_vars[i], rounded_x[ind]);
00326     }
00327   }
00328   
00329   delete[] rounded_x;
00330   cout << "generate_heuristic_solution() returns a solution.\n" 
00331        << endl;
00332   return(heur_sol);
00333 
00334 #else /* not HEUR_SOL */
00335 
00336   // return nothing
00337 
00338   return(0);
00339 #endif
00340 }
00341 
00342 /************************************************************************/
00343 void
00344 BB_lp::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
00345                     BCP_vec<BCP_cut*>& cuts,       // what to expand
00346                     BCP_vec<BCP_row*>& rows,       // the expanded rows
00347                     // things that the user can use for lifting cuts if allowed
00348                     const BCP_lp_result& lpres,
00349                     BCP_object_origin origin, bool allow_multiple)
00350 
00351   // Required function when indexed or algorithmic cuts are used.
00352   // Describes how to get a row of the matrix from the representation of the
00353   // cut.
00354 
00355 {
00356   const int cutnum = cuts.size();
00357   for (int i=0; i<cutnum; ++i) {
00358       const BB_indexed_cut *icut =
00359         dynamic_cast<const BB_indexed_cut*>(cuts[i]);
00360       if (icut) {
00361         const int ind = icut->index();
00362         rows.push_back(new BCP_row(p_desc->indexed->getVector(ind),
00363                                    p_desc->rlb_indexed[ind], 
00364                                    p_desc->rub_indexed[ind]));
00365         continue;
00366       }
00367       
00368       const OsiRowCut* bcut = dynamic_cast<const BB_cut*>(cuts[i]);
00369       if (bcut) {
00370         rows.push_back(new BCP_row(bcut->row(), bcut->lb(), bcut->ub()));
00371         continue;
00372       }
00373       
00374       throw BCP_fatal_error("Unknown cut type in cuts_to_rows.\n");
00375   }
00376 }
00377 
00378 /********************************************************************/
00379 BCP_branching_decision
00380 BB_lp::select_branching_candidates(const BCP_lp_result& lpres,
00381                                    const BCP_vec<BCP_var*>& vars,
00382                                    const BCP_vec<BCP_cut*>& cuts,
00383                                    const BCP_lp_var_pool& local_var_pool,
00384                                    const BCP_lp_cut_pool& local_cut_pool,
00385                                    BCP_vec<BCP_lp_branching_object*>& cands,
00386                                    bool force_branch)
00387 {
00388 #ifdef CUSTOM_BRANCH
00389 
00390   // Called at the end of each iteration. Possible return values are:
00391 
00392   // BCP_DoNotBranch_Fathomed : The node should be fathomed without branching
00393 
00394   // BCP_DoNotBranch :          BCP should continue to work on this node
00395 
00396   // BCP_DoBranch :             Branching must be done. In this case the 
00397   //                            method returns the branching object candidates 
00398   //                            in one of the arguments
00399 
00400   // Don't branch if cutting planes have been generated
00401 
00402   if(local_cut_pool.size() > 0) {
00403     
00404     cout << "select_branching_candidates() returns BCP_DoNotBranch" 
00405          << endl;
00406     
00407     return(BCP_DoNotBranch);
00408   }
00409   else {
00410    
00411     // Branch on the first fractional variable
00412 
00413     int i;
00414     const double *x = lpres.x();
00415     BCP_vec<int> select_pos;
00416     
00417     for(i=0; i<p_desc->colnum; i++) {
00418       
00419       if((x[i] > EPS) && (x[i] < 1-EPS)) {
00420         select_pos.push_back(i);
00421 
00422         // put in cands all variables whose index are in select_pos
00423 
00424         append_branching_vars(lpres.x(), vars, select_pos, cands);
00425 
00426         cout << "Branching on variable: " << i << " (" << x[i]
00427              << ")   depth: " << current_level() << endl;
00428         break;
00429       }
00430     }
00431   
00432     if (cands.size() == 0) {
00433       throw BCP_fatal_error("select_banching_candidate() : No cut in pool but couldn't select branching variable.\n");
00434     }
00435 
00436 
00437     cout << "select_branching_candidates() returns BCP_DoBranch" << endl;
00438 
00439     return BCP_DoBranch;
00440   }
00441 
00442 #else 
00443   return(BCP_lp_user::select_branching_candidates(lpres, vars, cuts, 
00444                                                   local_var_pool, 
00445                                                   local_cut_pool, cands,
00446                                                   force_branch));
00447 #endif
00448 }
00449 
00450 /**************************************************************************/
00451 void
00452 BB_lp::set_user_data_for_children(BCP_presolved_lp_brobj* best, 
00453                                   const int selected)
00454 
00455   // Given the branching decision (parameter "best"; "selected" is the
00456   // index of the chosen branching decision in the candidate list), 
00457   // set the user data for the children.
00458 
00459 {
00460 #ifdef USER_DATA
00461   BCP_lp_branching_object *cand = best->candidate();
00462   MY_user_data *curr_ud = dynamic_cast<MY_user_data*> (get_user_data());
00463   real_user_data *curr_rud = curr_ud->p_rud;
00464 
00465   for(int i=0; i<cand->child_num; i++) {
00466     MY_user_data *ud = new MY_user_data(curr_rud->max_card_set_zero);
00467     real_user_data *rud = ud->p_rud;
00468 
00469     rud->card_set_zero = curr_rud->card_set_zero;
00470 
00471     for(int j=0; j<curr_rud->card_set_zero; j++) {
00472       rud->set_zero[j] = curr_rud->set_zero[j];
00473     }
00474 
00475     int ind_br = (*(cand->forced_var_pos))[0];
00476 
00477     if((*(cand->forced_var_bd))[2*i + 1] < EPS) {
00478       rud->set_zero[curr_rud->card_set_zero] = ind_br;
00479       (rud->card_set_zero)++;
00480     }
00481     best->user_data()[i] = ud;
00482   }
00483 #endif /* USER_DATA */
00484 } /* set_user_data_for_children */
00485 
00486 
00487 
00488 
00489 

Generated on Wed Nov 30 03:03:46 2011 by  doxygen 1.4.7