/home/coin/SVN-release/OS-2.4.2/Bcp/src/LP/BCP_lp_user.cpp

Go to the documentation of this file.
00001 // Copyright (C) 2000, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 #include <cmath>
00004 #include <cstdarg>
00005 #include <cstdio>
00006 
00007 #include "CoinHelperFunctions.hpp"
00008 #include "CoinTime.hpp"
00009 
00010 #include "BCP_math.hpp"
00011 #include "BCP_error.hpp"
00012 #include "BCP_USER.hpp"
00013 #include "BCP_var.hpp"
00014 #include "BCP_cut.hpp"
00015 #include "BCP_problem_core.hpp"
00016 #include "BCP_lp_user.hpp"
00017 #include "BCP_lp.hpp"
00018 #include "BCP_lp_pool.hpp"
00019 #include "BCP_lp_result.hpp"
00020 #include "BCP_lp_node.hpp"
00021 #include "BCP_lp_branch.hpp"
00022 #include "BCP_problem_core.hpp"
00023 #include "BCP_solution.hpp"
00024 #include "BCP_functions.hpp"
00025 #include "BCP_lp_functions.hpp"
00026 
00027 //#############################################################################
00028 // Informational methods for the user
00029 double BCP_lp_user::upper_bound() const    { return p->ub(); }
00030 bool BCP_lp_user::over_ub(double lb) const { return p->over_ub(lb); }
00031 int BCP_lp_user::current_phase() const     { return p->phase; }
00032 int BCP_lp_user::current_level() const     { return p->node->level; }
00033 int BCP_lp_user::current_index() const     { return p->node->index; }
00034 int BCP_lp_user::current_iteration() const { return p->node->iteration_count; }
00035 double BCP_lp_user::start_time() const     { return p->start_time; }
00036 BCP_user_data* BCP_lp_user::get_user_data() { return p->node->user_data; }
00037 
00038 void BCP_lp_user::print(const bool ifprint, const char * format, ...) const {
00039   if (ifprint) {
00040     va_list valist;
00041     va_start(valist, format);
00042     vprintf(format, valist);
00043     va_end(valist);
00044   }
00045 }
00046 
00047 //#############################################################################
00048 // Informational methods for the user
00049 /* Methods to get/set BCP parameters on the fly */
00050 char
00051 BCP_lp_user::get_param(const BCP_lp_par::chr_params key) const
00052 { return p->par.entry(key); }
00053 int
00054 BCP_lp_user::get_param(const BCP_lp_par::int_params key) const
00055 { return p->par.entry(key); }
00056 double
00057 BCP_lp_user::get_param(const BCP_lp_par::dbl_params key) const
00058 { return p->par.entry(key); }
00059 const BCP_string&
00060 BCP_lp_user::get_param(const BCP_lp_par::str_params key) const
00061 { return p->par.entry(key); }
00062 
00063 void BCP_lp_user::set_param(const BCP_lp_par::chr_params key, const bool val)
00064 { p->par.set_entry(key, val); }
00065 void BCP_lp_user::set_param(const BCP_lp_par::chr_params key, const char val)
00066 { p->par.set_entry(key, val); }
00067 void BCP_lp_user::set_param(const BCP_lp_par::int_params key, const int val)
00068 { p->par.set_entry(key, val); }
00069 void BCP_lp_user::set_param(const BCP_lp_par::dbl_params key, const double val)
00070 { p->par.set_entry(key, val); }
00071 void BCP_lp_user::set_param(const BCP_lp_par::str_params key, const char * val)
00072 { p->par.set_entry(key, val); }
00073 
00074 //#############################################################################
00075 
00076 void
00077 BCP_lp_user::send_feasible_solution(const BCP_solution* sol)
00078 {
00079     // send the feas sol to the tree manager
00080     p->msg_buf.clear();
00081     pack_feasible_solution(p->msg_buf, sol);
00082     p->msg_env->send(p->get_parent() /*tree_manager*/,
00083                      BCP_Msg_FeasibleSolution, p->msg_buf);
00084 
00085     // update the UB if necessary
00086     const double obj = sol->objective_value();
00087     if (p->ub(obj) && p->node->colgen != BCP_GenerateColumns) {
00088         // FIXME: If we had a flag in the node that indicates not to
00089         // generate cols in it and in its descendants then the dual obj
00090         // limit could still be set...
00091         p->lp_solver->setDblParam(OsiDualObjectiveLimit, obj-p->granularity());
00092     }
00093 }
00094 
00095 //#############################################################################
00096 
00097 // These functions are member functions of the VIRTUAL class BCP_lp_user.
00098 // Therefore any concrete class derived from BCP_lp_user can override these
00099 // definitions. These are here as default functions.
00100 
00101 //#############################################################################
00102 // a few helper functions for selecting a subset of a double vector
00103 void
00104 BCP_lp_user::select_nonzeros(const double * first, const double * last,
00105                              const double etol,
00106                              BCP_vec<int>& nonzeros) const {
00107     nonzeros.reserve(last - first);
00108     BCP_vec<double>::const_iterator current = first;
00109     for ( ; current != last; ++current)
00110         if (CoinAbs(*current) > etol)
00111             nonzeros.unchecked_push_back(current - first);
00112 }
00113 //-----------------------------------------------------------------------------
00114 void
00115 BCP_lp_user::select_zeros(const double * first, const double * last,
00116                           const double etol,
00117                           BCP_vec<int>& zeros) const {
00118     zeros.reserve(last - first);
00119     BCP_vec<double>::const_iterator current = first;
00120     for ( ; current != last; ++current)
00121         if (CoinAbs(*current) <= etol)
00122             zeros.unchecked_push_back(current - first);
00123 }
00124 //-----------------------------------------------------------------------------
00125 void
00126 BCP_lp_user::select_positives(const double * first, const double * last,
00127                               const double etol,
00128                               BCP_vec<int>& positives) const {
00129     positives.reserve(last - first);
00130     BCP_vec<double>::const_iterator current = first;
00131     for ( ; current != last; ++current)
00132         if (*current > etol)
00133             positives.unchecked_push_back(current - first);
00134 }
00135 //-----------------------------------------------------------------------------
00136 void
00137 BCP_lp_user::select_fractions(const double * first, const double * last,
00138                               const double etol,
00139                               BCP_vec<int>& fractions) const {
00140     fractions.reserve(last - first);
00141     BCP_vec<double>::const_iterator current = first;
00142     for ( ; current != last; ++current)
00143         if (*current-floor(*current) > etol && ceil(*current)-*current > etol)
00144             fractions.unchecked_push_back(current - first);
00145 }
00146 
00147 //#############################################################################
00148 //#############################################################################
00149 // unpack the initial info for the appropriate process
00150 void
00151 BCP_lp_user::unpack_module_data(BCP_buffer & buf)
00152 {
00153   print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00154         "LP: Default unpack_module_data() executed.\n");
00155 }
00156 
00157 //#############################################################################
00158 
00160 int
00161 BCP_lp_user::process_id() const
00162 {
00163     return p->get_process_id();
00164 }
00165 
00167 int 
00168 BCP_lp_user::parent() const
00169 {
00170     return p->get_parent();
00171 }
00172 
00174 void
00175 BCP_lp_user::send_message(const int target, const BCP_buffer& buf,
00176                           BCP_message_tag tag)
00177 {
00178     p->msg_env->send(target, tag, buf);
00179 }    
00180 
00182 void
00183 BCP_lp_user::receive_message(const int sender, BCP_buffer& buf,
00184                              BCP_message_tag tag)
00185 {
00186     p->msg_env->receive(sender, tag, buf, -1);
00187 }
00188 
00190 void
00191 BCP_lp_user::broadcast_message(const BCP_process_t proc_type,
00192                                const BCP_buffer& buf)
00193 {
00194     throw BCP_fatal_error("\
00195 BCP_lp_user::broadcast_message: can't broadcast from an LP process.\n");
00196 }
00197 
00200 void
00201 BCP_lp_user::process_message(BCP_buffer& buf)
00202 {
00203     throw BCP_fatal_error("\
00204 BCP_lp_user::process_message() invoked but not overridden!\n");
00205 }
00206 
00207 //#############################################################################
00208 
00209 OsiSolverInterface * 
00210 BCP_lp_user::initialize_solver_interface()
00211 {
00212     throw BCP_fatal_error("\
00213 BCP_lp_user::initialize_solver_interface() invoked but not overridden!\n");
00214     return 0;
00215 }
00216 
00217 //#############################################################################
00218 
00219 void 
00220 BCP_lp_user::initialize_int_and_sos_list(std::vector<OsiObject *>& intAndSosObjects)
00221 {
00222     return;
00223 }
00224 
00225 //#############################################################################
00226 void
00227 BCP_lp_user::initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
00228                                              const BCP_vec<BCP_cut*>& cuts,
00229                                              const BCP_vec<BCP_obj_status>& vs,
00230                                              const BCP_vec<BCP_obj_status>& cs,
00231                                              BCP_vec<int>& var_changed_pos,
00232                                              BCP_vec<double>& var_new_bd,
00233                                              BCP_vec<int>& cut_changed_pos,
00234                                              BCP_vec<double>& cut_new_bd)
00235 {
00236   print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00237         "LP: Default initialize_new_search_tree_node() executed.\n");
00238 }
00239 
00240 //#############################################################################
00241 
00242 void
00243 BCP_lp_user::load_problem(OsiSolverInterface& osi, BCP_problem_core* core,
00244                           BCP_var_set& vars, BCP_cut_set& cuts)
00245 {
00246   const int varnum = vars.size();
00247   const int cutnum = cuts.size();
00248 
00249   if (varnum == 0) {
00250     // *FIXME* : kill all the processes
00251     throw BCP_fatal_error("There are no vars in the description for node %i!\n",
00252                           current_index());
00253   }
00254 
00255   if (cutnum == 0) {
00256     // *FIXME* : kill all the processes
00257     throw BCP_fatal_error("There are no cuts in the description for node %i!\n",
00258                           current_index());
00259   }
00260 
00261   BCP_vec<BCP_col*> cols;
00262   BCP_vec<BCP_row*> rows;
00263 
00264   int bvarnum = core->varnum();
00265   int bcutnum = core->cutnum();
00266 
00267   if (varnum == 0 || cutnum == 0){
00268     throw BCP_fatal_error("No rows or no cols to create a matrix from!\n");
00269   }
00270 
00271   BCP_lp_relax* m = 0;
00272   if (bvarnum == 0) {
00273     // no core vars. doesn't matter if there're any core cuts, the starting
00274     // matrix is computed the same way
00275     cols.reserve(varnum);
00276     vars_to_cols(cuts, vars, cols,
00277                  *p->lp_result, BCP_Object_FromTreeManager, false);
00278     BCP_vec<double> RLB;
00279     BCP_vec<double> RUB;
00280     RLB.reserve(cutnum);
00281     RUB.reserve(cutnum);
00282     BCP_cut_set::const_iterator ci = cuts.begin();
00283     BCP_cut_set::const_iterator lastci = cuts.end();
00284     for ( ; ci != lastci; ++ci) {
00285       RLB.unchecked_push_back((*ci)->lb());
00286       RUB.unchecked_push_back((*ci)->ub());
00287     }
00288     m = new BCP_lp_relax(cols, RLB, RUB);
00289     purge_ptr_vector(cols);
00290   } else {
00291     if (bcutnum == 0) {
00292       rows.reserve(cutnum);
00293       cuts_to_rows(vars, cuts, rows,
00294                    *p->lp_result, BCP_Object_FromTreeManager, false);
00295       BCP_vec<double> CLB;
00296       BCP_vec<double> CUB;
00297       BCP_vec<double> OBJ;
00298       CLB.reserve(varnum);
00299       CUB.reserve(varnum);
00300       OBJ.reserve(varnum);
00301       BCP_var_set::const_iterator vi = vars.begin();
00302       BCP_var_set::const_iterator lastvi = vars.end();
00303       for ( ; vi != lastvi; ++vi) {
00304         CLB.unchecked_push_back((*vi)->lb());
00305         CUB.unchecked_push_back((*vi)->ub());
00306         OBJ.unchecked_push_back((*vi)->obj());
00307       }
00308       m = new BCP_lp_relax(rows, CLB, CUB, OBJ);
00309       purge_ptr_vector(rows);
00310     } else {
00311       // has core vars and cuts, the starting matrix is the core matrix
00312       m = core->matrix;
00313     }
00314   }
00315    
00316   // We have the description in p.node. Load it into the lp solver.
00317   // First load the core matrix
00318   osi.loadProblem(*m,
00319                   m->ColLowerBound().begin(), m->ColUpperBound().begin(),
00320                   m->Objective().begin(),
00321                   m->RowLowerBound().begin(), m->RowUpperBound().begin());
00322 
00323   if (bvarnum > 0 && bcutnum > 0) {
00324     //-----------------------------------------------------------------------
00325     // We have to add the 'added' stuff only if we had a core matrix
00326     //-----------------------------------------------------------------------
00327     // Add the Named and Algo cols if there are any (and if we have any cols)
00328     if (varnum > bvarnum && bcutnum > 0) {
00329       BCP_vec<BCP_var*> added_vars(vars.entry(bvarnum), vars.end());
00330       BCP_vec<BCP_cut*> core_cuts(cuts.begin(), cuts.entry(bcutnum));
00331       cols.reserve(vars.size());
00332       vars_to_cols(core_cuts, added_vars, cols,
00333                    *p->lp_result, BCP_Object_FromTreeManager, false);
00334       BCP_lp_add_cols_to_lp(cols, &osi);
00335       purge_ptr_vector(cols);
00336     }
00337     //-----------------------------------------------------------------------
00338     // Add the Named and Algo rows if there are any (and if we have any such
00339     // rows AND if there are core vars)
00340     if (cutnum > bcutnum) {
00341       BCP_vec<BCP_cut*> added_cuts(cuts.entry(bcutnum), cuts.end());
00342       rows.reserve(added_cuts.size());
00343       BCP_fatal_error::abort_on_error = false;
00344       try {
00345          cuts_to_rows(vars, added_cuts, rows,
00346                       *p->lp_result, BCP_Object_FromTreeManager, false);
00347       }
00348       catch (...) {
00349       }
00350       BCP_fatal_error::abort_on_error = true;
00351       BCP_lp_add_rows_to_lp(rows, &osi);
00352       purge_ptr_vector(rows);
00353     }
00354   } else {
00355     // Otherwise (i.e., if we had no core matrix) we just have to get rid of
00356     // 'm'.
00357     delete m;
00358   }
00359   const bool dumpcuts = get_param(BCP_lp_par::Lp_DumpNodeDescCuts);
00360   const bool dumpvars = get_param(BCP_lp_par::Lp_DumpNodeDescVars);
00361   if (dumpcuts || dumpvars) {
00362     printf("LP: NEW ACTIVE NODE %i DUMP START ===========================\n",
00363            current_index());
00364     if (dumpvars) {
00365       printf("    Core Vars (bvarnum %i)\n", bvarnum);
00366       for (int i = 0; i < bvarnum; ++i) {
00367         printf("        var %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
00368                i, vars[i]->bcpind(), vars[i]->status(),
00369                vars[i]->lb(), vars[i]->ub());
00370       }
00371     }
00372     if (dumpcuts) {
00373       printf("    Core Cuts (bcutnum %i)\n", bcutnum);
00374       for (int i = 0; i < bcutnum; ++i) {
00375         printf("        cut %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
00376                i, cuts[i]->bcpind(), cuts[i]->status(),
00377                cuts[i]->lb(), cuts[i]->ub());
00378       }
00379     }
00380     if (dumpvars) {
00381       printf("    Extra Vars (extra varnum %i)\n", varnum - bvarnum);
00382       for (int i = bvarnum; i < varnum; ++i) {
00383         printf("        var %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
00384                i, vars[i]->bcpind(), vars[i]->status(),
00385                vars[i]->lb(), vars[i]->ub());
00386       }
00387     }
00388     if (dumpcuts) {
00389       printf("    Extra Cuts (extra cutnum %i)\n", cutnum - bcutnum);
00390       for (int i = bcutnum; i < cutnum; ++i) {
00391         printf("        cut %4i: bcpind %i,   status %i,   lb %f,   ub %f\n",
00392                i, cuts[i]->bcpind(), cuts[i]->status(),
00393                cuts[i]->lb(), cuts[i]->ub());
00394       }
00395     }
00396     printf("LP: NEW ACTIVE NODE %i DUMP END ===========================\n",
00397            current_index());
00398   }
00399 }
00400 
00401 //#############################################################################
00402 // Opportunity to reset things before optimization
00403 void
00404 BCP_lp_user::modify_lp_parameters(OsiSolverInterface* lp, const int changeType,
00405                                   bool in_strong_branching)
00406 {
00407   print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00408         "LP: Default prepare_for_optimization() executed.\n");
00409 
00410   if (changeType == 1) { // primal feas was affected
00411     lp->setHintParam(OsiDoDualInResolve, true, OsiHintTry);
00412     return;
00413   }
00414   if (changeType == 2) { // dual feas was affected
00415     lp->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
00416     return;
00417   }
00418   // Neither or both, so let the solver decide
00419   lp->setHintParam(OsiDoDualInResolve, true, OsiHintIgnore);
00420 }
00421 
00422 //#############################################################################
00423 void
00424 BCP_lp_user::process_lp_result(const BCP_lp_result& lpres,
00425                                const BCP_vec<BCP_var*>& vars,
00426                                const BCP_vec<BCP_cut*>& cuts,
00427                                const double old_lower_bound,
00428                                double& true_lower_bound,
00429                                BCP_solution*& sol,
00430                                BCP_vec<BCP_cut*>& new_cuts,
00431                                BCP_vec<BCP_row*>& new_rows,
00432                                BCP_vec<BCP_var*>& new_vars,
00433                                BCP_vec<BCP_col*>& new_cols)
00434 {
00435     p->user_has_lp_result_processing = false;
00436 }
00437 
00438 //#############################################################################
00439 // Generating a true lower bound
00440 double
00441 BCP_lp_user::compute_lower_bound(const double old_lower_bound,
00442                                  const BCP_lp_result& lpres,
00443                                  const BCP_vec<BCP_var*>& vars,
00444                                  const BCP_vec<BCP_cut*>& cuts)
00445 {
00446     // If columns are to be generated then we can't say anything, just return
00447     // the current lower bound
00448     if (p->node->colgen != BCP_DoNotGenerateColumns_Fathom)
00449         return old_lower_bound;
00450 
00451     // Otherwise we got the examine the termination code and the objective
00452     // value of the LP solution
00453     const int tc = lpres.termcode();
00454     if (tc & BCP_ProvenOptimal)
00455         return lpres.objval();
00456 
00457     // The limit (the upper bound) on the dual objective is proven to be
00458     // reached, but the objval might not reflect this! (the LP solver may not
00459     // make the last iteration that pushes objval over the limit). So we return
00460     // a high value ourselves.
00461     if (tc & BCP_DualObjLimReached)
00462         return p->ub() + 1e-5;
00463 
00464     // We can't say anything in any other case
00465     // (BCP_ProvenPrimalInf | BCP_ProvenDualInf | BCP_PrimalObjLimReached |
00466     //  BCP_TimeLimit | BCP_Abandoned), not to mention that some of these are
00467     //  impossible. Just return the current bound.
00468     return old_lower_bound;
00469 }
00470 
00471 //#############################################################################
00472 // Feasibility testing
00473 BCP_solution*
00474 BCP_lp_user::test_feasibility(const BCP_lp_result& lpres,
00475                               const BCP_vec<BCP_var*>& vars,
00476                               const BCP_vec<BCP_cut*>& cuts)
00477 {
00478     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00479           "LP: Default test_feasibility() executed.\n");
00480 
00481     const double etol = p->param(BCP_lp_par::IntegerTolerance);
00482     BCP_feasibility_test test =
00483         static_cast<BCP_feasibility_test>(p->param(BCP_lp_par::FeasibilityTest));
00484 
00485     BCP_solution* sol = NULL;
00486     switch (test){
00487     case BCP_Binary_Feasible:   sol = test_binary(lpres, vars, etol); break;
00488     case BCP_Integral_Feasible: sol = test_integral(lpres, vars, etol); break;
00489     case BCP_FullTest_Feasible: sol = test_full(lpres, vars, etol); break;
00490     default:
00491         throw BCP_fatal_error("LP: unknown test_feasibility() rule.\n");
00492     }
00493 
00494     return sol;
00495 }
00496 
00497 //-----------------------------------------------------------------------------
00498 
00499 BCP_solution_generic*
00500 BCP_lp_user::test_binary(const BCP_lp_result& lpres,
00501                          const BCP_vec<BCP_var*>& vars,
00502                          const double etol) const
00503 {
00504     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00505           "LP: Default test_binary() executed.\n");
00506 
00507     // Do anything only if the termination code is sensible
00508     const int tc = lpres.termcode();
00509     if (! (tc & BCP_ProvenOptimal))
00510         return 0;
00511 
00512     const double * x = lpres.x();
00513     const int varnum = vars.size();
00514     const double etol1 = 1 - etol;
00515     int i;
00516 
00517     for (i = 0 ; i < varnum; ++i) {
00518         const double xi = x[i];
00519         if (xi > etol && xi < etol1)
00520             return(NULL);
00521     }
00522 
00523     // This solution does not own the pointers to the variables
00524     BCP_solution_generic* sol = new BCP_solution_generic(false); 
00525     double obj = 0;
00526     for (i = 0 ; i < varnum; ++i) {
00527         if (x[i] > etol1) {
00528             sol->add_entry(vars[i], 1.0);
00529             obj += vars[i]->obj();
00530         }
00531     }
00532     sol->_objective = obj;
00533     return sol;
00534 }
00535 //-----------------------------------------------------------------------------
00536 BCP_solution_generic*
00537 BCP_lp_user::test_integral(const BCP_lp_result& lpres,
00538                            const BCP_vec<BCP_var*>& vars,
00539                            const double etol) const
00540 {
00541     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00542           "LP: Default test_integral() executed.\n");
00543 
00544     // Do anything only if the termination code is sensible
00545     const int tc = lpres.termcode();
00546     if (! (tc & BCP_ProvenOptimal))
00547         return 0;
00548 
00549     const double * x = lpres.x();
00550     const int varnum = vars.size();
00551     const double etol1 = 1 - etol;
00552     int i;
00553 
00554     for (i = 0 ; i < varnum; ++i) {
00555         const double val = x[i] - floor(x[i]);
00556         if (val > etol && val < etol1)
00557             return(NULL);
00558     }
00559   
00560     // This solution does not own the pointers to the variables
00561     BCP_solution_generic* sol = new BCP_solution_generic(false);
00562     double obj = 0;
00563     for (i = 0 ; i < varnum; ++i) {
00564         const double val = floor(x[i] + 0.5);
00565         if (val != 0.0) { // it's OK to test the double against 0,
00566                           // since we took the floor()
00567             sol->add_entry(vars[i], val);
00568             obj += val * vars[i]->obj();
00569         }
00570     }
00571     sol->_objective = obj;
00572     return sol;
00573 }
00574 //-----------------------------------------------------------------------------
00575 BCP_solution_generic*
00576 BCP_lp_user::test_full(const BCP_lp_result& lpres,
00577                        const BCP_vec<BCP_var*>& vars,
00578                        const double etol) const
00579 {
00580     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00581           "LP: Default test_full() executed.\n");
00582 
00583     // Do anything only if the termination code is sensible
00584     const int tc = lpres.termcode();
00585     if (! (tc & BCP_ProvenOptimal))
00586         return 0;
00587 
00588     const double * x = lpres.x();
00589     const int varnum = vars.size();
00590     const double etol1 = 1 - etol;
00591     int i;
00592 
00593     for (i = 0 ; i < varnum; ++i) {
00594         switch (vars[i]->var_type()){
00595         case BCP_BinaryVar:
00596             {
00597                 const double val = x[i];
00598                 if (val > etol && val < etol1)
00599                     return(NULL);
00600             }
00601             break;
00602         case BCP_IntegerVar:
00603             {
00604                 const double val = x[i] - floor(x[i]);
00605                 if (val > etol && val < etol1)
00606                     return(NULL);
00607             }
00608             break;
00609         case BCP_ContinuousVar:
00610             break;
00611         }
00612     }
00613 
00614     // This solution does not own the pointers to the variables
00615     BCP_solution_generic* sol = new BCP_solution_generic(false);
00616     double obj = 0;
00617     for (i = 0 ; i < varnum; ++i) {
00618         // round the integer vars
00619         const double val = (vars[i]->var_type() == BCP_ContinuousVar) ?
00620             x[i] : floor(x[i] + 0.5);
00621         if (val < -etol || val > etol) {
00622             sol->add_entry(vars[i], val);
00623             obj += val * vars[i]->obj();
00624         }
00625     }
00626     sol->_objective = obj;
00627     return sol;
00628 }
00629 
00630 //#############################################################################
00631 // Heuristic solution generation
00632 BCP_solution*
00633 BCP_lp_user::generate_heuristic_solution(const BCP_lp_result& lpres,
00634                                          const BCP_vec<BCP_var*>& vars,
00635                                          const BCP_vec<BCP_cut*>& cuts)
00636 {
00637     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00638           "LP: Default generate_heuristic_solution() executed.\n");
00639     return NULL;
00640 }
00641 //#############################################################################
00642 // Feasible solution packing
00643 void
00644 BCP_lp_user::pack_feasible_solution(BCP_buffer& buf, const BCP_solution* sol)
00645 {
00646     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00647           "LP: Default pack_feasible_solution() executed.\n");
00648 
00649     const BCP_solution_generic* gensol =
00650         dynamic_cast<const BCP_solution_generic*>(sol);
00651     const BCP_vec<BCP_var*> solvars = gensol->_vars;
00652     const BCP_vec<double> values = gensol->_values;
00653     const int size = solvars.size();
00654     buf.pack(size);
00655     for (int i = 0; i < size; ++i) {
00656         buf.pack(values[i]);
00657         p->pack_var(*solvars[i]);
00658     }
00659     buf.pack(gensol->objective_value());
00660 }
00661 
00662 //#############################################################################
00663 // pack message for CG/CP
00664 
00665 void
00666 BCP_lp_user::pack_primal_solution(BCP_buffer& buf,
00667                                   const BCP_lp_result& lpres,
00668                                   const BCP_vec<BCP_var*>& vars,
00669                                   const BCP_vec<BCP_cut*>& cuts)
00670 {
00671     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00672           "LP: Default pack_for_cg() executed.\n");
00673 
00674     BCP_vec<int> coll;
00675 
00676     const double petol = lpres.primalTolerance();
00677     const double * x = lpres.x();
00678     const int varnum = vars.size();
00679 
00680     switch (p->param(BCP_lp_par::InfoForCG)){
00681     case BCP_PrimalSolution_Nonzeros:
00682         select_nonzeros(x, x + varnum, petol, coll);
00683         buf.set_msgtag(BCP_Msg_ForCG_PrimalNonzeros);
00684         break;
00685     case BCP_PrimalSolution_Fractions:
00686         select_fractions(x, x+varnum, petol, coll);
00687         buf.set_msgtag(BCP_Msg_ForCG_PrimalFractions);
00688         break;
00689     case BCP_PrimalSolution_Full:
00690         coll.reserve(varnum);
00691         for (int i = 0; i < varnum; ++i) {
00692             coll.unchecked_push_back(i);
00693         }
00694         buf.set_msgtag(BCP_Msg_ForCG_PrimalFull);
00695         break;
00696     default:
00697         throw BCP_fatal_error("Incorrect msgtag in pack_for_cg() !\n");
00698     }
00699 
00700     const int size = coll.size();
00701     buf.pack(size);
00702     if (size > 0){
00703         BCP_vec<int>::const_iterator pos = coll.begin() - 1;
00704         const BCP_vec<int>::const_iterator last_pos = coll.end();
00705         while (++pos != last_pos) {
00706             buf.pack(x[*pos]);
00707             p->pack_var(*vars[*pos]);
00708         }
00709     }
00710 }
00711 
00712 //=============================================================================
00713 // pack message for VG/VP
00714 
00715 void
00716 BCP_lp_user::pack_dual_solution(BCP_buffer& buf,
00717                                 const BCP_lp_result& lpres,
00718                                 const BCP_vec<BCP_var*>& vars,
00719                                 const BCP_vec<BCP_cut*>& cuts)
00720 {
00721     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00722           "LP: Default pack_for_vg() executed.\n");
00723 
00724     BCP_vec<int> coll;
00725 
00726     const double detol = lpres.dualTolerance();
00727     const double * pi = lpres.pi();
00728     const int cutnum = cuts.size();
00729 
00730     switch (p->param(BCP_lp_par::InfoForVG)){
00731     case BCP_DualSolution_Nonzeros:
00732         select_nonzeros(pi, pi+cutnum, detol, coll);
00733         buf.set_msgtag(BCP_Msg_ForVG_DualNonzeros);
00734         break;
00735     case BCP_DualSolution_Full:
00736         coll.reserve(cutnum);
00737         for (int i = 0; i < cutnum; ++i) {
00738             coll.unchecked_push_back(i);
00739         }
00740         buf.set_msgtag(BCP_Msg_ForVG_DualFull);
00741         break;
00742     default:
00743         throw BCP_fatal_error("Incorrect msgtag in pack_lp_solution() !\n");
00744     }
00745 
00746     const int size = coll.size();
00747     buf.pack(size);
00748     if (size > 0){
00749         BCP_vec<int>::const_iterator pos = coll.begin() - 1;
00750         const BCP_vec<int>::const_iterator last_pos = coll.end();
00751         while (++pos != last_pos) {
00752             buf.pack(pi[*pos]);
00753             p->pack_cut(*cuts[*pos]);
00754         }
00755     }
00756 }
00757 
00758 //#############################################################################
00759 // lp solution displaying
00760 void
00761 BCP_lp_user::display_lp_solution(const BCP_lp_result& lpres,
00762                                  const BCP_vec<BCP_var*>& vars,
00763                                  const BCP_vec<BCP_cut*>& cuts,
00764                                  const bool final_lp_solution)
00765 {
00766     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00767           "LP: Default display_lp_solution() executed.\n");
00768 
00769     if (final_lp_solution) {
00770         if (! p->param(BCP_lp_par::LpVerb_FinalRelaxedSolution))
00771             return;
00772         print(true,"  LP : Displaying LP solution (FinalRelaxedSolution) :\n");
00773     } else {
00774         if (! p->param(BCP_lp_par::LpVerb_RelaxedSolution))
00775             return;
00776         print(true,"  LP : Displaying LP solution (RelaxedSolution) :\n");
00777     }
00778 
00779     const double ietol = p->param(BCP_lp_par::IntegerTolerance);
00780 
00781     print(true, "  LP : Displaying solution :\n");
00782     BCP_vec<int> coll;
00783     const double * x = lpres.x();
00784     select_nonzeros(x, x+vars.size(), ietol, coll);
00785     const int size = coll.size();
00786     for (int i = 0; i < size; ++i) {
00787         const int ind = coll[i];
00788         vars[ind]->display(x[ind]);
00789     }
00790 }
00791 
00792 //#############################################################################
00793 // restoring feasibility
00794 void
00795 BCP_lp_user::restore_feasibility(const BCP_lp_result& lpres,
00796                                  const std::vector<double*> dual_rays,
00797                                  const BCP_vec<BCP_var*>& vars,
00798                                  const BCP_vec<BCP_cut*>& cuts,
00799                                  BCP_vec<BCP_var*>& vars_to_add,
00800                                  BCP_vec<BCP_col*>& cols_to_add)
00801 {
00802     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00803           "LP: Default restore_feasibility() executed.\n");
00804 }
00805 
00806 //#############################################################################
00807 //#############################################################################
00808 
00809 void
00810 BCP_lp_user::cuts_to_rows(const BCP_vec<BCP_var*>& vars, // on what to expand
00811                           BCP_vec<BCP_cut*>& cuts,       // what to expand
00812                           BCP_vec<BCP_row*>& rows,       // the expanded rows
00813                           // things that the user can use for lifting cuts
00814                           // if allowed
00815                           const BCP_lp_result& lpres,
00816                           BCP_object_origin origin, bool allow_multiple)
00817 {
00818     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00819           "LP: Default cuts_to_rows() executed.\n");
00820     throw BCP_fatal_error("cuts_to_rows() missing.\n");
00821 }
00822 //-----------------------------------------------------------------------------
00823 void
00824 BCP_lp_user::vars_to_cols(const BCP_vec<BCP_cut*>& cuts, // on what to expand
00825                           BCP_vec<BCP_var*>& vars,       // what to expand
00826                           BCP_vec<BCP_col*>& cols,       // the expanded cols
00827                           // things that the user can use for lifting vars
00828                           // if allowed
00829                           const BCP_lp_result& lpres,
00830                           BCP_object_origin origin, bool allow_multiple)
00831 {
00832     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00833         "LP: Default vars_to_cols() executed.\n");
00834     throw BCP_fatal_error("vars_to_cols() missing.\n");
00835 }
00836 
00837 //#############################################################################
00838 
00839 void
00840 BCP_lp_user::generate_cuts_in_lp(const BCP_lp_result& lpres,
00841                                  const BCP_vec<BCP_var*>& vars,
00842                                  const BCP_vec<BCP_cut*>& cuts,
00843                                  BCP_vec<BCP_cut*>& new_cuts,
00844                                  BCP_vec<BCP_row*>& new_rows)
00845 {
00846     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00847           "LP: Default generate_cuts_in_lp() executed.\n");
00848 }
00849 //-----------------------------------------------------------------------------
00850 void
00851 BCP_lp_user::generate_vars_in_lp(const BCP_lp_result& lpres,
00852                                  const BCP_vec<BCP_var*>& vars,
00853                                  const BCP_vec<BCP_cut*>& cuts,
00854                                  const bool before_fathom,
00855                                  BCP_vec<BCP_var*>& new_vars,
00856                                  BCP_vec<BCP_col*>& new_cols)
00857 {
00858     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00859           "LP: Default generate_vars_in_lp() executed.\n");
00860 }
00861 //-----------------------------------------------------------------------------
00862 BCP_object_compare_result
00863 BCP_lp_user::compare_cuts(const BCP_cut* c0, const BCP_cut* c1)
00864 {
00865     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00866           "LP: Default compare_cuts() executed.\n");
00867     return BCP_DifferentObjs;
00868 }
00869 //-----------------------------------------------------------------------------
00870 BCP_object_compare_result
00871 BCP_lp_user::compare_vars(const BCP_var* v0, const BCP_var* v1)
00872 {
00873     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00874           "LP: Default compare_vars() executed.\n");
00875     return BCP_DifferentObjs;
00876 }
00877 
00878 //#############################################################################
00879 
00880 void
00881 BCP_lp_user::select_vars_to_delete(const BCP_lp_result& lpres,
00882                                    const BCP_vec<BCP_var*>& vars,
00883                                    const BCP_vec<BCP_cut*>& cuts,
00884                                    const bool before_fathom,
00885                                    BCP_vec<int>& deletable)
00886 {
00887     if (before_fathom && p->param(BCP_lp_par::NoCompressionAtFathom))
00888         return;
00889     const int varnum = vars.size();
00890     deletable.reserve(varnum);
00891     for (int i = p->core->varnum(); i < varnum; ++i) {
00892         BCP_var *var = vars[i];
00893         if (var->is_to_be_removed() ||
00894             (! var->is_non_removable() && var->lb() == 0 && var->ub() == 0)) {
00895             deletable.unchecked_push_back(i);
00896         }
00897     }
00898 }
00899 
00900 //=============================================================================
00901 
00902 void
00903 BCP_lp_user::select_cuts_to_delete(const BCP_lp_result& lpres,
00904                                    const BCP_vec<BCP_var*>& vars,
00905                                    const BCP_vec<BCP_cut*>& cuts,
00906                                    const bool before_fathom,
00907                                    BCP_vec<int>& deletable)
00908 {
00909     if (before_fathom && p->param(BCP_lp_par::NoCompressionAtFathom))
00910         return;
00911     const int cutnum = cuts.size();
00912     const int ineff_to_delete = p->param(BCP_lp_par::IneffectiveBeforeDelete);
00913     const double lb = lpres.objval();
00914     const BCP_vec<double>& lb_at_cutgen = p->node->lb_at_cutgen;
00915     deletable.reserve(cutnum);
00916     for (int i = p->core->cutnum(); i < cutnum; ++i) {
00917         BCP_cut *cut = cuts[i];
00918         if (cut->is_to_be_removed() ||
00919             (! cut->is_non_removable() &&
00920              cut->effective_count() <= -ineff_to_delete &&
00921              lb_at_cutgen[i] < lb - 0.0001)) {
00922             deletable.unchecked_push_back(i);
00923         }
00924     }
00925 }
00926 
00927 //#############################################################################
00928 
00929 void
00930 BCP_lp_user::logical_fixing(const BCP_lp_result& lpres,
00931                             const BCP_vec<BCP_var*>& vars,
00932                             const BCP_vec<BCP_cut*>& cuts,
00933                             const BCP_vec<BCP_obj_status>& var_status,
00934                             const BCP_vec<BCP_obj_status>& cut_status,
00935                             const int var_bound_changes_since_logical_fixing,
00936                             BCP_vec<int>& changed_pos,
00937                             BCP_vec<double>& new_bd)
00938 {
00939     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
00940           "LP: Default logical_fixing() executed.\n");
00941 }
00942 
00943 //#############################################################################
00944 void
00945 BCP_lp_user::reduced_cost_fixing(const double* dj, const double* x,
00946                                  const double gap,
00947                                  BCP_vec<BCP_var*>& vars, int& newly_changed)
00948 {
00949     OsiBabSolver* babSolver = getOsiBabSolver();
00950     if (babSolver && !babSolver->reducedCostsAccurate())
00951         return;
00952 
00953     newly_changed = 0;
00954     const bool atZero = get_param(BCP_lp_par::DoReducedCostFixingAtZero);
00955     const bool atAny = get_param(BCP_lp_par::DoReducedCostFixingAtAnything);
00956 
00957     if (! atZero && ! atAny)
00958         return;
00959 
00960     double petol = 0.0;
00961     p->lp_solver->getDblParam(OsiPrimalTolerance, petol);
00962     double detol = 0.0;
00963     p->lp_solver->getDblParam(OsiDualTolerance, detol);
00964 
00965     // If the gap is negative that means that we are above the limit, so
00966     // don't do anything.
00967     if (gap < 0)
00968         return;
00969 
00970     const int varnum = vars.size();
00971     BCP_vec<int> changed_indices;
00972     changed_indices.reserve(varnum);
00973     BCP_vec<double> changed_bounds;
00974     changed_bounds.reserve(2 * varnum);
00975 
00976     // Note that when this function is called, we must have a dual
00977     // feasible dual solution. Therefore we can use the lagrangean
00978     // relaxation to fix variables.
00979 
00980     // *FIXME* : If we knew that there are integral vars only, then
00981     // we could leave out the test for BCP_ContinuousVar...
00982     for (int i = 0; i < varnum; ++i) {
00983         BCP_var* var = vars[i];
00984         if (! var->is_fixed() && var->var_type() != BCP_ContinuousVar){
00985             if (dj[i] > detol) {
00986                 const double lb = var->lb();
00987                 const double new_ub = lb + floor(gap / dj[i]);
00988                 if (new_ub < var->ub() && (atAny || CoinAbs(x[i])<petol) ) {
00989                     vars[i]->set_ub(new_ub);
00990                     changed_indices.unchecked_push_back(i);
00991                     changed_bounds.unchecked_push_back(lb);
00992                     changed_bounds.unchecked_push_back(new_ub);
00993                 }
00994             } else if (dj[i] < -detol) {
00995                 const double ub = var->ub();
00996                 const double new_lb = ub - floor(gap / (-dj[i]));
00997                 if (new_lb > var->lb() && (atAny || CoinAbs(x[i])<petol) ) {
00998                     vars[i]->set_lb(new_lb);
00999                     changed_indices.unchecked_push_back(i);
01000                     changed_bounds.unchecked_push_back(new_lb);
01001                     changed_bounds.unchecked_push_back(ub);
01002                 }
01003             }
01004         }
01005     }
01006     newly_changed = changed_indices.size();
01007     if (newly_changed > 0) {
01008         p->lp_solver->setColSetBounds(changed_indices.begin(),
01009                                       changed_indices.end(),
01010                                       changed_bounds.begin());
01011     }
01012 }
01013 
01014 //#############################################################################
01015 //#############################################################################
01016 
01017 int
01018 BCP_lp_user::try_to_branch(OsiBranchingInformation& branchInfo,
01019                            OsiSolverInterface* solver,
01020                            OsiChooseVariable* choose,
01021                            OsiBranchingObject*& branchObject,
01022                            bool allowVarFix)
01023 {
01024     int returnStatus = 0;
01025     int numUnsatisfied = choose->setupList(&branchInfo, true);
01026     choose->setBestObjectIndex(-1);
01027     if (numUnsatisfied > 0) {
01028         if (choose->numberOnList() == 0) {
01029             // Nothing left for strong branching to evaluate
01030             if (choose->numberOnList() > 0 || choose->numberStrong() == 0) {
01031                 // There is something on the list
01032                 choose->setBestObjectIndex(choose->candidates()[0]);
01033             } else {
01034                 // There is nothing on the list
01035                 numUnsatisfied = choose->setupList(&branchInfo, false);
01036                 if (numUnsatisfied > 0) {
01037                     choose->setBestObjectIndex(choose->candidates()[0]);
01038                 }
01039             }
01040         } else {
01041             // Do the strong branching
01042             int ret = choose->chooseVariable(solver, &branchInfo, allowVarFix);
01043             /* Check if SB has fixed anything, and if so, apply the same
01044                changes to the vars */
01045             const double * clb = solver->getColLower();
01046             const double * cub = solver->getColUpper();
01047             BCP_vec<BCP_var*>& vars = p->node->vars;
01048             for (int i = numUnsatisfied - 1; i >= 0; --i) {
01049                 const int ind =
01050                   solver->object(choose->candidates()[i])->columnNumber();
01051                 if (ind >= 0) {
01052                   assert(vars[ind]->lb() <= clb[ind]);
01053                   assert(vars[ind]->ub() >= cub[ind]);
01054                   vars[ind]->set_lb_ub(clb[ind], cub[ind]);
01055                 }
01056             }
01057                 
01058             /* update number of strong iterations etc
01059             model->incrementStrongInfo(choose->numberStrongDone(),
01060                                        choose->numberStrongIterations(),
01061                                        ret==-1 ? 0:choose->numberStrongFixed(),
01062                                        ret==-1);
01063             */
01064             if (ret > 1) {
01065                 // has fixed some
01066                 returnStatus = -1;
01067             } else if (ret == -1) {
01068                 // infeasible
01069                 returnStatus = -2;
01070             } else if (ret == 0) {
01071                 // normal
01072                 returnStatus = 0;
01073                 numUnsatisfied = 1;
01074             } else {
01075                 // ones on list satisfied - double check
01076                 numUnsatisfied = choose->setupList(&branchInfo, false);
01077                 if (numUnsatisfied > 0) {
01078                     choose->setBestObjectIndex(choose->candidates()[0]);
01079                 }
01080             }
01081         }
01082     }
01083     if (! returnStatus) {
01084         if (numUnsatisfied > 0) {
01085             // create branching object
01086             /* FIXME: how the objects are created? */
01087             const OsiObject * obj = solver->object(choose->bestObjectIndex());
01088             branchObject = obj->createBranch(solver,
01089                                              &branchInfo,
01090                                              obj->whichWay());
01091         }
01092     }
01093 
01094     return returnStatus;
01095 }
01096 
01097 //#############################################################################
01098 
01099 BCP_branching_decision BCP_lp_user::
01100 select_branching_candidates(const BCP_lp_result& lpres,
01101                             const BCP_vec<BCP_var*>& vars,
01102                             const BCP_vec<BCP_cut*>& cuts,
01103                             const BCP_lp_var_pool& local_var_pool,
01104                             const BCP_lp_cut_pool& local_cut_pool,
01105                             BCP_vec<BCP_lp_branching_object*>& cands,
01106                             bool force_branch)
01107 {
01108     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
01109           "LP: Default select_branching_candidates() executed.\n");
01110 
01111     if (lpres.termcode() & BCP_Abandoned) {
01112       // *THINK*: maybe we should branch through...
01113       print(true, "\
01114 LP: ############ LP solver abandoned. Branching through...\n");
01115     }
01116 
01117     // *THINK* : this branching object selection is very primitive
01118     // *THINK* : should check for tail-off, could check for branching cuts, etc
01119     if (! force_branch &&
01120         (local_var_pool.size() > 0 || local_cut_pool.size() > 0))
01121         return BCP_DoNotBranch;
01122 
01123     OsiSolverInterface* lp = p->lp_solver;
01124 
01125     /* The last true: make a copy in brInfo of the sol of solver
01126        instead of just pointing into the solver's copy */
01127     OsiBranchingInformation brInfo(lp, true, true);
01128     lp->getDblParam(OsiDualObjectiveLimit, brInfo.cutoff_);
01129     brInfo.integerTolerance_ = p->param(BCP_lp_par::IntegerTolerance);
01130     brInfo.timeRemaining_ = get_param(BCP_lp_par::MaxRunTime) - CoinCpuTime();
01131     brInfo.numberSolutions_ = 0; /*FIXME*/
01132     brInfo.numberBranchingSolutions_ = 0; /*FIXME numBranchingSolutions_;*/
01133     brInfo.depth_ = current_level();
01134 
01135     OsiChooseStrong* strong = new OsiChooseStrong(lp);
01136     strong->setNumberBeforeTrusted(5); // the default in Cbc
01137     strong->setNumberStrong(p->param(BCP_lp_par::StrongBranchNum));
01138     strong->setTrustStrongForSolution(false);
01144     strong->setShadowPriceMode(0);
01145 
01146     OsiChooseVariable * choose = strong;
01147     OsiBranchingObject* brObj = NULL;
01148 
01149     bool allowVarFix = true;
01150     /*
01151       <li>  0: A branching object has been installed
01152       <li> -1: A monotone object was discovered
01153       <li> -2: An infeasible object was discovered
01154     */
01155     int brResult =
01156         try_to_branch(brInfo, lp, choose, brObj, allowVarFix);
01157     const int bestWhichWay = choose->bestWhichWay();
01158 
01159 #if 0
01160     /* FIXME:before doing anything check if we have found a new solution */
01161     if (choose->goodSolution()
01162         &&model->problemFeasibility()->feasible(model,-1)>=0) {
01163         // yes
01164         double objValue = choose->goodObjectiveValue();
01165         model->setBestSolution(CBC_STRONGSOL,
01166                                objValue,
01167                                choose->goodSolution()) ;
01168         model->setLastHeuristic(NULL);
01169         model->incrementUsed(choose->goodSolution());
01170         choose->clearGoodSolution();
01171     }
01172 #endif
01173 
01174     delete choose;
01175 
01176     switch (brResult) {
01177     case -2:
01178         // when doing strong branching a candidate has proved that the
01179         // problem is infeasible
01180         delete brObj;
01181         return BCP_DoNotBranch_Fathomed;
01182     case -1:
01183         // OsiChooseVariable::chooseVariable() returned 2, 3, or 4
01184         if (!brObj) {
01185             // just go back and resolve
01186             return BCP_DoNotBranch;
01187         }
01188         // otherwise might as well branch. The forced variable is
01189         // unlikely to jump up one more (though who knows...)
01190         break;
01191     case 0:
01192         if (!brObj) {
01193             // nothing got fixed, yet couldn't find something to branch on
01194             throw BCP_fatal_error("BM: Couldn't branch!\n");
01195         }
01196         // we've got a branching object
01197         break;
01198     default:
01199         throw BCP_fatal_error("\
01200 BCP: BCP_lp_user::try_to_branch returned with unknown return code.\n");
01201     }
01202 
01203     // If there were some fixings (brResult < 0) then propagate them where
01204     // needed
01205     if (allowVarFix && brResult < 0) {
01206         const double* clb = lp->getColLower();
01207         const double* cub = lp->getColUpper();
01208         /* There may or may not have been changes, but faster to set then to
01209            test... */
01210         BCP_vec<BCP_var*>& vars = p->node->vars;
01211         int numvars = vars.size();
01212         for (int i = 0; i < numvars; ++i) {
01213             vars[i]->change_bounds(clb[i], cub[i]);
01214         }
01215     }
01216     
01217     // all possibilities are 2-way branches
01218     int order[2] = {0, 1};
01219     if (bestWhichWay == 1) {
01220         order[0] = 1;
01221         order[1] = 0;
01222     }
01223 
01224     // Now interpret the result (at this point we must have a brObj
01225     OsiIntegerBranchingObject* intBrObj =
01226         dynamic_cast<OsiIntegerBranchingObject*>(brObj);
01227     if (intBrObj) {
01228         BCP_lp_integer_branching_object o(intBrObj);
01229         cands.push_back(new BCP_lp_branching_object(o, order));
01230     }
01231     OsiSOSBranchingObject* sosBrObj =
01232         dynamic_cast<OsiSOSBranchingObject*>(brObj);
01233     if (sosBrObj) {
01234         BCP_lp_sos_branching_object o(sosBrObj);
01235         cands.push_back(new BCP_lp_branching_object(lp, o, order));
01236     }
01237 
01238     delete brObj;
01239     
01240     if (cands.size() == 0) {
01241         throw BCP_fatal_error("\
01242  LP : No var/cut in pool but couldn't select branching object.\n");
01243     }
01244     return BCP_DoBranch;
01245 }
01246 
01247 //-----------------------------------------------------------------------------
01248 void
01249 BCP_lp_user::append_branching_vars(const double* x,
01250                                    const BCP_vec<BCP_var*>& vars,
01251                                    const BCP_vec<int>& select_pos,
01252                                    BCP_vec<BCP_lp_branching_object*>& cans)
01253 {
01254     BCP_vec<int>::const_iterator spi;
01255     BCP_vec<double> vbd(4, 0.0);
01256     BCP_vec<int> vpos(1, 0);
01257 
01258     for (spi = select_pos.begin(); spi != select_pos.end(); ++spi) {
01259         const int pos = *spi;
01260         vpos[0] = pos;
01261         vbd[0] = vars[pos]->lb();
01262         vbd[1] = floor(x[pos]);
01263         vbd[2] = ceil(x[pos]);
01264         vbd[3] = vars[pos]->ub();
01265         cans.push_back
01266             (new  BCP_lp_branching_object(2,
01267                                           0, 0, /* vars/cuts_to_add */
01268                                           &vpos, 0, &vbd, 0, /* forced parts */
01269                                           0, 0, 0, 0 /* implied parts */));
01270     }
01271 }
01272 
01273 //-----------------------------------------------------------------------------
01274 BCP_branching_object_relation BCP_lp_user::
01275 compare_branching_candidates(BCP_presolved_lp_brobj* new_presolved,
01276                              BCP_presolved_lp_brobj* old_presolved)
01277 {
01278     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
01279           "LP: Default compare_presolved_branching_objects() executed.\n");
01280 
01281     // change the objvals according to the termcodes (in order to be able to
01282     // make decisions based on objval, no matter what the termcode is).
01283     new_presolved->fake_objective_values(p->lp_result->objval());
01284 
01285     // If using this branching object we can fathom all children, then choose
01286     // this object
01287     if (new_presolved->fathomable(p->ub() - p->granularity()))
01288         return(BCP_NewPresolvedIsBetter_BranchOnIt);
01289 
01290     // If this is the first, keep it
01291     if (! old_presolved)
01292         return(BCP_NewPresolvedIsBetter);
01293 
01294     // If something had gone wrong with at least one descendant in the new then
01295     // we prefer to keep the old.
01296     if (new_presolved->had_numerical_problems())
01297         return(BCP_OldPresolvedIsBetter);
01298 
01299     // OK, so all descendants finished just fine. Do whatever the parameter
01300     // says
01301 
01302     if (p->param(BCP_lp_par::BranchingObjectComparison)&BCP_Comparison_Objval){
01303         const double objetol = 1e-7;
01304       
01305         BCP_vec<double> new_obj;
01306         new_presolved->get_lower_bounds(new_obj);
01307         std::sort(new_obj.begin(), new_obj.end());
01308         const int new_not_fathomed =
01309             new_obj.index(std::lower_bound(new_obj.begin(), new_obj.end(),
01310                                            BCP_DBL_MAX / 10));
01311 
01312         BCP_vec<double> old_obj;
01313         old_presolved->get_lower_bounds(old_obj);
01314         std::sort(old_obj.begin(), old_obj.end());
01315         const int old_not_fathomed =
01316             old_obj.index(std::lower_bound(old_obj.begin(), old_obj.end(),
01317                                            BCP_DBL_MAX / 10));
01318 
01319         if (new_not_fathomed < old_not_fathomed)
01320             return BCP_NewPresolvedIsBetter;
01321 
01322         if (new_not_fathomed > old_not_fathomed)
01323             return BCP_OldPresolvedIsBetter;
01324 
01325         BCP_vec<double>::const_iterator new_first = new_obj.begin();
01326         BCP_vec<double>::const_iterator new_last = new_obj.end();
01327         BCP_vec<double>::const_iterator old_first = old_obj.begin();
01328         BCP_vec<double>::const_iterator old_last = old_obj.end();
01329    
01330         switch(p->param(BCP_lp_par::BranchingObjectComparison)){
01331         case BCP_LowestAverageObjval:
01332         case BCP_HighestAverageObjval:
01333             {
01334                 double newavg = 0;
01335                 for ( ; new_first != new_last; ++new_first) {
01336                     if (*new_first < BCP_DBL_MAX / 10)
01337                         newavg += *new_first;
01338                 }
01339                 newavg /= new_not_fathomed;
01340                 double oldavg = 0;
01341                 for ( ; old_first != old_last; ++old_first) {
01342                     if (*old_first < BCP_DBL_MAX / 10)
01343                         oldavg += *old_first;
01344                 }
01345                 oldavg /= old_not_fathomed;
01346                 const bool high =
01347                     (p->param(BCP_lp_par::BranchingObjectComparison)
01348                      == BCP_HighestAverageObjval);
01349                 return (high && newavg > oldavg) || (!high && newavg < oldavg)?
01350                     BCP_NewPresolvedIsBetter : BCP_OldPresolvedIsBetter;
01351             }
01352 
01353         case BCP_LowestLowObjval:
01354             for ( ; new_first != new_last && old_first != old_last;
01355                   ++new_first, ++old_first) {
01356                 if (*new_first < *old_first - objetol)
01357                     return BCP_NewPresolvedIsBetter;
01358                 if (*new_first > *old_first + objetol)
01359                     return BCP_OldPresolvedIsBetter;
01360             }
01361             return old_first == old_last ?
01362                 BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;
01363 
01364         case BCP_HighestLowObjval:
01365             for ( ; new_first != new_last && old_first != old_last;
01366                   ++new_first, ++old_first) {
01367                 if (*new_first > *old_first + objetol)
01368                     return BCP_NewPresolvedIsBetter;
01369                 if (*new_first < *old_first - objetol)
01370                     return BCP_OldPresolvedIsBetter;
01371             }
01372             return old_first == old_last ?
01373                 BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;
01374 
01375         case BCP_LowestHighObjval:
01376             while (new_first != new_last && old_first != old_last) {
01377                 --new_last;
01378                 --old_last;
01379                 if (*new_last < *old_last - objetol)
01380                     return BCP_NewPresolvedIsBetter;
01381                 if (*new_last > *old_last + objetol)
01382                     return BCP_OldPresolvedIsBetter;
01383             }
01384             return old_first == old_last ?
01385                 BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;
01386 
01387         case BCP_HighestHighObjval:
01388             while (new_first != new_last && old_first != old_last) {
01389                 --new_last;
01390                 --old_last;
01391                 if (*new_last > *old_last + objetol)
01392                     return BCP_NewPresolvedIsBetter;
01393                 if (*new_last < *old_last - objetol)
01394                     return BCP_OldPresolvedIsBetter;
01395             }
01396             return old_first == old_last ?
01397                 BCP_OldPresolvedIsBetter : BCP_NewPresolvedIsBetter;
01398         default:
01399             throw BCP_fatal_error("\
01400 Unknown branching object comparison rule.\n");
01401         }
01402     }
01403 
01404     // shouldn't ever get here
01405     throw BCP_fatal_error("No branching object comparison rule is chosen.\n");
01406 
01407     // fake return
01408     return BCP_NewPresolvedIsBetter;
01409 }
01410 
01411 //-----------------------------------------------------------------------------
01412 void
01413 BCP_lp_user::set_actions_for_children(BCP_presolved_lp_brobj* best)
01414 {
01415     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
01416           "LP: Default set_actions_for_children() executed.\n");
01417 
01418     // by default every action is set to BCP_ReturnChild
01419     if (p->node->dive == BCP_DoNotDive)
01420         return;
01421    
01422     BCP_vec<BCP_child_action>& action = best->action();
01423     int i, ind;
01424 
01425     switch (p->param(BCP_lp_par::ChildPreference)){
01426     case BCP_PreferDiveDown:
01427         ind = 0;
01428         break;
01429 
01430     case BCP_PreferDiveUp:
01431         ind = best->candidate()->child_num - 1;
01432         break;
01433                 
01434     case BCP_PreferChild_LowBound:
01435         // NOTE: if the lowest objval child is fathomed then everything is
01436         for (ind = 0, i = best->candidate()->child_num - 1; i > 0; --i){
01437             if (best->lpres(i).objval() < best->lpres(ind).objval())
01438                 ind = i;
01439         }
01440         break;
01441 
01442     case BCP_PreferChild_HighBound:
01443         // NOTE: this selects the highest objval child NOT FATHOMED, thus if
01444         // the highest objval child is fathomed then everything is
01445         for (ind = 0, i = best->candidate()->child_num - 1; i > 0; --i) {
01446             if (! p->over_ub(best->lpres(i).objval()) &&
01447                 best->lpres(i).objval() < best->lpres(ind).objval())
01448                 ind = i;
01449         }
01450         break;
01451     default:
01452         // *THINK* : fractional branching
01453         throw BCP_fatal_error("LP : Unrecognized child preference.\n");
01454     }
01455     // mark the ind-th to be kept (if not over ub)
01456     if (! p->over_ub(best->lpres(ind).objval()))
01457         action[ind] = BCP_KeepChild;
01458 }
01459 
01460 //#############################################################################
01461 
01462 void
01463 BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best, 
01464                                         const int selected)
01465 {
01466     using_deprecated_set_user_data_for_children = true;
01467     set_user_data_for_children(best);
01468     print(using_deprecated_set_user_data_for_children,
01469           "\
01470 *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING ***\n\
01471 You have overridden\n\
01472   BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best)\n\
01473 which is a deprecated virtual function. Please override\n\
01474   BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best,\n\
01475                                           const int selected)\n\
01476 instead. The old version will go away, your code will still compile, but it\n\
01477 will not do what you intend it to be doing.\n\
01478 *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING ***\n"\
01479                );
01480 }
01481 
01482 //#############################################################################
01483 
01484 void
01485 BCP_lp_user::set_user_data_for_children(BCP_presolved_lp_brobj* best)
01486 {
01487     using_deprecated_set_user_data_for_children = false;
01488     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
01489           "LP: Default set_user_data_for_children() executed.\n");
01490 }
01491 
01492 //#############################################################################
01493 // purging the slack cut pool (candidates for branching on cut)
01494 void
01495 BCP_lp_user::purge_slack_pool(const BCP_vec<BCP_cut*>& slack_pool,
01496                               BCP_vec<int>& to_be_purged)
01497 {
01498     print(p->param(BCP_lp_par::ReportWhenDefaultIsExecuted),
01499           "LP: Default purge_slack_pool() executed.\n");
01500 
01501     switch (p->param(BCP_lp_par::SlackCutDiscardingStrategy)) {
01502     case BCP_DiscardSlackCutsAtNewNode:
01503         if (current_iteration() != 1)
01504             break;
01505     case BCP_DiscardSlackCutsAtNewIteration:
01506         {
01507             const int size = slack_pool.size();
01508             if (size > 0) {
01509                 to_be_purged.reserve(size);
01510                 for (int i = 0; i < size; ++i)
01511                     to_be_purged.unchecked_push_back(i);
01512             }
01513         }
01514         break;
01515     }
01516 }
01517 
01518 //#############################################################################

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