00001 #include "CoinHelperFunctions.hpp"
00002 #include "OsiClpSolverInterface.hpp"
00003 #include "MCF1_lp.hpp"
00004
00005
00006
00007 OsiSolverInterface* MCF1_lp::initialize_solver_interface()
00008 {
00009 return new OsiClpSolverInterface;
00010 }
00011
00012
00013
00014 void
00015 MCF1_adjust_bounds(const BCP_vec<BCP_var*>& vars,
00016 std::vector<MCF1_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 MCF1_var* v = dynamic_cast<MCF1_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<MCF1_branch_decision>& hist = history[v->commodity];
00029 for (int j = hist.size()-1; !violated && j >= 0; --j) {
00030 const MCF1_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 MCF1_adjust_bounds(const BCP_vec<BCP_var*>& vars,
00052 const int commodity,
00053 const MCF1_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 MCF1_var* v = dynamic_cast<MCF1_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 MCF1_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 MCF1_branching_var* v = dynamic_cast<MCF1_branching_var*>(vars[i]);
00102 if (v) {
00103 if (v->ub() == 0.0) {
00104
00105
00106 branch_history[v->commodity].
00107 push_back(MCF1_branch_decision(v->arc_index,
00108 v->lb_child0,
00109 v->ub_child0));
00110 } else {
00111 branch_history[v->commodity].
00112 push_back(MCF1_branch_decision(v->arc_index,
00113 v->lb_child1,
00114 v->ub_child1));
00115 }
00116 }
00117 }
00118
00119 MCF1_adjust_bounds(vars, branch_history, var_changed_pos, var_new_bd);
00120
00121
00122 purge_ptr_vector(gen_vars);
00123 }
00124
00125
00126
00127 BCP_solution* MCF1_lp::
00128 test_feasibility(const BCP_lp_result& lpres,
00129 const BCP_vec<BCP_var*>& vars,
00130 const BCP_vec<BCP_cut*>& cuts)
00131 {
00132 #if 0
00133 static int cnt = 0;
00134 char name[100];
00135 sprintf(name, "currlp-%i", cnt++);
00136
00137 getLpProblemPointer()->lp_solver->writeMps(name, "mps");
00138 printf("Current LP written in file %s.mps\n", name);
00139 getLpProblemPointer()->lp_solver->writeLp(name, "lp");
00140 printf("Current LP written in file %s.lp\n", name);
00141 #endif
00142
00143
00144
00145
00146
00147
00148 int i, j;
00149 for (i = data.numcommodities-1; i >= 0; --i) {
00150 flows[i].clear();
00151 }
00152
00153 const double* x = lpres.x();
00154 for (i = vars.size()-1; i >= 0; --i) {
00155 MCF1_var* v = dynamic_cast<MCF1_var*>(vars[i]);
00156 if (!v) continue;
00157 std::map<int,double>& f = flows[v->commodity];
00158 const int vsize = v->flow.getNumElements();
00159 const int* vind = v->flow.getIndices();
00160 const double* vval = v->flow.getElements();
00161 for (j = 0; j < vsize; ++j) {
00162 f[vind[j]] += vval[j]*x[i];
00163 }
00164 }
00165 for (i = data.numcommodities-1; i >= 0; --i) {
00166 std::map<int,double>& f = flows[i];
00167 for (std::map<int,double>::iterator fi=f.begin(); fi != f.end(); ++fi) {
00168 const double frac = fabs(fi->second - floor(fi->second) - 0.5);
00169 if (frac < 0.5-1e-6)
00170 return NULL;
00171 }
00172 }
00173
00174
00175 BCP_solution_generic* gsol = new BCP_solution_generic;
00176 for (i = data.numcommodities-1; i >= 0; --i) {
00177 std::map<int,double>& f = flows[i];
00178 CoinPackedVector flow(false);
00179 double weight = 0;
00180 for (std::map<int,double>::iterator fi=f.begin();
00181 fi != f.end();
00182 ++fi) {
00183 const double val = floor(fi->second + 0.5);
00184 if (val > 0) {
00185 flow.insert(fi->first, val);
00186 weight += val * data.arcs[fi->first].weight;
00187 }
00188 }
00189 gsol->add_entry(new MCF1_var(i, flow, weight), 1.0);
00190 }
00191 return gsol;
00192 }
00193
00194
00195
00196 double MCF1_lp::
00197 compute_lower_bound(const double old_lower_bound,
00198 const BCP_lp_result& lpres,
00199 const BCP_vec<BCP_var*>& vars,
00200 const BCP_vec<BCP_cut*>& cuts)
00201 {
00202
00203
00204
00205
00206
00207
00208 const int numarcs = data.numarcs;
00209 double* cost = new double[numarcs];
00210 const double* pi = lpres.pi();
00211 const double* nu = pi + numarcs;
00212
00213 int i, j;
00214
00215 for (i = numarcs-1; i >= 0; --i) {
00216 cost[i] = data.arcs[i].weight - pi[i];
00217 }
00218 cg_lp->setObjective(cost);
00219
00220
00221 int* ind = new int[numarcs];
00222 double* val = new double[numarcs];
00223 int cnt = 0;
00224
00225 for (i = data.numcommodities-1; i >= 0; --i) {
00226 const MCF1_data::commodity& comm = data.commodities[i];
00227 cg_lp->setRowBounds(comm.source, -comm.demand, -comm.demand);
00228 cg_lp->setRowBounds(comm.sink, comm.demand, comm.demand);
00229 const std::vector<MCF1_branch_decision>& hist = branch_history[i];
00230 for (j = hist.size() - 1; j >= 0; --j) {
00231 const MCF1_branch_decision& h = hist[j];
00232 cg_lp->setColBounds(h.arc_index, h.lb, h.ub);
00233 }
00234 cg_lp->initialSolve();
00235 if (cg_lp->isProvenOptimal() && cg_lp->getObjValue() < nu[i] - 1e-8) {
00236
00237
00238
00239
00240 const double* x = cg_lp->getColSolution();
00241 cnt = 0;
00242 double obj = 0.0;
00243 for (int j = 0; j < numarcs; ++j) {
00244 const double xval = floor(x[j] + 0.5);
00245 if (xval != 0.0) {
00246 ind[cnt] = j;
00247 val[cnt] = xval;
00248 ++cnt;
00249 obj += data.arcs[j].weight * xval;
00250 }
00251 }
00252 gen_vars.push_back(new MCF1_var(i, CoinPackedVector(cnt, ind, val,
00253 false), obj));
00254 }
00255 for (j = hist.size() - 1; j >= 0; --j) {
00256 const int ind = hist[j].arc_index;
00257 cg_lp->setColBounds(ind, data.arcs[ind].lb, data.arcs[ind].ub);
00258 }
00259 cg_lp->setRowBounds(comm.source, 0, 0);
00260 cg_lp->setRowBounds(comm.sink, 0, 0);
00261 }
00262 delete[] val;
00263 delete[] ind;
00264 delete[] cost;
00265
00266
00267
00268
00269
00270 double true_lower_bound = 0.0;
00271 generated_vars = (gen_vars.size() > 0);
00272
00273 if (generated_vars) {
00274 true_lower_bound = old_lower_bound;
00275 } else {
00276 true_lower_bound = lpres.objval();
00277 }
00278
00279
00280
00281
00282
00283
00284
00285 return true_lower_bound;
00286 }
00287
00288
00289
00290 void MCF1_lp::
00291 generate_vars_in_lp(const BCP_lp_result& lpres,
00292 const BCP_vec<BCP_var*>& vars,
00293 const BCP_vec<BCP_cut*>& cuts,
00294 const bool before_fathom,
00295 BCP_vec<BCP_var*>& new_vars,
00296 BCP_vec<BCP_col*>& new_cols)
00297 {
00298 new_vars.append(gen_vars);
00299 gen_vars.clear();
00300 }
00301
00302
00303
00304 void MCF1_lp::
00305 vars_to_cols(const BCP_vec<BCP_cut*>& cuts,
00306 BCP_vec<BCP_var*>& vars,
00307 BCP_vec<BCP_col*>& cols,
00308 const BCP_lp_result& lpres,
00309 BCP_object_origin origin, bool allow_multiple)
00310 {
00311 static const CoinPackedVector emptyVector(false);
00312 const int numvars = vars.size();
00313 for (int i = 0; i < numvars; ++i) {
00314 const MCF1_var* v =
00315 dynamic_cast<const MCF1_var*>(vars[i]);
00316 if (v) {
00317
00318
00319
00320
00321 BCP_col* col = new BCP_col(v->flow, v->weight, v->lb(), v->ub());
00322 col->insert(data.numarcs + v->commodity, 1.0);
00323 cols.push_back(col);
00324
00325
00326 continue;
00327 }
00328 const MCF1_branching_var* bv =
00329 dynamic_cast<const MCF1_branching_var*>(vars[i]);
00330 if (bv) {
00331 cols.push_back(new BCP_col(emptyVector, 0.0, bv->lb(), bv->ub()));
00332 }
00333 }
00334 }
00335
00336
00337
00338 BCP_branching_decision MCF1_lp::
00339 select_branching_candidates(const BCP_lp_result& lpres,
00340 const BCP_vec<BCP_var*>& vars,
00341 const BCP_vec<BCP_cut*>& cuts,
00342 const BCP_lp_var_pool& local_var_pool,
00343 const BCP_lp_cut_pool& local_cut_pool,
00344 BCP_vec<BCP_lp_branching_object*>& cands,
00345 bool force_branch)
00346 {
00347 if (generated_vars > 0) {
00348 return BCP_DoNotBranch;
00349 }
00350
00351 if (lpres.objval() > upper_bound() - 1e-6) {
00352 return BCP_DoNotBranch_Fathomed;
00353 }
00354
00355 int i, j;
00356
00357 const int dummyStart = par.entry(MCF1_par::AddDummySourceSinkArcs) ?
00358 data.numarcs - data.numcommodities : data.numarcs;
00359
00360
00361 for (i = data.numcommodities-1; i >= 0; --i) {
00362 std::map<int,double>& f = flows[i];
00363 int most_frac_ind = -1;
00364 double most_frac_val = 0.5-1e-6;
00365 double frac_val = 0.0;
00366 for (std::map<int,double>::iterator fi=f.begin(); fi != f.end(); ++fi){
00367 if (fi->first >= dummyStart)
00368 continue;
00369 const double frac = fabs(fi->second - floor(fi->second) - 0.5);
00370 if (frac < most_frac_val) {
00371 most_frac_ind = fi->first;
00372 most_frac_val = frac;
00373 frac_val = fi->second;
00374 }
00375 }
00376 if (most_frac_ind >= 0) {
00377 size_t pos;
00378 BCP_vec<BCP_var*> new_vars;
00379 BCP_vec<int> fvp;
00380 BCP_vec<double> fvb;
00381 int lb = data.arcs[most_frac_ind].lb;
00382 int ub = data.arcs[most_frac_ind].ub;
00383 for (j = branch_history[i].size() - 1; j >= 0; --j) {
00384
00385
00386 const MCF1_branch_decision& h = branch_history[i][j];
00387 if (h.arc_index == most_frac_ind) {
00388 lb = CoinMax(lb, h.lb);
00389 ub = CoinMin(ub, h.ub);
00390 }
00391 }
00392 const int mid = static_cast<int>(floor(frac_val));
00393 new_vars.push_back(new MCF1_branching_var(i, most_frac_ind,
00394 lb, mid, mid+1, ub));
00395
00396 BCP_vec<int> child0_pos, child1_pos;
00397 BCP_vec<double> child0_bd, child1_bd;
00398
00399 MCF1_branch_decision child0(most_frac_ind, lb, mid);
00400 MCF1_adjust_bounds(vars, i, child0, child0_pos, child0_bd);
00401
00402 MCF1_branch_decision child1(most_frac_ind, mid+1, ub);
00403 MCF1_adjust_bounds(vars, i, child1, child1_pos, child1_bd);
00404
00405
00406 fvp.push_back(-1);
00407 fvp.append(child0_pos);
00408 fvp.append(child1_pos);
00409 fvb.push_back(0.0);
00410 fvb.push_back(0.0);
00411 fvb.append(child0_bd);
00412 for (pos = 0; pos < child1_pos.size(); ++pos) {
00413 fvb.push_back(vars[child1_pos[pos]]->lb());
00414 fvb.push_back(vars[child1_pos[pos]]->ub());
00415 }
00416 fvb.push_back(1.0);
00417 fvb.push_back(1.0);
00418 for (pos = 0; pos < child0_pos.size(); ++pos) {
00419 fvb.push_back(vars[child0_pos[pos]]->lb());
00420 fvb.push_back(vars[child0_pos[pos]]->ub());
00421 }
00422 fvb.append(child1_bd);
00423 cands.push_back(new BCP_lp_branching_object(2,
00424 &new_vars,
00425 NULL,
00426 &fvp,NULL,&fvb,NULL,
00427 NULL,NULL,NULL,NULL));
00428 }
00429 }
00430 return BCP_DoBranch;
00431 }
00432
00433
00434
00435 void MCF1_lp::
00436 unpack_module_data(BCP_buffer& buf)
00437 {
00438 par.unpack(buf);
00439 data.unpack(buf);
00440
00441
00442 branch_history = new std::vector<MCF1_branch_decision>[data.numcommodities];
00443 flows = new std::map<int,double>[data.numcommodities];
00444
00445
00446 cg_lp = new OsiClpSolverInterface();
00447
00448 const int numCols = data.numarcs;
00449 const int numRows = data.numnodes;
00450 const int numNz = 2*numCols;
00451
00452 double *clb = new double[numCols];
00453 double *cub = new double[numCols];
00454 double *obj = new double[numCols];
00455 double *rlb = new double[numRows];
00456 double *rub = new double[numRows];
00457 CoinBigIndex *start = new int[numCols+1];
00458 int *index = new int[numNz];
00459 double *value = new double[numNz];
00460
00461
00462
00463 CoinZeroN(obj, numCols);
00464 CoinZeroN(clb, numCols);
00465 for (int i = numCols - 1; i >= 0; --i) {
00466 cub[i] = data.arcs[i].ub;
00467 }
00468
00469
00470 CoinZeroN(rlb, numRows);
00471 CoinZeroN(rub, numRows);
00472
00473 for (int i = 0; i < data.numarcs; ++i) {
00474 start[i] = 2*i;
00475 index[2*i] = data.arcs[i].tail;
00476 index[2*i+1] = data.arcs[i].head;
00477 value[2*i] = -1;
00478 value[2*i+1] = 1;
00479 }
00480 start[numCols] = 2*numCols;
00481
00482 cg_lp->loadProblem(numCols, numRows, start, index, value,
00483 clb, cub, obj, rlb, rub);
00484
00485 delete[] value;
00486 delete[] index;
00487 delete[] start;
00488 delete[] rub;
00489 delete[] rlb;
00490 delete[] obj;
00491 delete[] cub;
00492 delete[] clb;
00493 }