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);
00043 var_new_bd.push_back(0.0);
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);
00077 var_new_bd.push_back(0.0);
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
00095
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
00105
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
00120
00121
00122
00123
00124
00125
00126 MCF2_adjust_bounds(vars, branch_history, var_changed_pos, var_new_bd);
00127
00128
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
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
00151
00152
00153
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
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
00210
00211
00212
00213
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
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
00244
00245
00246
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
00274
00275
00276
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
00287
00288
00289
00290
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
00325
00326
00327
00328 BCP_col* col = new BCP_col(v->flow, v->weight, v->lb(), v->ub());
00329 col->insert(data.numarcs + v->commodity, 1.0);
00330 cols.push_back(col);
00331
00332
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
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 size_t pos;
00385 BCP_vec<BCP_var*> new_vars;
00386 BCP_vec<int> fvp;
00387 BCP_vec<double> fvb;
00388 BCP_vec<int> ivp;
00389 BCP_vec<double> ivb;
00390 int lb = data.arcs[most_frac_ind].lb;
00391 int ub = data.arcs[most_frac_ind].ub;
00392 for (j = branch_history[i].size() - 1; j >= 0; --j) {
00393
00394
00395 const MCF2_branch_decision& h = branch_history[i][j];
00396 if (h.arc_index == most_frac_ind) {
00397 lb = CoinMax(lb, h.lb);
00398 ub = CoinMin(ub, h.ub);
00399 }
00400 }
00401 const int mid = static_cast<int>(floor(frac_val));
00402 new_vars.push_back(new MCF2_branching_var(i, most_frac_ind,
00403 lb, mid, mid+1, ub));
00404 fvp.push_back(-1);
00405 fvb.push_back(0.0);
00406 fvb.push_back(0.0);
00407 fvb.push_back(1.0);
00408 fvb.push_back(1.0);
00409
00410 BCP_vec<int> child0_pos, child1_pos;
00411 BCP_vec<double> child0_bd, child1_bd;
00412
00413 MCF2_branch_decision child0(most_frac_ind, lb, mid);
00414 MCF2_adjust_bounds(vars, i, child0, child0_pos, child0_bd);
00415
00416 MCF2_branch_decision child1(most_frac_ind, mid+1, ub);
00417 MCF2_adjust_bounds(vars, i, child1, child1_pos, child1_bd);
00418
00419
00420
00421 ivp.append(child0_pos);
00422 ivp.append(child1_pos);
00423 ivb.append(child0_bd);
00424 for (pos = 0; pos < child1_pos.size(); ++pos) {
00425 ivb.push_back(vars[child1_pos[pos]]->lb());
00426 ivb.push_back(vars[child1_pos[pos]]->ub());
00427 }
00428 for (pos = 0; pos < child0_pos.size(); ++pos) {
00429 ivb.push_back(vars[child0_pos[pos]]->lb());
00430 ivb.push_back(vars[child0_pos[pos]]->ub());
00431 }
00432 ivb.append(child1_bd);
00433 cands.push_back(new BCP_lp_branching_object(2,
00434 &new_vars,
00435 NULL,
00436 &fvp,NULL,&fvb,NULL,
00437 &ivp,NULL,&ivb,NULL));
00438 }
00439 }
00440 return BCP_DoBranch;
00441 }
00442
00443
00444
00445 void MCF2_lp::
00446 unpack_module_data(BCP_buffer& buf)
00447 {
00448 par.unpack(buf);
00449 data.unpack(buf);
00450
00451
00452 branch_history = new std::vector<MCF2_branch_decision>[data.numcommodities];
00453 flows = new std::map<int,double>[data.numcommodities];
00454
00455
00456 cg_lp = new OsiClpSolverInterface();
00457
00458 const int numCols = data.numarcs;
00459 const int numRows = data.numnodes;
00460 const int numNz = 2*numCols;
00461
00462 double *clb = new double[numCols];
00463 double *cub = new double[numCols];
00464 double *obj = new double[numCols];
00465 double *rlb = new double[numRows];
00466 double *rub = new double[numRows];
00467 CoinBigIndex *start = new int[numCols+1];
00468 int *index = new int[numNz];
00469 double *value = new double[numNz];
00470
00471
00472
00473 CoinZeroN(obj, numCols);
00474 CoinZeroN(clb, numCols);
00475 for (int i = numCols - 1; i >= 0; --i) {
00476 cub[i] = data.arcs[i].ub;
00477 }
00478
00479
00480 CoinZeroN(rlb, numRows);
00481 CoinZeroN(rub, numRows);
00482
00483 for (int i = 0; i < data.numarcs; ++i) {
00484 start[i] = 2*i;
00485 index[2*i] = data.arcs[i].tail;
00486 index[2*i+1] = data.arcs[i].head;
00487 value[2*i] = -1;
00488 value[2*i+1] = 1;
00489 }
00490 start[numCols] = 2*numCols;
00491
00492 cg_lp->loadProblem(numCols, numRows, start, index, value,
00493 clb, cub, obj, rlb, rub);
00494
00495 delete[] value;
00496 delete[] index;
00497 delete[] start;
00498 delete[] rub;
00499 delete[] rlb;
00500 delete[] obj;
00501 delete[] cub;
00502 delete[] clb;
00503 }