/home/coin/SVN-release/OS-2.4.2/Bcp/examples/MCF-2/LP/MCF2_lp.cpp

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

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