00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include <stdio.h>
00055 #include <string.h>
00056 #include <map>
00057 #include <string>
00058 #include <utility>
00059 #include <vector>
00060
00061 #include "base/commandlineflags.h"
00062 #include "base/commandlineflags.h"
00063 #include "base/logging.h"
00064 #include "base/macros.h"
00065 #include "base/stringprintf.h"
00066 #include "linear_solver/linear_solver.h"
00067
00068 using std::string;
00069
00070 DEFINE_bool(colgen_verbose, false, "print verbosely");
00071 DEFINE_bool(colgen_complete, false, "generate all columns initially");
00072 DEFINE_int32(colgen_max_iterations, 500, "max iterations");
00073 DEFINE_string(colgen_solver, "clp", "solver - glpk or clp (default)");
00074 DEFINE_int32(colgen_instance, -1, "Which instance to solve (0 - 9)");
00075
00076 namespace operations_research {
00077
00078 struct Instance {
00079 int max_boxes;
00080 int width;
00081 int height;
00082 const char* grid;
00083 };
00084
00085 Instance kInstances[] = {
00086 { 4, 22, 6,
00087 "..@@@@@..............."
00088 "..@@@@@@........@@@..."
00089 ".....@@@@@......@@@..."
00090 ".......@@@@@@@@@@@@..."
00091 ".........@@@@@........"
00092 ".........@@@@@........"
00093 },
00094 { 3, 13, 10,
00095 "............."
00096 "............."
00097 "............."
00098 "...@@@@......"
00099 "...@@@@......"
00100 "...@@@@......"
00101 ".......@@@..."
00102 ".......@@@..."
00103 ".......@@@..."
00104 "............."
00105 },
00106 { 4, 13, 9,
00107 "............."
00108 "..@.@.@......"
00109 "...@.@.@....."
00110 "..@.@.@......"
00111 "..@.@.@......"
00112 "...@.@.@....."
00113 "....@.@......"
00114 "..........@@@"
00115 "..........@@@"
00116 },
00117 { 4, 13, 9,
00118 ".........@..."
00119 ".........@..."
00120 "@@@@@@@@@@..."
00121 "..@......@..."
00122 "..@......@..."
00123 "..@......@..."
00124 "..@@@@@@@@@@@"
00125 "..@.........."
00126 "..@.........."
00127 },
00128 { 7, 25, 14,
00129 "........................."
00130 "..@@@@@@@@@@@@@@@@@@@@..."
00131 "..@@@@@@@@@@@@@@@@@@@@..."
00132 "..@@.................@..."
00133 "..@@.................@..."
00134 "..@@.......@@@.......@.@."
00135 "..@@.......@@@.......@..."
00136 "..@@...@@@@@@@@@@@@@@@..."
00137 "..@@...@@@@@@@@@@@@@@@..."
00138 "..@@.......@@@.......@..."
00139 "..@@.......@@@.......@..."
00140 "..@@.................@..."
00141 "..@@.................@..."
00142 "........................."
00143 },
00144 { 6, 25, 16,
00145 "........................."
00146 "......@@@@@@@@@@@@@......"
00147 "........................."
00148 ".....@..........@........"
00149 ".....@..........@........"
00150 ".....@......@............"
00151 ".....@......@.@@@@@@@...."
00152 ".....@......@............"
00153 ".....@......@.@@@@@@@...."
00154 ".....@......@............"
00155 "....@@@@....@............"
00156 "....@@@@....@............"
00157 "..@@@@@@....@............"
00158 "..@@@.......@............"
00159 "..@@@...................."
00160 "..@@@@@@@@@@@@@@@@@@@@@@@"
00161 },
00162 { 5, 40, 18,
00163 "........................................"
00164 "........................................"
00165 "...@@@@@@..............................."
00166 "...@@@@@@..............................."
00167 "...@@@@@@..............................."
00168 "...@@@@@@.........@@@@@@@@@@............"
00169 "...@@@@@@.........@@@@@@@@@@............"
00170 "..................@@@@@@@@@@............"
00171 "..................@@@@@@@@@@............"
00172 ".............@@@@@@@@@@@@@@@............"
00173 ".............@@@@@@@@@@@@@@@............"
00174 "........@@@@@@@@@@@@...................."
00175 "........@@@@@@@@@@@@...................."
00176 "........@@@@@@.........................."
00177 "........@@@@@@.........................."
00178 "........................................"
00179 "........................................"
00180 "........................................"
00181 },
00182 { 8, 40, 18,
00183 "........................................"
00184 "..@@.@.@.@.............................."
00185 "..@@.@.@.@...............@.............."
00186 "..@@.@.@.@............@................."
00187 "..@@.@.@.@.............................."
00188 "..@@.@.@.@.................@............"
00189 "..@@.@..................@..............."
00190 "..@@.@.................................."
00191 "..@@.@.................................."
00192 "..@@.@................@@@@.............."
00193 "..@@.@..............@@@@@@@@............"
00194 "..@@.@.................................."
00195 "..@@.@..............@@@@@@@@............"
00196 "..@@.@.................................."
00197 "..@@.@................@@@@.............."
00198 "..@@.@.................................."
00199 "..@@.@.................................."
00200 "........................................"
00201 },
00202 { 10, 40, 19,
00203 "@@@@@..................................."
00204 "@@@@@..................................."
00205 "@@@@@..................................."
00206 "@@@@@..................................."
00207 "@@@@@..................................."
00208 "@@@@@...........@@@@@@@@@@@............."
00209 "@@@@@...........@@@@@@@@@@@............."
00210 "....................@@@@................"
00211 "....................@@@@................"
00212 "....................@@@@................"
00213 "....................@@@@................"
00214 "....................@@@@................"
00215 "...............@@@@@@@@@@@@@@..........."
00216 "...............@@@@@@@@@@@@@@..........."
00217 ".......@@@@@@@@@@@@@@@@@@@@@@..........."
00218 ".......@@@@@@@@@........................"
00219 "........................................"
00220 "........................................"
00221 "........................................"
00222 },
00223 { 10, 40, 25,
00224 "...................@...................."
00225 "...............@@@@@@@@@................"
00226 "............@@@.........@@@............."
00227 "...........@...............@............"
00228 "..........@.................@..........."
00229 ".........@...................@.........."
00230 ".........@...................@.........."
00231 ".........@.....@@......@@....@.........."
00232 "........@.....@@@@....@@@@....@........."
00233 "........@.....................@........."
00234 "........@.....................@........."
00235 "........@..........@@.........@........."
00236 ".......@@..........@@.........@@........"
00237 "........@.....................@........."
00238 "........@.....................@........."
00239 "........@......@@@@@@@@@......@........."
00240 "........@......@@@@@@@@@......@........."
00241 ".........@...................@.........."
00242 ".........@...................@.........."
00243 ".........@...................@.........."
00244 "..........@.................@..........."
00245 "...........@...............@............"
00246 "............@@@.........@@@............."
00247 "...............@@@@@@@@@................"
00248 "...................@...................."
00249 }
00250 };
00251
00252 const int kInstanceCount = 10;
00253
00254
00255
00256 class Box {
00257 public:
00258 static const int kAreaCost = 1;
00259 static const int kFixedCost = 10;
00260
00261 Box() {}
00262 Box(int x_min, int x_max, int y_min, int y_max)
00263 : x_min_(x_min), x_max_(x_max), y_min_(y_min), y_max_(y_max) {
00264 CHECK_GE(x_max, x_min);
00265 CHECK_GE(y_max, y_min);
00266 }
00267
00268 int x_min() const { return x_min_; }
00269 int x_max() const { return x_max_; }
00270 int y_min() const { return y_min_; }
00271 int y_max() const { return y_max_; }
00272
00273
00274 int Compare(const Box& box) const {
00275 int c;
00276 if ((c = (x_min() - box.x_min())) != 0) return c;
00277 if ((c = (x_max() - box.x_max())) != 0) return c;
00278 if ((c = (y_min() - box.y_min())) != 0) return c;
00279 return y_max() - box.y_max();
00280 }
00281
00282 bool Contains(int x, int y) const {
00283 return x >= x_min() && x <= x_max() && y >= y_min() && y <= y_max();
00284 }
00285
00286 int Cost() const {
00287 return kAreaCost * (x_max() - x_min() + 1) * (y_max() - y_min() + 1)
00288 + kFixedCost;
00289 }
00290
00291 string DebugString() const {
00292 return StringPrintf("[%d,%dx%d,%d]c%d",
00293 x_min(), y_min(), x_max(), y_max(), Cost());
00294 }
00295
00296 private:
00297 int x_min_;
00298 int x_max_;
00299 int y_min_;
00300 int y_max_;
00301 };
00302
00303 struct BoxLessThan {
00304 bool operator()(const Box& b1, const Box& b2) const {
00305 return b1.Compare(b2) < 0;
00306 }
00307 };
00308
00309
00310
00311 class CoveringProblem {
00312 public:
00313
00314
00315
00316 CoveringProblem(MPSolver* const solver, const Instance& instance)
00317 : solver_(solver),
00318 max_boxes_(instance.max_boxes),
00319 width_(instance.width),
00320 height_(instance.height),
00321 grid_(instance.grid) {}
00322
00323
00324
00325 bool Init() {
00326
00327 int size = strlen(grid_);
00328 if (size != area()) {
00329 return false;
00330 }
00331 for (int i = 0; i < size; ++i) {
00332 char c = grid_[i];
00333 if ((c != '@') && (c != '.')) return false;
00334 }
00335
00336 AddCellConstraints();
00337 AddMaxBoxesConstraint();
00338 if (!FLAGS_colgen_complete) {
00339 AddBox(Box(0, width()-1, 0, height()-1));
00340 } else {
00341
00342
00343 for (int y_min = 0; y_min < height(); ++y_min) {
00344 for (int y_max = y_min; y_max < height(); ++y_max) {
00345 for (int x_min = 0; x_min < width(); ++x_min) {
00346 for (int x_max = x_min; x_max < width(); ++x_max) {
00347 AddBox(Box(x_min, x_max, y_min, y_max));
00348 }
00349 }
00350 }
00351 }
00352 }
00353 return true;
00354 }
00355
00356 int width() const { return width_; }
00357 int height() const { return height_; }
00358 int area() const { return width() * height(); }
00359 int max_boxes() const { return max_boxes_; }
00360
00361 bool IsCellOccupied(int x, int y) const { return grid_[index(x, y)] == '@'; }
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 double GetOptimalBox(Box* const target) {
00384
00385 const double kCostChangeThreshold = -.01;
00386
00387
00388
00389 std::vector<double> upper_left_sums(area());
00390 ComputeUpperLeftSums(&upper_left_sums);
00391
00392 const double max_boxes_dual = max_boxes_constraint_->dual_value();
00393 double best_reduced_cost = kCostChangeThreshold;
00394 Box best_box;
00395 for (int y_min = 0; y_min < height(); ++y_min) {
00396 for (int y_max = y_min; y_max < height(); ++y_max) {
00397 for (int x_min = 0; x_min < width(); ++x_min) {
00398 for (int x_max = x_min; x_max < width(); ++x_max) {
00399 Box box(x_min, x_max, y_min, y_max);
00400 const double cell_coverage_dual =
00401 + zero_access(upper_left_sums, x_max, y_max)
00402 - zero_access(upper_left_sums, x_max, y_min - 1)
00403 - zero_access(upper_left_sums, x_min - 1, y_max)
00404 + zero_access(upper_left_sums, x_min - 1, y_min - 1);
00405
00406
00407
00408
00409 const double reduced_cost = box.Cost()
00410 - (cell_coverage_dual + max_boxes_dual);
00411
00412 if (reduced_cost < best_reduced_cost) {
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422 if (boxes_.find(box) == boxes_.end()) {
00423 best_reduced_cost = reduced_cost;
00424 best_box = box;
00425 }
00426 }
00427 }
00428 }
00429 }
00430 }
00431
00432 if (best_reduced_cost < kCostChangeThreshold) {
00433 if (target) {
00434 *target = best_box;
00435 }
00436 return best_reduced_cost;
00437 }
00438 return 0;
00439 }
00440
00441
00442
00443 MPVariable* AddBox(const Box& box) {
00444 CHECK(boxes_.find(box) == boxes_.end());
00445 MPVariable* const var = solver_->MakeNumVar(0., 1., box.DebugString());
00446 solver_->SetObjectiveCoefficient(var, box.Cost());
00447 max_boxes_constraint_->SetCoefficient(var, 1.0);
00448 for (int y = box.y_min(); y <= box.y_max(); ++y) {
00449 for (int x = box.x_min(); x <= box.x_max(); ++x) {
00450 cell(x, y)->SetCoefficient(var, 1.0);
00451 }
00452 }
00453 boxes_[box] = var;
00454 return var;
00455 }
00456
00457 string PrintGrid() const {
00458 string output = StringPrintf("width = %d, height = %d, max_boxes = %d\n",
00459 width_,
00460 height_,
00461 max_boxes_);
00462 for (int y = 0; y < height_; ++y) {
00463 StringAppendF(&output,
00464 "%s\n",
00465 string(grid_ + width_ * y, width_).c_str());
00466 }
00467 return output;
00468 }
00469
00470
00471
00472
00473
00474 string PrintCovering() const {
00475 static const double kTolerance = 1e-5;
00476 string output = StringPrintf("cost = %lf\n", solver_->objective_value());
00477 scoped_array<char> display(new char[(width_ + 1) * height_ + 1]);
00478 for (int y = 0; y < height_; ++y) {
00479 memcpy(display.get() + y * (width_ + 1),
00480 grid_ + width_ * y,
00481 width_);
00482 display[y * (width_ + 1) + width_] = '\n';
00483 }
00484 display[height_ * (width_ + 1)] = '\0';
00485 int active_box_index = 0;
00486 for (BoxTable::const_iterator i = boxes_.begin(); i != boxes_.end(); ++i) {
00487 const double value = i->second->solution_value();
00488 if (value > kTolerance) {
00489 const char box_character =
00490 (i->second->solution_value() >= (1. - kTolerance) ? 'A' : 'a') +
00491 active_box_index++;
00492 StringAppendF(&output, "%c: box %s with value %lf\n",
00493 box_character,
00494 i->first.DebugString().c_str(),
00495 value);
00496 const Box& box = i->first;
00497 for (int x = box.x_min(); x <= box.x_max(); ++x) {
00498 for (int y = box.y_min(); y <= box.y_max(); ++y) {
00499 display[x + y * (width_ + 1)] = box_character;
00500 }
00501 }
00502 }
00503 }
00504 StringAppendF(&output, "%s", display.get());
00505 return output;
00506 }
00507
00508 protected:
00509 int index(int x, int y) const { return width_ * y + x; }
00510 MPConstraint* cell(int x, int y) { return cells_[index(x, y)]; }
00511 const MPConstraint* cell(int x, int y) const { return cells_[index(x, y)]; }
00512
00513
00514
00515 void AddCellConstraints() {
00516 cells_.resize(area());
00517 for (int y = 0; y < height(); ++y) {
00518 for (int x = 0; x < width(); ++x) {
00519 cells_[index(x, y)] = solver_->MakeRowConstraint
00520 (IsCellOccupied(x, y) ? 1. : 0., 1.);
00521 }
00522 }
00523 }
00524
00525
00526 void AddMaxBoxesConstraint() {
00527 max_boxes_constraint_ = solver_->MakeRowConstraint(0., max_boxes());
00528 }
00529
00530
00531 double zero_access(const std::vector<double>& array, int x, int y) const {
00532 if (x < 0 || y < 0) {
00533 return 0;
00534 }
00535 return array[index(x, y)];
00536 }
00537
00538
00539
00540 void ComputeUpperLeftSums(std::vector<double>* upper_left_sums) const {
00541 for (int y = 0; y < height(); ++y) {
00542 for (int x = 0; x < width(); ++x) {
00543 upper_left_sums->operator[](index(x, y)) =
00544 cell(x, y)->dual_value()
00545 + zero_access(*upper_left_sums, x - 1, y)
00546 + zero_access(*upper_left_sums, x, y - 1)
00547 - zero_access(*upper_left_sums, x - 1, y - 1);
00548 }
00549 }
00550 }
00551
00552 typedef std::map<Box, MPVariable*, BoxLessThan> BoxTable;
00553 MPSolver* const solver_;
00554 const int max_boxes_;
00555 const int width_;
00556 const int height_;
00557 const char* const grid_;
00558 std::vector<MPConstraint*> cells_;
00559 BoxTable boxes_;
00560 MPConstraint* max_boxes_constraint_;
00561 };
00562
00563
00564
00565
00566
00567 void SolveInstance(const Instance& instance,
00568 MPSolver::OptimizationProblemType solver_type) {
00569
00570 MPSolver solver("ColumnGeneration", solver_type);
00571 solver.SuppressOutput();
00572 solver.SetMinimization();
00573
00574
00575 CoveringProblem problem(&solver, instance);
00576 CHECK(problem.Init());
00577 LOG(INFO) << "Initial problem:\n" << problem.PrintGrid();
00578
00579 int step_number = 0;
00580 while (step_number < FLAGS_colgen_max_iterations) {
00581 if (FLAGS_colgen_verbose) {
00582 LOG(INFO) << "Step number " << step_number;
00583 }
00584
00585
00586 CHECK_EQ(MPSolver::OPTIMAL, solver.Solve());
00587 if (FLAGS_colgen_verbose) {
00588 LOG(INFO) << problem.PrintCovering();
00589 }
00590
00591
00592 Box box;
00593 const double reduced_cost = problem.GetOptimalBox(&box);
00594 if (reduced_cost >= 0) {
00595 break;
00596 }
00597
00598
00599 if (FLAGS_colgen_verbose) {
00600 LOG(INFO) << "Adding " << box.DebugString()
00601 << ", reduced_cost =" << reduced_cost;
00602 }
00603 problem.AddBox(box);
00604
00605 ++step_number;
00606 }
00607
00608 if (step_number >= FLAGS_colgen_max_iterations) {
00609
00610 CHECK_EQ(MPSolver::OPTIMAL, solver.Solve());
00611 }
00612
00613 LOG(INFO) << step_number << " columns added";
00614 LOG(INFO) << "Final coverage: " << problem.PrintCovering();
00615 }
00616 }
00617
00618 int main(int argc, char** argv) {
00619 string usage = "column_generation\n";
00620 usage += " --colgen_verbose print verbosely\n";
00621 usage += " --colgen_max_iterations <n> max columns to generate\n";
00622 usage += " --colgen_complete generate all columns at start\n";
00623
00624
00625 operations_research::MPSolver::OptimizationProblemType solver_type;
00626 bool found = false;
00627 #if defined(USE_CLP)
00628 if (FLAGS_colgen_solver == "clp") {
00629 solver_type = operations_research::MPSolver::CLP_LINEAR_PROGRAMMING;
00630 found = true;
00631 }
00632 #endif // USE_CLP
00633 #if defined(USE_GLPK)
00634 if (FLAGS_colgen_solver == "glpk") {
00635 solver_type = operations_research::MPSolver::GLPK_LINEAR_PROGRAMMING;
00636 found = true;
00637 }
00638 #endif // USE_GLPK
00639 if (!found) {
00640 LOG(ERROR) << "Unknown solver " << FLAGS_colgen_solver;
00641 return 1;
00642 }
00643
00644 if (FLAGS_colgen_instance == -1) {
00645 for (int i = 0; i < operations_research::kInstanceCount; ++i) {
00646 const operations_research::Instance& instance =
00647 operations_research::kInstances[i];
00648 operations_research::SolveInstance(instance, solver_type);
00649 }
00650 } else {
00651 CHECK_GE(FLAGS_colgen_instance, 0);
00652 CHECK_LT(FLAGS_colgen_instance, operations_research::kInstanceCount);
00653 const operations_research::Instance& instance =
00654 operations_research::kInstances[FLAGS_colgen_instance];
00655 operations_research::SolveInstance(instance, solver_type);
00656 }
00657 return 0;
00658 }