/home/coin/SVN-release/OS-2.0.0/Bcp/examples/MCF-3/LP/MCF3_lp.cpp

Go to the documentation of this file.
00001 #include "CoinHelperFunctions.hpp"
00002 #include "OsiClpSolverInterface.hpp"
00003 #include "MCF3_lp.hpp"
00004 
00005 //#############################################################################
00006 
00007 OsiSolverInterface* MCF3_lp::initialize_solver_interface()
00008 {
00009   return new OsiClpSolverInterface;
00010 }
00011 
00012 //#############################################################################
00013 
00014 void
00015 MCF3_adjust_bounds(const BCP_vec<BCP_var*>& vars,
00016                    std::vector<MCF3_branch_decision>* history,
00017                    BCP_vec<int>& var_changed_pos,
00018                    BCP_vec<double>& var_new_bd)
00019 {
00020   for (int i = vars.size()-1; i >= 0; --i) {
00021     MCF3_var* v = dynamic_cast<MCF3_var*>(vars[i]);
00022     if (!v) continue;
00023     const int vsize = v->flow.getNumElements();
00024     const int* vind = v->flow.getIndices();
00025     const double* vval = v->flow.getElements();
00026 
00027     bool violated = false;
00028     const std::vector<MCF3_branch_decision>& hist = history[v->commodity];
00029     for (int j = hist.size()-1; !violated && j >= 0; --j) {
00030       const MCF3_branch_decision& h = hist[j];
00031       double f = 0.0;
00032       for (int k = 0; k < vsize; ++k) {
00033         if (vind[k] == h.arc_index) {
00034           f = vval[k];
00035           break;
00036         }
00037       }
00038       violated = (f < h.lb || f > h.ub);
00039     }
00040     if (violated) {
00041       var_changed_pos.push_back(i);
00042       var_new_bd.push_back(0.0); // the new lower bound on the var
00043       var_new_bd.push_back(0.0); // the new upper bound on the var
00044     }
00045   }
00046 }
00047 
00048 /*---------------------------------------------------------------------------*/
00049 
00050 void
00051 MCF3_adjust_bounds(const BCP_vec<BCP_var*>& vars,
00052                    const int commodity,
00053                    const MCF3_branch_decision& decision,
00054                    BCP_vec<int>& var_changed_pos,
00055                    BCP_vec<double>& var_new_bd)
00056 {
00057   for (int i = vars.size()-1; i >= 0; --i) {
00058     MCF3_var* v = dynamic_cast<MCF3_var*>(vars[i]);
00059     if (!v)
00060       continue;
00061     if (v->commodity != commodity)
00062       continue;
00063     const int vsize = v->flow.getNumElements();
00064     const int* vind = v->flow.getIndices();
00065     const double* vval = v->flow.getElements();
00066 
00067     double f = 0.0;
00068     for (int k = 0; k < vsize; ++k) {
00069       if (vind[k] == decision.arc_index) {
00070         f = vval[k];
00071         break;
00072       }
00073     }
00074     if (f < decision.lb || f > decision.ub) {
00075       var_changed_pos.push_back(i);
00076       var_new_bd.push_back(0.0); // the new lower bound on the var
00077       var_new_bd.push_back(0.0); // the new upper bound on the var
00078     }
00079   }
00080 }
00081 
00082 /*---------------------------------------------------------------------------*/
00083 
00084 void MCF3_lp::
00085 initialize_new_search_tree_node(const BCP_vec<BCP_var*>& vars,
00086                                 const BCP_vec<BCP_cut*>& cuts,
00087                                 const BCP_vec<BCP_obj_status>& var_status,
00088                                 const BCP_vec<BCP_obj_status>& cut_status,
00089                                 BCP_vec<int>& var_changed_pos,
00090                                 BCP_vec<double>& var_new_bd,
00091                                 BCP_vec<int>& cut_changed_pos,
00092                                 BCP_vec<double>& cut_new_bd)
00093 {
00094   branch_history = dynamic_cast<MCF3_user*>(get_user_data())->branch_history;
00095 
00096   /* The branching decisions may set bounds that are violated by some
00097      of the variables (more precisely by the flow they represent). Find them
00098      and set their upper bounds to 0, since it's highly unlikely that we'd
00099      be able to find a fractional \lambda vector such that:
00100      1) such a flow has nonzero coefficient;
00101      2) the convex combination of all flows with nonzero \lambda is integer.
00102   */
00103   MCF3_adjust_bounds(vars, branch_history, var_changed_pos, var_new_bd);
00104 
00105   // clear out our local pool
00106   purge_ptr_vector(gen_vars);
00107 }
00108 
00109 /*---------------------------------------------------------------------------*/
00110 
00111 BCP_solution* MCF3_lp::
00112 test_feasibility(const BCP_lp_result& lpres,
00113                  const BCP_vec<BCP_var*>& vars,
00114                  const BCP_vec<BCP_cut*>& cuts)
00115 {
00116 #if 0
00117   static int cnt = 0;
00118   char name[100];
00119   sprintf(name, "currlp-%i", cnt++);
00120   //  sprintf(name, "currlp");
00121   getLpProblemPointer()->lp_solver->writeMps(name, "mps");
00122   printf("Current LP written in file %s.mps\n", name);
00123   getLpProblemPointer()->lp_solver->writeLp(name, "lp");
00124   printf("Current LP written in file %s.lp\n", name);
00125 #endif
00126 
00127   // Feasibility testing: we need to test whether the convex combination of
00128   // the current columns (according to \lambda, the current primal solution)
00129   // is integer feasible for the original problem. The only thing that can
00130   // be violated is integrality.
00131   
00132   int i, j;
00133   for (i = data.numcommodities-1; i >= 0; --i) {
00134     flows[i].clear();
00135   }
00136 
00137   const double* x = lpres.x();
00138   for (i = vars.size()-1; i >= 0; --i) {
00139     MCF3_var* v = dynamic_cast<MCF3_var*>(vars[i]);
00140     if (!v) continue;
00141     std::map<int,double>& f = flows[v->commodity];
00142     const int vsize = v->flow.getNumElements();
00143     const int* vind = v->flow.getIndices();
00144     const double* vval = v->flow.getElements();
00145     for (j = 0; j < vsize; ++j) {
00146       f[vind[j]] += vval[j]*x[i];
00147     }
00148   }
00149   for (i = data.numcommodities-1; i >= 0; --i) {
00150     std::map<int,double>& f = flows[i];
00151     for (std::map<int,double>::iterator fi=f.begin(); fi != f.end(); ++fi) {
00152       const double frac = fabs(fi->second - floor(fi->second) - 0.5);
00153       if (frac < 0.5-1e-6)
00154         return NULL;
00155     }
00156   }
00157 
00158   // Found an integer solution
00159   BCP_solution_generic* gsol = new BCP_solution_generic;
00160   for (i = data.numcommodities-1; i >= 0; --i) {
00161     std::map<int,double>& f = flows[i];
00162     CoinPackedVector flow(false);
00163     double weight = 0;
00164     for (std::map<int,double>::iterator fi=f.begin();
00165          fi != f.end();
00166          ++fi) {
00167       const double val = floor(fi->second + 0.5);
00168       if (val > 0) {
00169         flow.insert(fi->first, val);
00170         weight += val * data.arcs[fi->first].weight;
00171       }
00172     }
00173     gsol->add_entry(new MCF3_var(i, flow, weight), 1.0);
00174   }
00175   return gsol;
00176 }
00177 
00178 /*---------------------------------------------------------------------------*/
00179 
00180 double MCF3_lp::
00181 compute_lower_bound(const double old_lower_bound,
00182                     const BCP_lp_result& lpres,
00183                     const BCP_vec<BCP_var*>& vars,
00184                     const BCP_vec<BCP_cut*>& cuts)
00185 {
00186   // To compute a true lower bound we need to generate variables first (if
00187   // we can). These are saved so that we can return them in
00188   // generate_vars_in_lp.
00189 
00190   // generate variables:  for each commodity generate a flow.
00191 
00192   const int numarcs = data.numarcs;
00193   double* cost = new double[numarcs];
00194   const double* pi = lpres.pi();
00195   const double* nu = pi + numarcs;
00196 
00197   int i, j;
00198 
00199   for (i = numarcs-1; i >= 0; --i) {
00200     cost[i] = data.arcs[i].weight - pi[i];
00201   }
00202   cg_lp->setObjective(cost);
00203 
00204   // This will hold generated variables
00205   int* ind = new int[numarcs];
00206   double* val = new double[numarcs];
00207   int cnt = 0;
00208 
00209   for (i = data.numcommodities-1; i >= 0; --i) {
00210     const MCF3_data::commodity& comm = data.commodities[i];
00211     cg_lp->setRowBounds(comm.source, -comm.demand, -comm.demand);
00212     cg_lp->setRowBounds(comm.sink, comm.demand, comm.demand);
00213     const std::vector<MCF3_branch_decision>& hist = branch_history[i];
00214     for (j = hist.size() - 1; j >= 0; --j) {
00215       const MCF3_branch_decision& h = hist[j];
00216       cg_lp->setColBounds(h.arc_index, h.lb, h.ub);
00217     }
00218     cg_lp->initialSolve();
00219     if (cg_lp->isProvenOptimal() && cg_lp->getObjValue() < nu[i] - 1e-8) {
00220       // we have generated a column. Create a var out of it. Round the
00221       // double values while we are here, after all, they should be
00222       // integer. there can only be some tiny roundoff error by the LP
00223       // solver
00224       const double* x = cg_lp->getColSolution();
00225       cnt = 0;
00226       double obj = 0.0;
00227       for (int j = 0; j < numarcs; ++j) {
00228         const double xval = floor(x[j] + 0.5);
00229         if (xval != 0.0) {
00230           ind[cnt] = j;
00231           val[cnt] = xval;
00232           ++cnt;
00233           obj += data.arcs[j].weight * xval;
00234         }
00235       }
00236       gen_vars.push_back(new MCF3_var(i, CoinPackedVector(cnt, ind, val,
00237                                                           false), obj));
00238     }
00239     for (j = hist.size() - 1; j >= 0; --j) {
00240       const int ind = hist[j].arc_index;
00241       cg_lp->setColBounds(ind, data.arcs[ind].lb, data.arcs[ind].ub);
00242     }
00243     cg_lp->setRowBounds(comm.source, 0, 0);
00244     cg_lp->setRowBounds(comm.sink, 0, 0);
00245   }
00246   delete[] val;
00247   delete[] ind;
00248   delete[] cost;
00249 
00250   // Excercise: do the same in a random order and apply dual discounting
00251   // Not yet implemented.
00252 
00253   // Now get a true lower bound
00254   double true_lower_bound = 0.0;
00255   generated_vars = (gen_vars.size() > 0);
00256 
00257   if (generated_vars) {
00258     true_lower_bound = old_lower_bound;
00259   } else {
00260     true_lower_bound = lpres.objval();
00261   }
00262 
00263   // Excercise: Get a better true lower bound
00264   // Hint: lpres.objval() + The sum of the reduced costs of the
00265   // variables with the most negative reduced cost in each subproblem
00266   // yield a true lower bound
00267   // Not yet implemented.
00268 
00269   return true_lower_bound;
00270 }
00271 
00272 /*---------------------------------------------------------------------------*/
00273 
00274 void MCF3_lp::
00275 generate_vars_in_lp(const BCP_lp_result& lpres,
00276                     const BCP_vec<BCP_var*>& vars,
00277                     const BCP_vec<BCP_cut*>& cuts,
00278                     const bool before_fathom,
00279                     BCP_vec<BCP_var*>& new_vars,
00280                     BCP_vec<BCP_col*>& new_cols)
00281 {
00282   new_vars.append(gen_vars);
00283   gen_vars.clear();
00284 }
00285 
00286 /*---------------------------------------------------------------------------*/
00287 
00288 void MCF3_lp::
00289 vars_to_cols(const BCP_vec<BCP_cut*>& cuts,
00290              BCP_vec<BCP_var*>& vars,
00291              BCP_vec<BCP_col*>& cols,
00292              const BCP_lp_result& lpres,
00293              BCP_object_origin origin, bool allow_multiple)
00294 {
00295   static const CoinPackedVector emptyVector(false);
00296   const int numvars = vars.size();
00297   for (int i = 0; i < numvars; ++i) {
00298     const MCF3_var* v =
00299       dynamic_cast<const MCF3_var*>(vars[i]);
00300     if (v) {
00301       // Since we do not generate cuts, we can just disregard the "cuts"
00302       // argument, since the column corresponding to the var is exactly
00303       // the flow (plus the entry in the appropriate convexity
00304       // constraint)
00305       BCP_col* col = new BCP_col(v->flow, v->weight, v->lb(), v->ub());
00306       col->insert(data.numarcs + v->commodity, 1.0);
00307       cols.push_back(col);
00308       // Excercise: if we had generated cuts, then the coefficients for
00309       // those rows would be appended to the end of each column
00310       continue;
00311     }
00312   }
00313 }
00314 
00315 /*---------------------------------------------------------------------------*/
00316 
00317 BCP_branching_decision MCF3_lp::
00318 select_branching_candidates(const BCP_lp_result& lpres,
00319                             const BCP_vec<BCP_var*>& vars,
00320                             const BCP_vec<BCP_cut*>& cuts,
00321                             const BCP_lp_var_pool& local_var_pool,
00322                             const BCP_lp_cut_pool& local_cut_pool,
00323                             BCP_vec<BCP_lp_branching_object*>& cands,
00324                             bool force_branch)
00325 {
00326   if (generated_vars > 0) {
00327     return BCP_DoNotBranch;
00328   }
00329 
00330   if (lpres.objval() > upper_bound() - 1e-6) {
00331     return BCP_DoNotBranch_Fathomed;
00332   }
00333 
00334   int i, j;
00335 
00336   const int dummyStart = par.entry(MCF3_par::AddDummySourceSinkArcs) ?
00337     data.numarcs - data.numcommodities : data.numarcs;
00338                 
00339   // find a few fractional original variables and do strong branching on them
00340   commodities_with_candidate.clear();
00341   for (i = data.numcommodities-1; i >= 0; --i) {
00342     new_branch[i].clear();
00343     std::map<int,double>& f = flows[i];
00344     int most_frac_ind = -1;
00345     double most_frac_val = 0.5-1e-6;
00346     double frac_val = 0.0;
00347     for (std::map<int,double>::iterator fi=f.begin(); fi != f.end(); ++fi){
00348       if (fi->first >= dummyStart)
00349         continue;
00350       const double frac = fabs(fi->second - floor(fi->second) - 0.5);
00351       if (frac < most_frac_val) {
00352         most_frac_ind = fi->first;
00353         most_frac_val = frac;
00354         frac_val = fi->second;
00355       }
00356     }
00357     if (most_frac_ind >= 0) {
00358       size_t pos;
00359       BCP_vec<int> ivp;
00360       BCP_vec<double> ivb;
00361       int lb = data.arcs[most_frac_ind].lb;
00362       int ub = data.arcs[most_frac_ind].ub;
00363       for (j = branch_history[i].size() - 1; j >= 0; --j) {
00364         // To correctly set lb/ub we need to check whether we have
00365         // already branched on this arc
00366         const MCF3_branch_decision& h = branch_history[i][j];
00367         if (h.arc_index == most_frac_ind) {
00368           lb = CoinMax(lb, h.lb);
00369           ub = CoinMin(ub, h.ub);
00370         }
00371       }
00372       const int mid = static_cast<int>(floor(frac_val));
00373       // Now the implied changes
00374       BCP_vec<int> child0_pos, child1_pos;
00375       BCP_vec<double> child0_bd, child1_bd;
00376 
00377       MCF3_branch_decision child0(most_frac_ind, lb, mid);
00378       new_branch[i].push_back(child0);
00379       MCF3_adjust_bounds(vars, i, child0, child0_pos, child0_bd);
00380 
00381       MCF3_branch_decision child1(most_frac_ind, mid+1, ub);
00382       new_branch[i].push_back(child1);
00383       MCF3_adjust_bounds(vars, i, child1, child1_pos, child1_bd);
00384 
00385       // child0_pos and child1_pos are going to be disjoint, thus
00386       // creating a merged list (need not be ordered!) is easy
00387       ivp.append(child0_pos);
00388       ivp.append(child1_pos);
00389       ivb.append(child0_bd);
00390       for (pos = 0; pos < child1_pos.size(); ++pos) {
00391         ivb.push_back(vars[child1_pos[pos]]->lb());
00392         ivb.push_back(vars[child1_pos[pos]]->ub());
00393       }
00394       for (pos = 0; pos < child0_pos.size(); ++pos) {
00395         ivb.push_back(vars[child0_pos[pos]]->lb());
00396         ivb.push_back(vars[child0_pos[pos]]->ub());
00397       }
00398       ivb.append(child1_bd);
00399       cands.push_back(new BCP_lp_branching_object(2, // num of children
00400                                                   NULL,
00401                                                   NULL, // no new cuts
00402                                                   NULL,NULL,NULL,NULL,
00403                                                   &ivp,NULL,&ivb,NULL));
00404       commodities_with_candidate.push_back(i);
00405     }
00406   }
00407   return BCP_DoBranch;
00408 }
00409 
00410 /*---------------------------------------------------------------------------*/
00411 
00412 void MCF3_lp::
00413 set_user_data_for_children(BCP_presolved_lp_brobj* best, 
00414                            const int selected)
00415 {
00416   const int commodity_branched_on = commodities_with_candidate[selected];
00417   BCP_vec<BCP_user_data*>& ud = best->user_data();
00418   MCF3_user* node_user = dynamic_cast<MCF3_user*>(get_user_data());
00419   MCF3_user* u;
00420 
00421   u = new MCF3_user(*node_user);
00422   u->branch_history[commodity_branched_on].
00423     push_back(new_branch[commodity_branched_on][0]);
00424   ud[0] = u;
00425     
00426   u = new MCF3_user(*node_user);
00427   u->branch_history[commodity_branched_on].
00428     push_back(new_branch[commodity_branched_on][1]);
00429   ud[1] = u;
00430 }
00431 
00432 /*===========================================================================*/
00433 
00434 void MCF3_lp::
00435 unpack_module_data(BCP_buffer& buf)
00436 {
00437   par.unpack(buf);
00438   data.unpack(buf);
00439 
00440   // This is the place where we can preallocate some data structures
00441   flows = new std::map<int,double>[data.numcommodities];
00442   new_branch = new std::vector<MCF3_branch_decision>[data.numcommodities];
00443 
00444   // Create the LP that will be used to generate columns
00445   cg_lp = new OsiClpSolverInterface();
00446 
00447   const int numCols = data.numarcs;
00448   const int numRows = data.numnodes;
00449   const int numNz = 2*numCols;
00450 
00451   double *clb = new double[numCols];
00452   double *cub = new double[numCols];
00453   double *obj = new double[numCols];
00454   double *rlb = new double[numRows];
00455   double *rub = new double[numRows];
00456   CoinBigIndex *start = new int[numCols+1];
00457   int *index = new int[numNz];
00458   double *value = new double[numNz];
00459 
00460   // all these will be properly set for the search tree node in the
00461   // initialize_new_search_tree_node method
00462   CoinZeroN(obj, numCols);
00463   CoinZeroN(clb, numCols);
00464   for (int i = numCols - 1; i >= 0; --i) {
00465     cub[i] = data.arcs[i].ub;
00466   }
00467   // and these will be properly set for the subproblem in the
00468   // generate_vars_in_lp method
00469   CoinZeroN(rlb, numRows);
00470   CoinZeroN(rub, numRows);
00471 
00472   for (int i = 0; i < data.numarcs; ++i) {
00473     start[i] = 2*i;
00474     index[2*i] = data.arcs[i].tail;
00475     index[2*i+1] = data.arcs[i].head;
00476     value[2*i] = -1;
00477     value[2*i+1] = 1;
00478   }
00479   start[numCols] = 2*numCols;
00480 
00481   cg_lp->loadProblem(numCols, numRows, start, index, value,
00482                      clb, cub, obj, rlb, rub);
00483 
00484   delete[] value;
00485   delete[] index;
00486   delete[] start;
00487   delete[] rub;
00488   delete[] rlb;
00489   delete[] obj;
00490   delete[] cub;
00491   delete[] clb;
00492 }                       

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