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

Generated on Mon Aug 3 03:02:14 2009 by  doxygen 1.4.7