always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

reorganize all examples in the examples directory

Browse Source
This commit is contained in:
lperron@google.com
2012-03-28 14:23:23 +00:00
parent dd602fde3c
commit 64752231d6
786 changed files with 213 additions and 206 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
examples/cpp/README Normal file
View File

@@ -0,0 +1,55 @@
<summary> C++ examples to demonstrate usage of the different
Operations Research libraries. </summary>
- %Constraint Solver examples:
- cryptarithm.cc Demonstrates the use of basic modeling objects
(integer variables, arithmetic constraints and expressions,
simple search).
- golomb.cc Demonstrates how to handle objective functions and collect
solutions found during the search.
- magic_square.cc Shows how to use the automatic search to solve your
problem.
- costas_array.cc Solves the problem of Costas Array (a constrained
assignment problem used for radars) with two version. On version is
a feasibility version with hard constraints, the other version is
an optimization version with soft constraints and violation costs.
- jobshop.cc Demonstrates scheduling of jobs on different machines.
- jobshop_ls.cc Demonstrates scheduling of jobs on different machines with
a search using Local Search and Large Neighorhood Search.
- nqueens.cc Solves the n-queen problem. It also demonstrates how to break
symmetries during search.
- network_routing.cc Solves a multicommodity mono-routing
problem with capacity constraints and a max usage cost structure.
- sports_scheduling.cc Finds a soccer championship schedule. Its uses an
original approach where all constraints attached to either one team,
or one week are regrouped into one global 'AllowedAssignment' constraints.
- dobble_ls.cc Shows how to write Local Search operators and Local Search
filters in a context of an assignment/partitioning problem. It also
shows how to write a simple constraint.
- Routing examples:
- tsp.cc Travelling Salesman Problem.
- cvrptw.cc Capacitated Vehicle Routing Problem with Time Windows.
- pdptw.cc Pickup and Delivery Problem with Time Windows.
- Graph examples:
- flow_api.cc Demonstrates how to use Min-Cost Flow and Max-Flow api.
- linear_assignment_api.cc Demonstrates how to use the Linear Sum
Assignment solver.
- dimacs_assignment.cc Solves DIMACS challenge on assignment
problems.
- Linear and integer programming examples:
- linear_programming.cc Demonstrates how to use the linear solver
wrapper API to solve Linear Programming problems.
- integer_programming.cc Demonstrates how to use the linear solver
wrapper API to solve Integer Programming problems.
- linear_solver_protocol_buffers.cc Demonstrates how protocol
buffers can be used as input and output to the linear solver wrapper.
- strawberry_fields_with_column_generation.cc Complex example that
demonstrates how to use dynamic column generation to solve a 2D
covering problem.
- Utilities
- model_util.cc A utility to manipulate model files (.cp) dumped by the
solver.

476
examples/cpp/costas_array.cc Normal file
View File

@@ -0,0 +1,476 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Costas Array Problem
//
// Finds an NxN matrix of 0s and 1s, with only one 1 per row,
// and one 1 per column, such that all displacement vectors
// between each pair of 1s are distinct.
//
// This example contains two separate implementations. CostasHard()
// uses hard constraints, whereas CostasSoft() uses a minimizer to
// minimize the number of duplicates.
#include <time.h>
#include <set>
#include <utility>
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "constraint_solver/constraint_solver.h"
#include "constraint_solver/constraint_solveri.h"
#include "base/random.h"
DEFINE_int32(minsize, 0, "Minimum degree of Costas matrix.");
DEFINE_int32(maxsize, 0, "Maximum degree of Costas matrix.");
DEFINE_int32(freevar, 5, "Number of free variables.");
DEFINE_int32(freeorderedvar, 4, "Number of variables in ordered subset.");
DEFINE_int32(sublimit, 20, "Number of attempts per subtree.");
DEFINE_int32(timelimit, 120000, "Time limit for local search.");
DEFINE_bool(soft_constraints, false, "Use soft constraints.");
DEFINE_string(export_profile, "", "filename to save the profile overview");
namespace operations_research {
// Checks that all pairwise distances are unique and returns all violators
void CheckConstraintViolators(const std::vector<int64>& vars,
std::vector<int>* const violators) {
int dim = vars.size();
// Check that all indices are unique
for (int i = 0; i < dim; ++i) {
for (int k = i + 1; k < dim; ++k) {
if (vars[i] == vars[k]) {
violators->push_back(i);
violators->push_back(k);
}
}
}
// Check that all differences are unique for each level
for (int level = 1; level < dim; ++level) {
for (int i = 0; i < dim - level; ++i) {
const int difference = vars[i + level] - vars[i];
for (int k = i + 1; k < dim - level; ++k) {
if (difference == vars[k + level] - vars[k]) {
violators->push_back(k + level);
violators->push_back(k);
violators->push_back(i + level);
violators->push_back(i);
}
}
}
}
}
// Check that all pairwise differences are unique
bool CheckCostas(const std::vector<int64>& vars) {
std::vector<int> violators;
CheckConstraintViolators(vars, &violators);
return violators.empty();
}
// Cycles all possible permutations
class OrderedLNS: public BaseLNS {
public:
OrderedLNS(const std::vector<IntVar*>& vars, int free_elements) :
BaseLNS(vars.data(), vars.size()),
free_elements_(free_elements) {
index_ = 0;
// Start of with the first free_elements_ as a permutations, eg. 0,1,2,3,...
for (int i = 0; i < free_elements; ++i) {
index_ += i * pow(static_cast<double>(vars.size()),
free_elements - i - 1);
}
}
virtual bool NextFragment(std::vector<int>* const fragment) {
int dim = Size();
std::set<int> fragment_set;
do {
int work_index = index_;
fragment->clear();
fragment_set.clear();
for (int i = 0; i < free_elements_; ++i) {
int current_index = work_index % dim;
work_index = work_index / dim;
std::pair<std::set<int>::iterator, bool> ret =
fragment_set.insert(current_index);
// Check if element has been used before
if (ret.second) {
fragment->push_back(current_index);
} else {
break;
}
}
// Go to next possible permutation
++index_;
// Try again if a duplicate index is used
} while (fragment_set.size() < free_elements_);
return true;
}
private:
int index_;
const int free_elements_;
};
// RandomLNS is used for the local search and frees the
// number of elements specified in 'free_elements' randomly.
class RandomLNS: public BaseLNS {
public:
RandomLNS(const std::vector<IntVar*>& vars, int free_elements) :
BaseLNS(vars.data(), vars.size()),
free_elements_(free_elements),
rand_(ACMRandom::HostnamePidTimeSeed()) {
}
virtual bool NextFragment(std::vector<int>* const fragment) {
std::vector<int> weighted_elements;
std::vector<int64> values;
// Create weighted vector for randomizer. Add one of each possible elements
// to the weighted vector and then add more elements depending on the
// number of conflicts
for (int i = 0; i < Size(); ++i) {
values.push_back(Value(i));
// Insert each variable at least once.
weighted_elements.push_back(i);
}
CheckConstraintViolators(values, &weighted_elements);
int size = weighted_elements.size();
// Randomly insert elements in vector until no more options remain
while (fragment->size() < std::min(free_elements_, size)) {
const int index = weighted_elements[rand_.Next() % size];
fragment->push_back(index);
// Remove all elements with this index from weighted_elements
for (std::vector<int>::iterator pos = weighted_elements.begin();
pos != weighted_elements.end(); ) {
if (*pos == index) {
// Try to erase as many elements as possible at the same time
std::vector<int>::iterator end = pos;
while ((end + 1) != weighted_elements.end() && *(end + 1) == *pos) {
++end;
}
pos = weighted_elements.erase(pos, end);
} else {
++pos;
}
}
size = weighted_elements.size();
}
return true;
}
private:
const int free_elements_;
ACMRandom rand_;
};
class Evaluator {
public:
explicit Evaluator(const std::vector<IntVar*>& vars) : vars_(vars) {}
// Prefer the value with the smallest domain
int64 VarEvaluator(int64 index) {
return vars_[index]->Size();
}
// Penalize for each time the value appears in a different domain,
// as values have to be unique
int64 ValueEvaluator(int64 id, int64 value) {
int appearance = 0;
for (int i = 0; i < vars_.size(); ++i) {
if (i != id && vars_[i]->Contains(value)) {
++appearance;
}
}
return appearance;
}
private:
std::vector<IntVar*> vars_;
};
// Computes a Costas Array using soft constraints.
// Instead of enforcing that all distance vectors are unique, we
// minimize the number of duplicate distance vectors.
void CostasSoft(const int dim) {
Solver solver("Costas");
// For the matrix and for the count of occurrences of each possible distance
// for each stage
const int num_elements = dim + (2 * dim + 1) * (dim);
// create the variables
std::vector<IntVar*> vars;
solver.MakeIntVarArray(num_elements, -dim, dim, "var_", &vars);
// the costas matrix
std::vector<IntVar*> matrix(dim);
// number of occurrences per stage
std::vector<IntVar*> occurences;
// All possible values of the distance vector
// used to count the occurrence of all these values
std::vector<int64> possible_values(2 * dim + 1);
for (int64 i = -dim; i <= dim; ++i) {
possible_values[i + dim] = i;
}
int index = 0;
// Initialize the matrix that contains the coordinates of all '1' per row
for (; index < dim; ++index) {
matrix[index] = vars[index];
vars[index]->SetMin(1);
}
// First constraint for the elements in the Costas Matrix. We want
// them to be unique.
std::vector<IntVar*> matrix_count;
solver.MakeIntVarArray(2 * dim + 1, 0, dim, "matrix_count_", &matrix_count);
solver.AddConstraint(solver.MakeDistribute(matrix, possible_values,
matrix_count));
// Here we only consider the elements from 1 to dim.
for (int64 j = dim + 1; j <= 2 * dim; ++j) {
// Penalize if an element occurs more than once.
vars[index]
= solver.MakeSemiContinuousExpr(solver.MakeSum(matrix_count[j], -1),
0, 1)->Var();
occurences.push_back(vars[index++]);
}
// Count the number of duplicates for each stage
for (int i = 1; i < dim; ++i) {
std::vector<IntVar*> subset(dim - i);
// Initialize each stage
for (int j = 0; j < dim - i; ++j) {
IntVar* const diff = solver.MakeDifference(vars[j + i], vars[j])->Var();
subset[j] = diff;
}
// Count the number of occurrences for all possible values
std::vector<IntVar*> domain_count;
solver.MakeIntVarArray(2 * dim + 1, 0, dim, "domain_count_", &domain_count);
solver.AddConstraint(solver.MakeDistribute(subset,
possible_values,
domain_count));
// Penalize occurrences of more than one
for (int64 j = 0; j <= 2 * dim; ++j) {
vars[index] =
solver.MakeSemiContinuousExpr(solver.MakeSum(domain_count[j], -1),
0, dim - i)->Var();
occurences.push_back(vars[index++]);
}
}
// We would like to minimize the penalties that we just computed
IntVar* const objective_var = solver.MakeSum(occurences)->Var();
OptimizeVar* const total_duplicates = solver.MakeMinimize(objective_var, 1);
SearchMonitor* const log = solver.MakeSearchLog(1000, objective_var);
// Out of all solutions, we just want to store the last one.
SolutionCollector* const collector = solver.MakeLastSolutionCollector();
collector->Add(vars);
// The first solution that the local optimization is based on
Evaluator evaluator(matrix);
DecisionBuilder* const first_solution =
solver.MakePhase(matrix,
NewPermanentCallback(&evaluator,
&Evaluator::VarEvaluator),
NewPermanentCallback(&evaluator,
&Evaluator::ValueEvaluator));
SearchLimit* const search_time_limit =
solver.MakeLimit(FLAGS_timelimit, kint64max, kint64max, kint64max);
// Locally optimize solutions for LNS
SearchLimit* const fail_limit = solver.MakeLimit(kint64max, kint64max,
FLAGS_sublimit, kint64max);
DecisionBuilder* const subdecision_builder =
solver.MakeSolveOnce(first_solution, fail_limit);
std::vector<LocalSearchOperator*> localSearchOperators;
// Apply RandomLNS to free FLAGS_freevar variables at each stage
localSearchOperators.push_back(
solver.RevAlloc(new OrderedLNS(matrix, FLAGS_freevar)));
// Go through all possible permutations one by one
localSearchOperators.push_back(
solver.RevAlloc(new OrderedLNS(matrix, FLAGS_freeorderedvar)));
LocalSearchPhaseParameters* const ls_params =
solver.MakeLocalSearchPhaseParameters(
solver.ConcatenateOperators(localSearchOperators),
subdecision_builder);
DecisionBuilder* const second_phase =
solver.MakeLocalSearchPhase(matrix.data(),
matrix.size(),
first_solution,
ls_params);
// Try to find a solution
solver.Solve(second_phase,
collector,
log,
total_duplicates,
search_time_limit);
if (collector->solution_count() > 0) {
std::vector<int64> costas_matrix;
string output;
for (int n = 0; n < dim; ++n) {
const int64 v = collector->Value(0, vars[n]);
costas_matrix.push_back(v);
StringAppendF(&output, "%3lld", v);
}
if (!CheckCostas(costas_matrix)) {
LOG(INFO) << "No Costas Matrix found, closest solution displayed.";
}
LOG(INFO) << output;
} else {
LOG(INFO) << "No solution found";
}
}
// Computes a Costas Array.
void CostasHard(const int dim) {
SolverParameters parameters;
if (!FLAGS_export_profile.empty()) {
parameters.profile_level = SolverParameters::NORMAL_PROFILING;
}
Solver solver("costas", parameters);
// we need space for the matrix and for each pair
const int num_elements = dim + dim * (dim - 1) / 2;
// create the variables
std::vector<IntVar*> vars;
solver.MakeIntVarArray(num_elements, -dim, dim, "var", &vars);
std::vector<IntVar*> matrix(dim);
// Initialize the matrix that contains the coordinates of all '1' per row
for (int m = 0; m < dim; ++m) {
matrix[m] = vars[m];
vars[m]->SetMin(1);
}
solver.AddConstraint(solver.MakeAllDifferent(matrix));
int index = dim;
// Check that the pairwise difference is unique
for (int i = 1; i < dim; ++i) {
std::vector<IntVar*> subset(dim - i);
for (int j = 0; j < dim - i; ++j) {
IntVar* const diff = solver.MakeDifference(vars[j + i], vars[j])->Var();
vars[index++] = diff;
subset[j] = diff;
}
solver.AddConstraint(solver.MakeAllDifferent(subset));
}
DecisionBuilder* const db = solver.MakePhase(vars,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
solver.NewSearch(db);
if (solver.NextSolution()) {
std::vector<int64> costas_matrix;
string output;
for (int n = 0; n < dim; ++n) {
const int64 v = vars[n]->Value();
costas_matrix.push_back(v);
StringAppendF(&output, "%3lld", v);
}
LOG(INFO) << output << " (" << solver.wall_time() << "ms)";
CHECK(CheckCostas(costas_matrix)) <<
": Solution is not a valid Costas Matrix.";
} else {
LOG(INFO) << "No solution found.";
}
if (!FLAGS_export_profile.empty()) {
solver.ExportProfilingOverview(FLAGS_export_profile);
}
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
int min = 1;
int max = 10;
if (FLAGS_minsize != 0) {
min = FLAGS_minsize;
if (FLAGS_maxsize != 0) {
max = FLAGS_maxsize;
} else {
max = min;
}
}
for (int i = min; i <= max; ++i) {
LOG(INFO) << "Computing Costas Array for dim = " << i;
if (FLAGS_soft_constraints) {
operations_research::CostasSoft(i);
} else {
operations_research::CostasHard(i);
}
}
return 0;
}

111
examples/cpp/cryptarithm.cc Normal file
View File

@@ -0,0 +1,111 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Cryptoarithmetics problem
//
// Solves equation SEND + MORE = MONEY among numbers
// where each digit is represented by a letter.
//
// Solution:
// S=9; M=1; O=0; E=5; N=6; D=7; R=8; Y=2.
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "constraint_solver/constraint_solver.h"
namespace operations_research {
void Cryptoarithmetics() {
Solver solver("cryptarithm");
// model with carry
IntVar* s = solver.MakeIntVar(1, 9, "s");
IntVar* m = solver.MakeIntVar(1, 9, "m");
IntVar* o = solver.MakeIntVar(0, 9, "o");
IntVar* e = solver.MakeIntVar(0, 9, "e");
IntVar* n = solver.MakeIntVar(0, 9, "n");
IntVar* d = solver.MakeIntVar(0, 9, "d");
IntVar* r = solver.MakeIntVar(0, 9, "r");
IntVar* y = solver.MakeIntVar(0, 9, "y");
std::vector<IntVar*> letters;
letters.push_back(s);
letters.push_back(m);
letters.push_back(o);
letters.push_back(e);
letters.push_back(n);
letters.push_back(d);
letters.push_back(r);
letters.push_back(y);
solver.AddConstraint(solver.MakeAllDifferent(letters));
// carry variables
IntVar* c1 = solver.MakeIntVar(0, 1, "c1");
IntVar* c2 = solver.MakeIntVar(0, 1, "c2");
IntVar* c3 = solver.MakeIntVar(0, 1, "c3");
// initial constraint is separated into four small constraints
IntVar* v1 = solver.MakeSum(e, d)->Var();
IntVar* v2 = solver.MakeSum(y, solver.MakeProd(c1, 10))->Var();
solver.AddConstraint(solver.MakeEquality(v1, v2));
v1 = solver.MakeSum(solver.MakeSum(c1, n), r)->Var();
v2 = solver.MakeSum(e, solver.MakeProd(c2, 10))->Var();
solver.AddConstraint(solver.MakeEquality(v1, v2));
v1 = solver.MakeSum(solver.MakeSum(c2, e), o)->Var();
v2 = solver.MakeSum(n, solver.MakeProd(c3, 10))->Var();
solver.AddConstraint(solver.MakeEquality(v1, v2));
v1 = solver.MakeSum(solver.MakeSum(c3, s), m)->Var();
v2 = solver.MakeSum(o, solver.MakeProd(m, 10))->Var();
solver.AddConstraint(solver.MakeEquality(v1, v2));
DecisionBuilder* const db = solver.MakePhase(letters,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
solver.NewSearch(db);
if (solver.NextSolution()) {
CHECK_EQ(s->Value(), 9);
CHECK_EQ(m->Value(), 1);
CHECK_EQ(o->Value(), 0);
CHECK_EQ(e->Value(), 5);
CHECK_EQ(n->Value(), 6);
CHECK_EQ(d->Value(), 7);
CHECK_EQ(r->Value(), 8);
CHECK_EQ(y->Value(), 2);
LG << "S=" << s->Value();
LG << "M=" << m->Value();
LG << "O=" << o->Value();
LG << "E=" << e->Value();
LG << "N=" << n->Value();
LG << "D=" << d->Value();
LG << "R=" << r->Value();
LG << "Y=" << y->Value();
} else {
LG << "Cannot solve problem: number of failures " << solver.failures();
}
solver.EndSearch();
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::Cryptoarithmetics();
return 0;
}

282
examples/cpp/cvrptw.cc Normal file
View File

@@ -0,0 +1,282 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Capacitated Vehicle Routing Problem with Time Windows (and optional orders).
// A description of the problem can be found here:
// http://en.wikipedia.org/wiki/Vehicle_routing_problem.
// The variant which is tackled by this model includes a capacity dimension,
// time windows and optional orders, with a penalty cost if orders are not
// performed. For the sake of simplicty, orders are randomly located and
// distances are computed using the Manhattan distance. Distances are assumed
// to be in meters and times in seconds.
#include <vector>
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "base/stringprintf.h"
#include "constraint_solver/routing.h"
#include "base/random.h"
using operations_research::Assignment;
using operations_research::IntVar;
using operations_research::RoutingModel;
using operations_research::Solver;
using operations_research::ACMRandom;
using operations_research::StringAppendF;
using operations_research::StringPrintf;
using operations_research::scoped_array;
using operations_research::scoped_ptr;
DEFINE_int32(vrp_orders, 100, "Nodes in the problem.");
DEFINE_int32(vrp_vehicles, 20, "Size of Traveling Salesman Problem instance.");
DEFINE_bool(vrp_use_deterministic_random_seed, false,
"Use deterministic random seeds.");
const char* kTime = "Time";
const char* kCapacity = "Capacity";
// Random seed generator.
int32 GetSeed() {
if (FLAGS_vrp_use_deterministic_random_seed) {
return ACMRandom::DeterministicSeed();
} else {
return ACMRandom::HostnamePidTimeSeed();
}
}
// Location container, contains positions of orders and can be used to obtain
// Manhattan distances/times between locations.
class LocationContainer {
public:
explicit LocationContainer(int64 speed)
: randomizer_(GetSeed()), speed_(speed) {
CHECK_LT(0, speed_);
}
void AddLocation(int64 x, int64 y) {
locations_.push_back(Location(x, y));
}
void AddRandomLocation(int64 x_max, int64 y_max) {
AddLocation(randomizer_.Uniform(x_max + 1), randomizer_.Uniform(y_max + 1));
}
int64 ManhattanDistance(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return locations_[from].DistanceTo(locations_[to]);
}
int64 ManhattanTime(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return ManhattanDistance(from, to) / speed_;
}
private:
class Location {
public:
Location() : x_(0), y_(0) {}
Location(int64 x, int64 y) : x_(x), y_(y) {}
int64 DistanceTo(const Location& location) const {
return Abs(x_ - location.x_) + Abs(y_ - location.y_);
}
private:
static int64 Abs(int64 value) { return std::max(value, -value); }
int64 x_;
int64 y_;
};
ACMRandom randomizer_;
const int64 speed_;
ITIVector<RoutingModel::NodeIndex, Location> locations_;
};
// Random demand.
class RandomDemand {
public:
RandomDemand(int size, RoutingModel::NodeIndex depot)
: size_(size), depot_(depot) {
CHECK_LT(0, size_);
}
void Initialize() {
const int64 kDemandMax = 5;
const int64 kDemandMin = 1;
demand_.reset(new int64[size_]);
ACMRandom randomizer(GetSeed());
for (int order = 0; order < size_; ++order) {
if (order == depot_) {
demand_[order] = 0;
} else {
demand_[order] =
kDemandMin + randomizer.Uniform(kDemandMax - kDemandMin + 1);
}
}
}
int64 Demand(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return demand_[from.value()];
}
private:
scoped_array<int64> demand_;
const int size_;
const RoutingModel::NodeIndex depot_;
};
// Service time (proportional to demand) + transition time callback.
class ServiceTimePlusTransition {
public:
ServiceTimePlusTransition(int64 time_per_demand_unit,
RoutingModel::NodeEvaluator2* demand,
RoutingModel::NodeEvaluator2* transition_time)
: time_per_demand_unit_(time_per_demand_unit),
demand_(demand),
transition_time_(transition_time) {}
int64 Compute(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return time_per_demand_unit_ * demand_->Run(from, to)
+ transition_time_->Run(from, to);
}
private:
const int64 time_per_demand_unit_;
scoped_ptr<RoutingModel::NodeEvaluator2> demand_;
scoped_ptr<RoutingModel::NodeEvaluator2> transition_time_;
};
// Route plan displayer.
// TODO(user): Move the display code to the routing library.
void DisplayPlan(const RoutingModel& routing, const Assignment& plan) {
// Display plan cost.
string plan_output = StringPrintf("Cost %lld\n", plan.ObjectiveValue());
// Display dropped orders.
string dropped;
for (int order = 1; order < routing.nodes(); ++order) {
if (plan.Value(routing.NextVar(order)) == order) {
if (dropped.empty()) {
StringAppendF(&dropped, " %d", order);
} else {
StringAppendF(&dropped, ", %d", order);
}
}
}
if (!dropped.empty()) {
plan_output += "Dropped orders:" + dropped + "\n";
}
// Display actual output for each vehicle.
for (int route_number = 0;
route_number < routing.vehicles();
++route_number) {
int64 order = routing.Start(route_number);
StringAppendF(&plan_output, "Route %d: ", route_number);
if (routing.IsEnd(plan.Value(routing.NextVar(order)))) {
plan_output += "Empty\n";
} else {
while (!routing.IsEnd(order)) {
IntVar* load_var = routing.CumulVar(order, kCapacity);
IntVar* time_var = routing.CumulVar(order, kTime);
StringAppendF(&plan_output, "%lld Load(%lld) Time(%lld, %lld) -> ",
order,
plan.Value(load_var),
plan.Min(time_var),
plan.Max(time_var));
order = plan.Value(routing.NextVar(order));
}
IntVar* load_var = routing.CumulVar(order, kCapacity);
IntVar* time_var = routing.CumulVar(order, kTime);
StringAppendF(&plan_output, "%lld Load(%lld) Time(%lld, %lld)\n",
order,
plan.Value(load_var),
plan.Min(time_var),
plan.Max(time_var));
}
}
LG << plan_output;
}
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
CHECK_LT(0, FLAGS_vrp_orders) << "Specify an instance size greater than 0.";
CHECK_LT(0, FLAGS_vrp_vehicles) << "Specify a non-null vehicle fleet size.";
// VRP of size FLAGS_vrp_size.
// Nodes are indexed from 0 to FLAGS_vrp_orders, the starts and ends of
// the routes are at node 0.
const RoutingModel::NodeIndex kDepot(0);
RoutingModel routing(FLAGS_vrp_orders + 1, FLAGS_vrp_vehicles);
routing.SetDepot(kDepot);
// Setting first solution heuristic (cheapest addition).
routing.SetCommandLineOption("routing_first_solution", "PathCheapestArc");
// Disabling Large Neighborhood Search, comment out to activate it.
routing.SetCommandLineOption("routing_no_lns", "true");
// Setting up locations.
const int64 kXMax = 100000;
const int64 kYMax = 100000;
const int64 kSpeed = 10;
LocationContainer locations(kSpeed);
for (int location = 0; location <= FLAGS_vrp_orders; ++location) {
locations.AddRandomLocation(kXMax, kYMax);
}
// Setting the cost function.
routing.SetCost(NewPermanentCallback(&locations,
&LocationContainer::ManhattanDistance));
// Adding capacity dimension constraints.
const int64 kVehicleCapacity = 40;
const int64 kNullCapacitySlack = 0;
RandomDemand demand(routing.nodes(), kDepot);
demand.Initialize();
routing.AddDimension(NewPermanentCallback(&demand, &RandomDemand::Demand),
kNullCapacitySlack, kVehicleCapacity, kCapacity);
// Adding time dimension constraints.
const int64 kTimePerDemandUnit = 300;
const int64 kHorizon = 24 * 3600;
ServiceTimePlusTransition time(
kTimePerDemandUnit,
NewPermanentCallback(&demand, &RandomDemand::Demand),
NewPermanentCallback(&locations, &LocationContainer::ManhattanTime));
routing.AddDimension(
NewPermanentCallback(&time,
&ServiceTimePlusTransition::Compute),
kHorizon, kHorizon, kTime);
// Adding time windows.
ACMRandom randomizer(GetSeed());
const int64 kTWDuration = 5 * 3600;
for (int order = 1; order < routing.nodes(); ++order) {
const int64 start = randomizer.Uniform(kHorizon - kTWDuration);
routing.CumulVar(order, kTime)->SetRange(start, start + kTWDuration);
}
// Adding penalty costs to allow skipping orders.
const int64 kPenalty = 100000;
const RoutingModel::NodeIndex kFirstNodeAfterDepot(1);
for (RoutingModel::NodeIndex order = kFirstNodeAfterDepot;
order < routing.nodes(); ++order) {
std::vector<RoutingModel::NodeIndex> orders(1, order);
routing.AddDisjunction(orders, kPenalty);
}
// Solve, returns a solution if any (owned by RoutingModel).
const Assignment* solution = routing.Solve();
if (solution != NULL) {
DisplayPlan(routing, *solution);
} else {
LG << "No solution found.";
}
return 0;
}

184
examples/cpp/dimacs_assignment.cc Normal file
View File

@@ -0,0 +1,184 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include <algorithm>
#include <cstdlib>
#include "base/hash.h"
#include <string>
#include <vector>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/timer.h"
#include "algorithms/hungarian.h"
#include "graph/ebert_graph.h"
#include "graph/linear_assignment.h"
#include "cpp/parse_dimacs_assignment.h"
#include "cpp/print_dimacs_assignment.h"
DEFINE_bool(assignment_compare_hungarian, false,
"Compare result and speed against Hungarian method.");
DEFINE_string(assignment_problem_output_file, "",
"Print the problem to this file in DIMACS format (after layout "
"is optimized, if applicable).");
namespace operations_research {
CostValue BuildAndSolveHungarianInstance(
const LinearSumAssignment<ForwardStarGraph>& assignment) {
const ForwardStarGraph& graph = assignment.Graph();
typedef std::vector<double> HungarianRow;
typedef std::vector<HungarianRow> HungarianProblem;
HungarianProblem hungarian_cost;
hungarian_cost.resize(assignment.NumLeftNodes());
// First we have to find the biggest cost magnitude so we can
// initialize the arc costs that aren't really there.
CostValue largest_cost_magnitude = 0;
for (ForwardStarGraph::ArcIterator arc_it(graph);
arc_it.Ok();
arc_it.Next()) {
ArcIndex arc = arc_it.Index();
CostValue cost_magnitude = ::std::abs(assignment.ArcCost(arc));
largest_cost_magnitude = ::std::max(largest_cost_magnitude, cost_magnitude);
}
double missing_arc_cost = static_cast<double>((assignment.NumLeftNodes() *
largest_cost_magnitude) +
1);
for (HungarianProblem::iterator row = hungarian_cost.begin();
row != hungarian_cost.end();
++row) {
row->resize(assignment.NumNodes() - assignment.NumLeftNodes(),
missing_arc_cost);
}
// We're using a graph representation without forward arcs, so in
// order to use the generic ForwardStarGraph::ArcIterator we would
// need to increase our memory footprint by building the array of
// arc tails (since we need tails to build the input to the
// hungarian algorithm). We opt for the alternative of iterating
// over hte arcs via adjacency lists, which gives us the arc tails
// implicitly.
for (ForwardStarGraph::NodeIterator node_it(graph);
node_it.Ok();
node_it.Next()) {
NodeIndex node = node_it.Index();
NodeIndex tail = (node - ForwardStarGraph::kFirstNode);
for (ForwardStarGraph::OutgoingArcIterator arc_it(graph, node);
arc_it.Ok();
arc_it.Next()) {
ArcIndex arc = arc_it.Index();
NodeIndex head = (graph.Head(arc) - assignment.NumLeftNodes() -
ForwardStarGraph::kFirstNode);
double cost = static_cast<double>(assignment.ArcCost(arc));
hungarian_cost[tail][head] = cost;
}
}
hash_map<int, int> result;
hash_map<int, int> wish_this_could_be_null;
WallTimer timer;
VLOG(1) << "Beginning Hungarian method.";
timer.Start();
MinimizeLinearAssignment(hungarian_cost, &result, &wish_this_could_be_null);
double elapsed = timer.GetInMs() / 1000.0;
LOG(INFO) << "Hungarian result computed in " << elapsed << " seconds.";
double result_cost = 0.0;
for (int i = 0; i < assignment.NumLeftNodes(); ++i) {
int mate = result[i];
result_cost += hungarian_cost[i][mate];
}
return static_cast<CostValue>(result_cost);
}
void DisplayAssignment(
const LinearSumAssignment<ForwardStarGraph>& assignment) {
for (LinearSumAssignment<ForwardStarGraph>::BipartiteLeftNodeIterator
node_it(assignment);
node_it.Ok();
node_it.Next()) {
const NodeIndex left_node = node_it.Index();
const ArcIndex matching_arc = assignment.GetAssignmentArc(left_node);
const NodeIndex right_node = assignment.Head(matching_arc);
VLOG(5) << "assigned (" << left_node << ", " << right_node << "): "
<< assignment.ArcCost(matching_arc);
}
}
static const char* const kUsageTemplate = "usage: %s <filename>";
int solve_dimacs_assignment(int argc, char* argv[]) {
string usage;
if (argc < 1) {
usage = StringPrintf(kUsageTemplate, "solve_dimacs_assignment");
} else {
usage = StringPrintf(kUsageTemplate, argv[0]);
}
google::SetUsageMessage(usage);
google::ParseCommandLineFlags(&argc, &argv, true);
if (argc < 2) {
LOG(FATAL) << usage;
}
string error_message;
// Handle on the graph we will need to delete because the
// LinearSumAssignment object does not take ownership of it.
ForwardStarGraph* graph = NULL;
LinearSumAssignment<ForwardStarGraph>* assignment =
ParseDimacsAssignment(argv[1], &error_message, &graph);
if (assignment == NULL) {
LOG(FATAL) << error_message;
}
if (!FLAGS_assignment_problem_output_file.empty()) {
// The following tail array management stuff is done in a generic
// way so we can plug in different types of graphs for which the
// TailArrayManager template can be instantiated, even though we
// know the type of the graph explicitly. In this way, the type of
// the graph can be switched just by changing the graph type in
// this file and making no other changes to the code.
TailArrayManager<ForwardStarGraph> tail_array_manager(graph);
PrintDimacsAssignmentProblem(*assignment, tail_array_manager,
FLAGS_assignment_problem_output_file);
tail_array_manager.ReleaseTailArrayIfForwardGraph();
}
CostValue hungarian_cost = 0.0;
bool hungarian_solved = false;
if (FLAGS_assignment_compare_hungarian) {
hungarian_cost = BuildAndSolveHungarianInstance(*assignment);
hungarian_solved = true;
}
WallTimer timer;
timer.Start();
bool success = assignment->ComputeAssignment();
double elapsed = timer.GetInMs() / 1000.0;
if (success) {
CostValue cost = assignment->GetCost();
DisplayAssignment(*assignment);
LOG(INFO) << "Cost of optimum assignment: " << cost;
LOG(INFO) << "Computed in " << elapsed << " seconds.";
LOG(INFO) << assignment->StatsString();
if (hungarian_solved && (cost != hungarian_cost)) {
LOG(ERROR) << "Optimum cost mismatch: " << cost << " vs. "
<< hungarian_cost << ".";
}
} else {
LOG(WARNING) << "Given problem is infeasible.";
}
delete assignment;
delete graph;
return 0;
}
} // namespace operations_research
int main(int argc, char* argv[]) {
return ::operations_research::solve_dimacs_assignment(argc, argv);
}

814
examples/cpp/dobble_ls.cc Normal file
View File

@@ -0,0 +1,814 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// This problem is inspired by the Dobble game (aka Spot-It in the
// USA). In this game, we have 57 cards, 57 symbols, and 8 symbols
// per card. We want to assign symbols per card such that any two
// cards have exactly one symbol in common. These numbers can be
// generalized: we have N cards, each with K different symbols, and
// there are N different symbols overall.
//
// This is a feasability problem. We transform that into an
// optimization problem where we penalize cards whose intersection is
// of cardinality different from 1. A feasible solution of the
// original problem is a solution with a zero cost.
//
// Furthermore, we solve this problem using local search, and with a
// dedicated constraint.
//
// The purpose of the example is to demonstrates how to write local
// search operators and local search filters.
#include <algorithm>
#include <vector>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/concise_iterator.h"
#include "base/map-util.h"
#include "constraint_solver/constraint_solveri.h"
#include "util/bitset.h"
#include "base/random.h"
DEFINE_int32(symbols_per_card, 8, "Number of symbols per card.");
DEFINE_int32(ls_seed, 1, "Seed for the random number generator (used by "
"the Local Neighborhood Search).");
DEFINE_bool(use_filter, true, "Use filter in the local search to prune moves.");
DEFINE_int32(num_swaps, 4, "If num_swap > 0, the search for an optimal "
"solution will be allowed to use an operator that swaps the "
"symbols of up to num_swap pairs ((card1, symbol on card1), "
"(card2, symbol on card2)).");
DEFINE_int32(time_limit_in_ms,
60000,
"Time limit for the global search in ms.");
namespace operations_research {
// ----- Dedicated constraint to count the symbols shared by two cards -----
// This constraint maintains:
// sum_i(card1_symbol_vars[i]*card2_symbol_vars[i]) == count_var.
// with all card_symbol_vars[i] being boolean variables.
class SymbolsSharedByTwoCardsConstraint : public Constraint {
public:
// This constructor does not take any ownership on its arguments.
SymbolsSharedByTwoCardsConstraint(Solver* const solver,
const std::vector<IntVar*>& card1_symbol_vars,
const std::vector<IntVar*>& card2_symbol_vars,
IntVar* const num_symbols_in_common_var)
: Constraint(solver),
card1_symbol_vars_(card1_symbol_vars),
card2_symbol_vars_(card2_symbol_vars),
num_symbols_(card1_symbol_vars.size()),
num_symbols_in_common_var_(num_symbols_in_common_var) {
// Checks that cards have the same size.
CHECK_EQ(card1_symbol_vars.size(), card2_symbol_vars.size());
// Verify that we are really dealing with boolean variables.
for (int i = 0; i < num_symbols_; ++i) {
CHECK_GE(card1_symbol_vars_[i]->Min(), 0);
CHECK_GE(card2_symbol_vars_[i]->Min(), 0);
CHECK_LE(card1_symbol_vars_[i]->Max(), 1);
CHECK_LE(card2_symbol_vars_[i]->Max(), 1);
}
}
virtual ~SymbolsSharedByTwoCardsConstraint() {}
// Adds observers (named Demon) to variable events. These demons are
// responsible for implementing the propagation algorithm of the
// constraint.
virtual void Post() {
// Create a demon 'global_demon' that will bind events on
// variables to the calling of the 'InitialPropagate()' method. As
// this method is expensive, 'global_demon' has a low priority. As
// such, InitialPropagate will be called after all normal demons
// and constraints have reached a fixed point. Note
// that ownership of the 'global_demon' belongs to the solver.
Demon* const global_demon =
solver()->MakeDelayedConstraintInitialPropagateCallback(this);
// Attach to all variables.
for (int i = 0; i < num_symbols_; ++i) {
card1_symbol_vars_[i]->WhenBound(global_demon);
card2_symbol_vars_[i]->WhenBound(global_demon);
}
// Attach to cardinality variable.
num_symbols_in_common_var_->WhenBound(global_demon);
}
// This is the main propagation method.
//
// It scans all card1_symbol_vars * card2_symbol_vars and increments 3
// counters:
// - min_symbols_in_common if both booleans variables are bound to true.
// - max_symbols_in_common if both booleans are not bound to false.
//
// Then we know that num_symbols_in_common_var is in the range
// [min_symbols_in_common .. max_symbols_in_common].
//
// Now, if num_symbols_in_common_var->Max() ==
// min_symbols_in_common, then all products that contribute to
// max_symbols_in_common but not to min_symbols_in_common should be
// set to 0.
//
// Conversely, if num_symbols_in_common_var->Min() ==
// max_symbols_in_common, then all products that contribute to
// max_symbols_in_common should be set to 1.
virtual void InitialPropagate() {
int max_symbols_in_common = 0;
int min_symbols_in_common = 0;
for (int i = 0; i < num_symbols_; ++i) {
if (card1_symbol_vars_[i]->Min() == 1 &&
card2_symbol_vars_[i]->Min() == 1) {
min_symbols_in_common++;
}
if (card1_symbol_vars_[i]->Max() == 1 &&
card2_symbol_vars_[i]->Max() == 1) {
max_symbols_in_common++;
}
}
num_symbols_in_common_var_->SetRange(min_symbols_in_common,
max_symbols_in_common);
// If min_symbols_in_common == max_symbols_in_common, it means
// that num_symbols_in_common_var_ is already fully determined: we
// have nothing to do.
if (min_symbols_in_common == max_symbols_in_common) {
DCHECK_EQ(min_symbols_in_common, num_symbols_in_common_var_->Max());
DCHECK_EQ(min_symbols_in_common, num_symbols_in_common_var_->Min());
return;
}
if (num_symbols_in_common_var_->Max() == min_symbols_in_common) {
// All undecided product terms should be forced to 0.
for (int i = 0; i < num_symbols_; ++i) {
// If both Min() are 0, then we can't force either variable to
// be zero (even if we know that their product is zero),
// because either variable could be 1 as long as the other is
// 0.
if (card1_symbol_vars_[i]->Min() == 1 &&
card2_symbol_vars_[i]->Min() == 0) {
card2_symbol_vars_[i]->SetValue(0);
} else if (card2_symbol_vars_[i]->Min() == 1 &&
card1_symbol_vars_[i]->Min() == 0) {
card1_symbol_vars_[i]->SetValue(0);
}
}
} else if (num_symbols_in_common_var_->Min() == max_symbols_in_common) {
// All undecided product terms should be forced to 1.
for (int i = 0; i < num_symbols_; ++i) {
if (card1_symbol_vars_[i]->Max() == 1 &&
card2_symbol_vars_[i]->Max() == 1) {
// Note that we also force already-decided product terms,
// but this doesn't change anything.
card1_symbol_vars_[i]->SetValue(1);
card2_symbol_vars_[i]->SetValue(1);
}
}
}
}
private:
std::vector<IntVar*> card1_symbol_vars_;
std::vector<IntVar*> card2_symbol_vars_;
const int num_symbols_;
IntVar* const num_symbols_in_common_var_;
};
// Creates two integer variables: one that counts the number of
// symbols common to two cards, and one that counts the absolute
// difference between the first var and 1 (i.e. the violation of the
// objective). Returns the latter (both vars are owned by the Solver
// anyway).
IntVar* CreateViolationVar(Solver* const solver,
const std::vector<IntVar*>& card1_symbol_vars,
const std::vector<IntVar*>& card2_symbol_vars,
int num_symbols_per_card) {
IntVar* const num_symbols_in_common_var =
solver->MakeIntVar(0, num_symbols_per_card);
// RevAlloc transfers the ownership of the constraint to the solver.
solver->AddConstraint(
solver->RevAlloc(
new SymbolsSharedByTwoCardsConstraint(solver,
card1_symbol_vars,
card2_symbol_vars,
num_symbols_in_common_var)));
return solver->MakeAbs(solver->MakeSum(num_symbols_in_common_var, -1))->Var();
}
// ---------- Local Search ----------
// The "local search", or "local neighborhood search", works like
// this: starting from a given solution to the problem, other
// solutions in its neighborhood are generated from it, some of them
// might be selected (because they're better, for example) to become a
// starting point for other neighborhood searches, and so on.. The
// detailed search algorithm can vary and depends on the problem to
// solve.
//
// The fundamental building block for the local search is the "search
// operator", which has three fundamental methods in its API:
//
// - Its constructor, which keeps (mutable) references to the
// solver's internal variables (here, the card symbol variables).
//
// - OnStart(), which is called at the start of a local search, and
// after each solution (i.e. when the local search starts again from
// that new solution). The solver variables are supposed to represent
// a valid solution at this point. This method is used by the search
// operator to initialize its state and be ready to start the
// exploration of the neighborhood of the given solution.
//
// - MakeOneNeighbor(), which picks a neighbor of the initial
// solution, and changes the solver's internal variables accordingly
// to represent that new state.
//
// All local search operators on this problem will derive from the
// parent class below, which contains some shared code to store a
// compact representation of which symbols appeal on each cards.
class DobbleOperator : public IntVarLocalSearchOperator {
public:
DobbleOperator(const IntVar* const* card_symbol_vars,
int num_vars,
int num_cards,
int num_symbols,
int num_symbols_per_card)
: IntVarLocalSearchOperator(card_symbol_vars, num_vars),
num_cards_(num_cards),
num_symbols_(num_symbols),
num_symbols_per_card_(num_symbols_per_card),
symbols_per_card_(num_cards) {
CHECK_GT(num_cards, 0);
CHECK_GT(num_vars, 0);
CHECK_GT(num_symbols, 0);
CHECK_GT(num_symbols_per_card, 0);
for (int card = 0; card < num_cards; ++card) {
symbols_per_card_[card].assign(num_symbols_per_card, -1);
}
}
virtual ~DobbleOperator() {}
protected:
// OnStart() simply stores the current symbols per card in
// symbols_per_card_, and defers further initialization to the
// subclass's InitNeighborhoodSearch() method.
virtual void OnStart() {
for (int card = 0; card < num_cards_; ++card) {
int found = 0;
for (int symbol = 0; symbol < num_symbols_; ++symbol) {
if (Value(VarIndex(card, symbol)) == 1) {
symbols_per_card_[card][found++] = symbol;
}
}
DCHECK_EQ(num_symbols_per_card_, found);
}
InitNeighborhoodSearch();
}
virtual void InitNeighborhoodSearch() = 0;
// Find the index of the variable corresponding to the given symbol
// on the given card.
int VarIndex(int card, int symbol) {
return card * num_symbols_ + symbol;
}
// Move symbol1 from card1 to card2, and symbol2 from card2 to card1.
void SwapTwoSymbolsOnCards(int card1, int symbol1, int card2, int symbol2) {
SetValue(VarIndex(card1, symbol1), 0);
SetValue(VarIndex(card2, symbol2), 0);
SetValue(VarIndex(card1, symbol2), 1);
SetValue(VarIndex(card2, symbol1), 1);
}
const int num_cards_;
const int num_symbols_;
const int num_symbols_per_card_;
std::vector<std::vector<int> > symbols_per_card_;
};
// ----- Swap 2 symbols -----
// This operator explores *all* pairs (card1, some symbol on card1),
// (card2, some symbol on card2) and swaps the symbols between the two
// cards.
//
// Note that this could create invalid moves (for example, by adding a
// symbol to a card that already had it); see the DobbleFilter class
// below to see how we filter those out.
class SwapSymbols : public DobbleOperator {
public:
SwapSymbols(const IntVar* const* card_symbol_vars,
int num_vars,
int num_cards,
int num_symbols,
int num_symbols_per_card)
: DobbleOperator(card_symbol_vars,
num_vars,
num_cards,
num_symbols,
num_symbols_per_card),
current_card1_(-1),
current_card2_(-1),
current_symbol1_(-1),
current_symbol2_(-1) {
}
virtual ~SwapSymbols() {}
// Finds the next swap, returns false when it has finished.
virtual bool MakeOneNeighbor() {
if (!PickNextSwap()) {
VLOG(1) << "finished neighborhood";
return false;
}
const int symbol1 = symbols_per_card_[current_card1_][current_symbol1_];
const int symbol2 = symbols_per_card_[current_card2_][current_symbol2_];
SwapTwoSymbolsOnCards(current_card1_, symbol1, current_card2_, symbol2);
return true;
}
private:
// Reinit the exploration loop.
virtual void InitNeighborhoodSearch() {
current_card1_ = 0;
current_card2_ = 1;
current_symbol1_ = 0;
current_symbol2_ = -1;
}
// Compute the next move. It returns false when there are none.
bool PickNextSwap() {
current_symbol2_++;
if (current_symbol2_ == num_symbols_per_card_) {
current_symbol2_ = 0;
current_symbol1_++;
if (current_symbol1_ == num_symbols_per_card_) {
current_symbol1_ = 0;
current_card2_++;
if (current_card2_ == num_cards_) {
current_card1_++;
current_card2_ = current_card1_ + 1;
}
}
}
return current_card1_ < num_cards_ - 1;
}
int current_card1_;
int current_card2_;
int current_symbol1_;
int current_symbol2_;
};
// Multiple swaps of two symbols. This operator is an expanded version
// of the previous operator.
//
// At each step, it will pick a number num_swaps at random in
// [2 .. max_num_swaps], and then pick num_swaps random pairs (card1,
// some symbol on card1), (card2, some symbol on card2), and swap the
// symbols of each pair.
//
// As the search space (the "neighborhood") is huge, we use a
// randomized "infinite" version instead of an iterative, exhaustive
// one.
class SwapSymbolsOnCardPairs : public DobbleOperator {
public:
SwapSymbolsOnCardPairs(const IntVar* const* card_symbol_vars,
int num_vars,
int num_cards,
int num_symbols,
int num_symbols_per_card,
int max_num_swaps)
: DobbleOperator(card_symbol_vars,
num_vars,
num_cards,
num_symbols,
num_symbols_per_card),
rand_(FLAGS_ls_seed),
max_num_swaps_(max_num_swaps) {
CHECK_GE(max_num_swaps, 2);
}
virtual ~SwapSymbolsOnCardPairs() {}
protected:
virtual bool MakeOneNeighbor() {
const int num_swaps = rand_.Uniform(max_num_swaps_ - 1) + 2;
for (int i = 0; i < num_swaps; ++i) {
const int card_1 = rand_.Uniform(num_cards_);
const int symbol_index_1 = rand_.Uniform(num_symbols_per_card_);
const int symbol_1 = symbols_per_card_[card_1][symbol_index_1];
const int card_2 = rand_.Uniform(num_cards_);
const int symbol_index_2 = rand_.Uniform(num_symbols_per_card_);
const int symbol_2 = symbols_per_card_[card_2][symbol_index_2];
SwapTwoSymbolsOnCards(card_1, symbol_1, card_2, symbol_2);
}
return true;
}
virtual void InitNeighborhoodSearch() {}
private:
ACMRandom rand_;
const int max_num_swaps_;
};
// ----- Local Search Filter -----
// A filter is responsible for rejecting a local search move faster
// than what the propagation of the constraint solver would do.
// Its API consists in:
// - The constructor, which takes as input a reference to all the
// variables relevant to the filter.
// - OnSynchronize(), called at the beginning of the search and
// after each move to a new solution (when the local search
// restarts from it).
// - Accept(), which takes as input an attempted move (in the form
// of a Delta to tentatively apply to the variables), and returns
// true iff this move is found valid.
//
// To decide if a move is valid, first this DobbleFilter builds a
// bitmap of symbols per card. Then for each move, it updates the
// bitmap according to the move and checks the following constraints:
// - First, each card still has num_symbols_per_card symbols. - The
// cost of the assignment described by the move is better than the
// current one.
// After the check is done, the original bitmap is restored if the
// move was rejected, so as to be ready for the next evaluation.
//
// Please note that this filter uses a fixed size bitset and
// effectively limits the number of cards to 63, and thus the number
// of symbols per card to 8.
class DobbleFilter : public IntVarLocalSearchFilter {
public:
DobbleFilter(const IntVar* const* card_symbol_vars,
int num_vars,
int num_cards,
int num_symbols,
int num_symbols_per_card)
: IntVarLocalSearchFilter(card_symbol_vars, num_vars),
num_cards_(num_cards),
num_symbols_(num_symbols),
num_symbols_per_card_(num_symbols_per_card),
temporary_bitset_(0),
symbol_bitmask_per_card_(num_cards, 0),
violation_costs_(num_cards_, std::vector<int>(num_cards_, 0)) {}
// We build the current bitmap and the matrix of violation cost
// between any two cards.
virtual void OnSynchronize() {
symbol_bitmask_per_card_.assign(num_cards_, 0);
for (int card = 0; card < num_cards_; ++card) {
for (int symbol = 0; symbol < num_symbols_; ++symbol) {
if (Value(VarIndex(card, symbol))) {
SetBit64(&symbol_bitmask_per_card_[card], symbol);
}
}
}
for (int card1 = 0; card1 < num_cards_; ++card1) {
for (int card2 = 0; card2 < num_cards_; ++card2) {
violation_costs_[card1][card2] =
ViolationCost(BitCount64(symbol_bitmask_per_card_[card1] &
symbol_bitmask_per_card_[card2]));
}
}
DCHECK(CheckCards());
}
// The LocalSearchFilter::Accept() API also takes a deltadelta,
// which is the difference between the current delta and the last
// delta that was given to Accept() -- but we don't use it here.
virtual bool Accept(const Assignment* delta,
const Assignment* unused_deltadelta) {
const Assignment::IntContainer& solution_delta = delta->IntVarContainer();
const int solution_delta_size = solution_delta.Size();
// The input const Assignment* delta given to Accept() may
// actually contain "Deactivated" elements, which represent
// variables that have been freed -- they are not bound to a
// single value anymore. This happens with LNS-type (Large
// Neighborhood Search) LocalSearchOperator, which are not used in
// this example as of 2012-01; and we refer the reader to
// ./routing.cc for an example of such LNS-type operators.
//
// For didactical purposes, we will assume for a moment that a
// LNS-type operator might be applied. The Filter will still be
// called, but our DobbleFilter here won't be able to work, since
// it needs every variable to be bound (i.e. have a fixed value),
// in the assignment that it considers. Therefore, we include here
// a snippet of code that will detect if the input assignment is
// not fully bound. For further details, read ./routing.cc -- but
// we strongly advise the reader to first try and understand all
// of this file.
for (int i = 0; i < solution_delta_size; ++i) {
if (!solution_delta.Element(i).Activated()) {
VLOG(1)
<< "Element #" << i << " of the delta assignment given to"
<< " DobbleFilter::Accept() is not activated (i.e. its variable"
<< " is not bound to a single value anymore). This means that"
<< " we are in a LNS phase, and the DobbleFilter won't be able"
<< " to filter anything. Returning true.";
return true;
}
}
VLOG(1) << "No LNS, size = " << solution_delta_size;
// Collect the set of cards that have been modified by this move.
std::vector<int> touched_cards;
ComputeTouchedCards(solution_delta, &touched_cards);
// Check basic metrics to fail fast.
if (!CheckCards()) {
RestoreBitsetPerCard();
DCHECK(CheckCards());
VLOG(1) << "reject by size";
return false;
}
// Compute new cost.
const int cost_delta = ComputeNewCost(touched_cards);
// Reset the data structure to the state before the evaluation.
RestoreBitsetPerCard();
// And exit (this is only valid for a greedy descent and would
// reject valid moves in tabu search for instance).
if (cost_delta >= 0) {
VLOG(1) << "reject";
}
return cost_delta < 0;
}
private:
// Undo information after an evaluation.
struct UndoChange {
UndoChange(int c, uint64 b) : card(c), bitset(b) {}
int card;
uint64 bitset;
};
int VarIndex(int card, int symbol) {
return card * num_symbols_ + symbol;
}
void ClearBitset() {
temporary_bitset_ = 0;
}
// For each touched card, compare against all others to compute the
// delta in term of cost. We use an bitset to avoid counting twice
// between two cards appearing in the local search move.
int ComputeNewCost(const std::vector<int>& touched_cards) {
ClearBitset();
int cost_delta = 0;
for (int i = 0; i < touched_cards.size(); ++i) {
const int touched = touched_cards[i];
SetBit64(&temporary_bitset_, touched);
const uint64 card_bitset = symbol_bitmask_per_card_[touched];
const std::vector<int>& row_cost = violation_costs_[touched];
for (int other_card = 0; other_card < num_cards_; ++other_card) {
if (!IsBitSet64(&temporary_bitset_, other_card)) {
cost_delta +=
ViolationCost(BitCount64(card_bitset &
symbol_bitmask_per_card_[other_card]));
cost_delta -= row_cost[other_card];
}
}
}
return cost_delta;
}
// Collects all card indices appearing in the local search move.
void ComputeTouchedCards(const Assignment::IntContainer& solution_delta,
std::vector<int>* const touched_cards) {
ClearBitset();
const int solution_delta_size = solution_delta.Size();
const int kUnassigned = -1;
for (int index = 0; index < solution_delta_size; ++index) {
int64 touched_var = kUnassigned;
FindIndex(solution_delta.Element(index).Var(), &touched_var);
CHECK_NE(touched_var, kUnassigned);
const int card = touched_var / num_symbols_;
const int symbol = touched_var % num_symbols_;
const int new_value = solution_delta.Element(index).Value();
if (!IsBitSet64(&temporary_bitset_, card)) {
SaveRestoreInformation(card);
touched_cards->push_back(card);
SetBit64(&temporary_bitset_, card);
}
if (new_value) {
SetBit64(&symbol_bitmask_per_card_[card], symbol);
} else {
ClearBit64(&symbol_bitmask_per_card_[card], symbol);
}
}
}
// Undo all modifications done when evaluating a move.
void RestoreBitsetPerCard() {
for (int i = 0; i < restore_information_.size(); ++i) {
symbol_bitmask_per_card_[restore_information_[i].card] =
restore_information_[i].bitset;
}
restore_information_.clear();
}
// Stores undo information for a given card.
void SaveRestoreInformation(int card) {
restore_information_.push_back(UndoChange(card,
symbol_bitmask_per_card_[card]));
}
// Checks that after the local search move, each card would still have
// num_symbols_per_card symbols on it.
bool CheckCards() {
for (int i = 0; i < num_cards_; ++i) {
if (num_symbols_per_card_ != BitCount64(symbol_bitmask_per_card_[i])) {
VLOG(1) << "card " << i << " has bitset of size "
<< BitCount64(symbol_bitmask_per_card_[i]);
return false;
}
}
return true;
}
int ViolationCost(uint64 cardinality) const {
return (cardinality > 0 ? cardinality - 1 : 1);
}
const int num_cards_;
const int num_symbols_;
const int num_symbols_per_card_;
uint64 temporary_bitset_;
std::vector<uint64> symbol_bitmask_per_card_;
std::vector<std::vector<int> > violation_costs_;
std::vector<UndoChange> restore_information_;
};
// ----- Main Method -----
void SolveDobble(int num_cards, int num_symbols, int num_symbols_per_card) {
LOG(INFO) << "Solving dobble assignment problem:";
LOG(INFO) << " - " << num_cards << " cards";
LOG(INFO) << " - " << num_symbols << " symbols";
LOG(INFO) << " - " << num_symbols_per_card << " symbols per card";
// Creates the solver.
Solver solver("dobble");
// Creates the matrix of boolean variables (cards * symbols).
std::vector<std::vector<IntVar*> > card_symbol_vars(num_cards);
std::vector<IntVar*> all_card_symbol_vars;
for (int card_index = 0; card_index < num_cards; ++card_index) {
solver.MakeBoolVarArray(num_symbols,
StringPrintf("card_%i_", card_index),
&card_symbol_vars[card_index]);
for (int symbol_index = 0;
symbol_index < num_symbols;
++symbol_index) {
all_card_symbol_vars.push_back(
card_symbol_vars[card_index][symbol_index]);
}
}
// Creates cardinality intersection variables and remember the
// violation variables.
std::vector<IntVar*> violation_vars;
for (int card1 = 0; card1 < num_cards; ++card1) {
for (int card2 = 0; card2 < num_cards; ++card2) {
if (card1 != card2) {
violation_vars.push_back(CreateViolationVar(&solver,
card_symbol_vars[card1],
card_symbol_vars[card2],
num_symbols_per_card));
}
}
}
// Create the objective variable.
IntVar* const objective_var = solver.MakeSum(violation_vars)->Var();
// Add constraint: there must be exactly num_symbols_per_card
// symbols per card.
for (int card = 0; card < num_cards; ++card) {
solver.AddConstraint(solver.MakeSumEquality(card_symbol_vars[card],
num_symbols_per_card));
}
// IMPORTANT OPTIMIZATION:
// Add constraint: each symbol appears on exactly
// num_symbols_per_card cards (i.e. symbols are evenly
// distributed). This constraint is actually redundant, because it
// is a (non-trivial) consequence of the other constraints and of
// the model. But adding it makes the search go faster.
for (int symbol_index = 0; symbol_index < num_symbols; ++symbol_index) {
std::vector<IntVar*> tmp;
for (int card_index = 0; card_index < num_cards; ++card_index) {
tmp.push_back(card_symbol_vars[card_index][symbol_index]);
}
solver.AddConstraint(solver.MakeSumEquality(tmp, num_symbols_per_card));
}
// Search.
LOG(INFO) << "Solving with Local Search";
LOG(INFO) << " - time limit = " << FLAGS_time_limit_in_ms << " ms";
// Start a DecisionBuilder phase to find a first solution, using the
// strategy "Pick some random, yet unassigned card symbol variable
// and set its value to 1".
DecisionBuilder* const build_db = solver.MakePhase(
all_card_symbol_vars,
Solver::CHOOSE_RANDOM, // Solver::IntVarStrategy
Solver::ASSIGN_MAX_VALUE); // Solver::IntValueStrategy
// Creates local search operators.
std::vector<LocalSearchOperator*> operators;
LocalSearchOperator* const switch_operator =
solver.RevAlloc(new SwapSymbols(all_card_symbol_vars.data(),
all_card_symbol_vars.size(),
num_cards,
num_symbols,
num_symbols_per_card));
operators.push_back(switch_operator);
LOG(INFO) << " - add switch operator";
if (FLAGS_num_swaps > 0) {
LocalSearchOperator* const swaps_operator =
solver.RevAlloc(new SwapSymbolsOnCardPairs(all_card_symbol_vars.data(),
all_card_symbol_vars.size(),
num_cards,
num_symbols,
num_symbols_per_card,
FLAGS_num_swaps));
operators.push_back(swaps_operator);
LOG(INFO) << " - add swaps operator with at most "
<< FLAGS_num_swaps << " swaps";
}
// Creates filter.
std::vector<LocalSearchFilter*> filters;
if (FLAGS_use_filter) {
filters.push_back(solver.RevAlloc(
new DobbleFilter(all_card_symbol_vars.data(),
all_card_symbol_vars.size(),
num_cards,
num_symbols,
num_symbols_per_card)));
}
// Main decision builder that regroups the first solution decision
// builder and the combination of local search operators and
// filters.
DecisionBuilder* const final_db =
solver.MakeLocalSearchPhase(
all_card_symbol_vars,
build_db,
solver.MakeLocalSearchPhaseParameters(
solver.ConcatenateOperators(operators, true),
NULL, // Sub decision builder, not needed here.
NULL, // Limit the search for improving move, we will stop
// the exploration of the local search at the first
// improving solution (first accept).
filters));
std::vector<SearchMonitor*> monitors;
// Optimize var search monitor.
OptimizeVar* const optimize = solver.MakeMinimize(objective_var, 1);
monitors.push_back(optimize);
// Search log.
SearchMonitor* const log = solver.MakeSearchLog(100000, optimize);
monitors.push_back(log);
// Search limit.
SearchLimit* const time_limit =
solver.MakeLimit(FLAGS_time_limit_in_ms, kint64max, kint64max, kint64max);
monitors.push_back(time_limit);
// And solve!
solver.Solve(final_db, monitors);
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
// These constants comes directly from the dobble game.
// There are actually 55 cards, but we can create up to 57 cards.
const int kSymbolsPerCard = FLAGS_symbols_per_card;
const int kCards = kSymbolsPerCard * (kSymbolsPerCard - 1) + 1;
const int kSymbols = kCards;
operations_research::SolveDobble(kCards,
kSymbols,
kSymbolsPerCard);
return 0;
}

91
examples/cpp/flow_api.cc Normal file
View File

@@ -0,0 +1,91 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "graph/ebert_graph.h"
#include "graph/max_flow.h"
#include "graph/min_cost_flow.h"
namespace operations_research {
// ----- Min Cost Flow -----
// Test on a 4x4 matrix. Example taken from
// http://www.ee.oulu.fi/~mpa/matreng/eem1_2-1.htm
void MinCostFlowOn4x4Matrix() {
LOG(INFO) << "Min Cost Flow on 4x4 Matrix";
const int kNumSources = 4;
const int kNumTargets = 4;
const CostValue kCost[kNumSources][kNumTargets] = {
{ 90, 75, 75, 80 },
{ 35, 85, 55, 65 },
{ 125, 95, 90, 105 },
{ 45, 110, 95, 115 }
};
const CostValue kExpectedCost = 275;
StarGraph graph(kNumSources + kNumTargets, kNumSources * kNumTargets);
MinCostFlow min_cost_flow(&graph);
for (NodeIndex source = 0; source < kNumSources; ++source) {
for (NodeIndex target = 0; target < kNumTargets; ++target) {
ArcIndex arc = graph.AddArc(source, kNumSources + target);
min_cost_flow.SetArcUnitCost(arc, kCost[source][target]);
min_cost_flow.SetArcCapacity(arc, 1);
}
}
for (NodeIndex source = 0; source < kNumSources; ++source) {
min_cost_flow.SetNodeSupply(source, 1);
}
for (NodeIndex target = 0; target < kNumTargets; ++target) {
min_cost_flow.SetNodeSupply(kNumSources + target, -1);
}
CHECK(min_cost_flow.Solve());
CHECK_EQ(MinCostFlow::OPTIMAL, min_cost_flow.status());
CostValue total_flow_cost = min_cost_flow.GetOptimalCost();
CHECK_EQ(kExpectedCost, total_flow_cost);
}
// ----- Max Flow -----
void MaxFeasibleFlow() {
LOG(INFO) << "Max Feasible Flow";
const int kNumNodes = 6;
const int kNumArcs = 9;
const NodeIndex kTail[kNumArcs] = { 0, 0, 0, 0, 1, 2, 3, 3, 4 };
const NodeIndex kHead[kNumArcs] = { 1, 2, 3, 4, 3, 4, 4, 5, 5 };
const FlowQuantity kCapacity[kNumArcs] = { 5, 8, 5, 3, 4, 5, 6, 6, 4 };
const FlowQuantity kExpectedFlow[kNumArcs] = { 4, 4, 2, 0, 4, 4, 0, 6, 4 };
const FlowQuantity kExpectedTotalFlow = 10;
StarGraph graph(kNumNodes, kNumArcs);
MaxFlow max_flow(&graph, 0, kNumNodes - 1);
for (int i = 0; i < kNumArcs; ++i) {
ArcIndex arc = graph.AddArc(kTail[i], kHead[i]);
max_flow.SetArcCapacity(arc, kCapacity[i]);
}
CHECK(max_flow.Solve());
CHECK_EQ(MaxFlow::OPTIMAL, max_flow.status());
FlowQuantity total_flow = max_flow.GetOptimalFlow();
CHECK_EQ(total_flow, kExpectedTotalFlow);
for (int i = 0; i < kNumArcs; ++i) {
CHECK_EQ(kExpectedFlow[i], max_flow.Flow(i)) << " i = " << i;
}
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::MinCostFlowOn4x4Matrix();
operations_research::MaxFeasibleFlow();
return 0;
}

582
examples/cpp/global_arith.cc Normal file
View File

@@ -0,0 +1,582 @@
// Copyright 2010 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "base/stl_util.h"
#include "constraint_solver/constraint_solveri.h"
#include "examples/global_arith.h"
namespace operations_research {
class ArithmeticPropagator;
// ----- SubstitutionMap -----
class SubstitutionMap {
public:
struct Offset { // to_replace = var_index + offset
Offset() : var_index(-1), offset(0) {}
Offset(int v, int64 o) : var_index(v), offset(o) {}
int var_index;
int64 offset;
};
void AddSubstitution(int left_var, int right_var, int64 right_offset) {
// TODO(lperron) : Perform transitive closure.
substitutions_[left_var] = Offset(right_var, right_offset);
}
void ProcessAllSubstitutions(Callback3<int, int, int64>* const hook) {
for (hash_map<int, Offset>::const_iterator it = substitutions_.begin();
it != substitutions_.end();
++it) {
hook->Run(it->first, it->second.var_index, it->second.offset);
}
}
private:
hash_map<int, Offset> substitutions_;
};
// ----- Bounds -----
struct Bounds {
Bounds() : lb(kint64min), ub(kint64max) {}
Bounds(int64 l, int64 u) : lb(l), ub(u) {}
void Intersect(int64 new_lb, int64 new_ub) {
lb = std::max(lb, new_lb);
ub = std::min(ub, new_ub);
}
void Intersect(const Bounds& other) {
Intersect(other.lb, other.ub);
}
void Union(int64 new_lb, int64 new_ub) {
lb = std::min(lb, new_lb);
ub = std::max(ub, new_ub);
}
void Union(const Bounds& other) {
Union(other.lb, other.ub);
}
bool IsEqual(const Bounds& other) {
return (ub == other.ub && lb == other.lb);
}
bool IsIncluded(const Bounds& other) {
return (ub <= other.ub && lb >= other.lb);
}
int64 lb;
int64 ub;
};
// ----- BoundsStore -----
class BoundsStore {
public:
BoundsStore(vector<Bounds>* initial_bounds)
: initial_bounds_(initial_bounds) {}
void SetRange(int var_index, int64 lb, int64 ub) {
hash_map<int, Bounds>::iterator it = modified_bounds_.find(var_index);
if (it == modified_bounds_.end()) {
Bounds new_bounds(lb, ub);
const Bounds& initial = (*initial_bounds_)[var_index];
new_bounds.Intersect(initial);
if (!new_bounds.IsEqual(initial)) {
modified_bounds_.insert(make_pair(var_index, new_bounds));
}
} else {
it->second.Intersect(lb, ub);
}
}
void Clear() {
modified_bounds_.clear();
}
const hash_map<int, Bounds>& modified_bounds() const {
return modified_bounds_;
}
vector<Bounds>* initial_bounds() const { return initial_bounds_; }
void Apply() {
for (hash_map<int, Bounds>::const_iterator it = modified_bounds_.begin();
it != modified_bounds_.end();
++it) {
(*initial_bounds_)[it->first] = it->second;
}
}
private:
vector<Bounds>* initial_bounds_;
hash_map<int, Bounds> modified_bounds_;
};
// ----- ArithmeticConstraint -----
class ArithmeticConstraint {
public:
virtual ~ArithmeticConstraint() {}
const vector<int>& vars() const { return vars_; }
virtual bool Propagate(BoundsStore* const store) = 0;
virtual void Replace(int to_replace, int var, int64 offset) = 0;
virtual bool Deduce(ArithmeticPropagator* const propagator) const = 0;
virtual string DebugString() const = 0;
private:
const vector<int> vars_;
};
// ----- ArithmeticPropagator -----
class ArithmeticPropagator : PropagationBaseObject {
public:
ArithmeticPropagator(Solver* const solver, Demon* const demon)
: PropagationBaseObject(solver), demon_(demon) {}
void ReduceProblem() {
for (int constraint_index = 0;
constraint_index < constraints_.size();
++constraint_index) {
if (constraints_[constraint_index]->Deduce(this)) {
protected_constraints_.insert(constraint_index);
}
}
scoped_ptr<Callback3<int, int, int64> > hook(
NewPermanentCallback(this,
&ArithmeticPropagator::ProcessOneSubstitution));
substitution_map_.ProcessAllSubstitutions(hook.get());
}
void Post() {
for (int constraint_index = 0;
constraint_index < constraints_.size();
++constraint_index) {
const vector<int>& vars = constraints_[constraint_index]->vars();
for (int var_index = 0; var_index < vars.size(); ++var_index) {
dependencies_[vars[var_index]].push_back(constraint_index);
}
}
}
void InitialPropagate() {
}
void Update(int var_index) {
Enqueue(demon_);
}
void AddConstraint(ArithmeticConstraint* const ct) {
constraints_.push_back(ct);
}
void AddVariable(int64 lb, int64 ub) {
bounds_.push_back(Bounds(lb, ub));
}
const vector<IntVar*> vars() const { return vars_; }
int VarIndex(IntVar* const var) {
hash_map<IntVar*, int>::const_iterator it = var_map_.find(var);
if (it == var_map_.end()) {
const int index = var_map_.size();
var_map_[var] = index;
return index;
} else {
return it->second;
}
}
void AddSubstitution(int left_var, int right_var, int64 right_offset) {
substitution_map_.AddSubstitution(left_var, right_var, right_offset);
}
void AddNewBounds(int var_index, int64 lb, int64 ub) {
bounds_[var_index].Intersect(lb, ub);
}
void ProcessOneSubstitution(int left_var, int right_var, int64 right_offset) {
for (int constraint_index = 0;
constraint_index < constraints_.size();
++constraint_index) {
if (!ContainsKey(protected_constraints_, constraint_index)) {
ArithmeticConstraint* const constraint = constraints_[constraint_index];
constraint->Replace(left_var, right_var, right_offset);
}
}
}
void PrintModel() {
LOG(INFO) << "Vars:";
for (int i = 0; i < bounds_.size(); ++i) {
LOG(INFO) << " var<" << i << "> = [" << bounds_[i].lb
<< " .. " << bounds_[i].ub << "]";
}
LOG(INFO) << "Constraints";
for (int i = 0; i < constraints_.size(); ++i) {
LOG(INFO) << " " << constraints_[i]->DebugString();
}
}
private:
Demon* const demon_;
vector<IntVar*> vars_;
hash_map<IntVar*, int> var_map_;
vector<ArithmeticConstraint*> constraints_;
vector<Bounds> bounds_;
vector<vector<int> > dependencies_; // from var indices to constraints.
SubstitutionMap substitution_map_;
hash_set<int> protected_constraints_;
};
// ----- Custom Constraints -----
class RowConstraint : public ArithmeticConstraint {
public:
RowConstraint(int64 lb, int64 ub) : lb_(lb), ub_(ub) {}
virtual ~RowConstraint() {}
void AddTerm(int var_index, int64 coefficient) {
// TODO(lperron): Check not present.
coefficients_[var_index] = coefficient;
}
virtual bool Propagate(BoundsStore* const store) {
return true;
}
virtual void Replace(int to_replace, int var, int64 offset) {
hash_map<int, int64>::iterator find_other = coefficients_.find(to_replace);
if (find_other != coefficients_.end()) {
hash_map<int, int64>::iterator find_var = coefficients_.find(var);
const int64 other_coefficient = find_other->second;
if (lb_ != kint64min) {
lb_ += other_coefficient * offset;
}
if (ub_ != kint64max) {
ub_ += other_coefficient * offset;
}
coefficients_.erase(find_other);
if (find_var == coefficients_.end()) {
coefficients_[var] = other_coefficient;
} else {
find_var->second += other_coefficient;
if (find_var->second == 0) {
coefficients_.erase(find_var);
}
}
}
}
virtual bool Deduce(ArithmeticPropagator* const propagator) const {
// Deduce Simple translation from one var to another.
if (lb_ == ub_ && coefficients_.size() == 2) {
hash_map<int, int64>::const_iterator it = coefficients_.begin();
const int var1 = it->first;
const int64 coeff1 = it->second;
++it;
const int var2 = it->first;
const int64 coeff2 = it->second;
++it;
CHECK(it == coefficients_.end());
if (coeff1 == 1 && coeff2 == -1) {
propagator->AddSubstitution(var1, var2, lb_);
return true;
} else if (coeff1 == -1 && coeff2 && 1) {
propagator->AddSubstitution(var2, var1, lb_);
return true;
}
}
return false;
}
virtual string DebugString() const {
string output = "(";
bool first = true;
for (hash_map<int, int64>::const_iterator it = coefficients_.begin();
it != coefficients_.end();
++it) {
if (it->second != 0) {
if (first) {
first = false;
if (it->second == 1) {
output += StringPrintf("var<%d>", it->first);
} else if (it->second == -1) {
output += StringPrintf("-var<%d>", it->first);
} else {
output += StringPrintf("%lld*var<%d>", it->second, it->first);
}
} else if (it->second == 1) {
output += StringPrintf(" + var<%d>", it->first);
} else if (it->second == -1) {
output += StringPrintf(" - var<%d>", it->first);
} else if (it->second > 0) {
output += StringPrintf(" + %lld*var<%d>", it->second, it->first);
} else {
output += StringPrintf(" - %lld*var<%d>", -it->second, it->first);
}
}
}
if (lb_ == ub_) {
output += StringPrintf(" == %lld)", ub_);
} else if (lb_ == kint64min) {
output += StringPrintf(" <= %lld)", ub_);
} else if (ub_ == kint64max) {
output += StringPrintf(" >= %lld)", lb_);
} else {
output += StringPrintf(" in [%lld .. %lld])", lb_, ub_);
}
return output;
}
private:
hash_map<int, int64> coefficients_;
int64 lb_;
int64 ub_;
};
class OrConstraint : public ArithmeticConstraint {
public:
OrConstraint(ArithmeticConstraint* const left,
ArithmeticConstraint* const right)
: left_(left), right_(right) {}
virtual ~OrConstraint() {}
virtual bool Propagate(BoundsStore* const store) {
return true;
}
virtual void Replace(int to_replace, int var, int64 offset) {
left_->Replace(to_replace, var, offset);
right_->Replace(to_replace, var, offset);
}
virtual bool Deduce(ArithmeticPropagator* const propagator) const {
return false;
}
virtual string DebugString() const {
return StringPrintf("Or(%s, %s)",
left_->DebugString().c_str(),
right_->DebugString().c_str());
}
private:
ArithmeticConstraint* const left_;
ArithmeticConstraint* const right_;
};
// ----- GlobalArithmeticConstraint -----
GlobalArithmeticConstraint::GlobalArithmeticConstraint(Solver* const solver)
: Constraint(solver),
propagator_(NULL) {
propagator_.reset(new ArithmeticPropagator(
solver,
solver->MakeDelayedConstraintInitialPropagateCallback(this)));
}
GlobalArithmeticConstraint::~GlobalArithmeticConstraint() {
STLDeleteElements(&constraints_);
}
void GlobalArithmeticConstraint::Post() {
const vector<IntVar*>& vars = propagator_->vars();
for (int var_index = 0; var_index < vars.size(); ++var_index) {
Demon* const demon =
MakeConstraintDemon1(solver(),
this,
&GlobalArithmeticConstraint::Update,
"Update",
var_index);
vars[var_index]->WhenRange(demon);
}
LOG(INFO) << "----- Before reduction -----";
propagator_->PrintModel();
LOG(INFO) << "----- After reduction -----";
propagator_->ReduceProblem();
propagator_->PrintModel();
LOG(INFO) << "---------------------------";
propagator_->Post();
}
void GlobalArithmeticConstraint::InitialPropagate() {
propagator_->InitialPropagate();
}
void GlobalArithmeticConstraint::Update(int var_index) {
propagator_->Update(var_index);
}
ConstraintRef GlobalArithmeticConstraint::MakeScalProdGreaterOrEqualConstant(
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant) {
RowConstraint* const constraint = new RowConstraint(constant, kint64max);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), coefficients[index]);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeScalProdLessOrEqualConstant(
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant) {
RowConstraint* const constraint = new RowConstraint(kint64min, constant);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), coefficients[index]);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeScalProdEqualConstant(
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant) {
RowConstraint* const constraint = new RowConstraint(constant, constant);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), coefficients[index]);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeSumGreaterOrEqualConstant(
const vector<IntVar*> vars,
int64 constant) {
RowConstraint* const constraint = new RowConstraint(constant, kint64max);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), 1);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeSumLessOrEqualConstant(
const vector<IntVar*> vars, int64 constant) {
RowConstraint* const constraint = new RowConstraint(kint64min, constant);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), 1);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeSumEqualConstant(
const vector<IntVar*> vars, int64 constant) {
RowConstraint* const constraint = new RowConstraint(constant, constant);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), 1);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeRowConstraint(
int64 lb,
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 ub) {
RowConstraint* const constraint = new RowConstraint(lb, ub);
for (int index = 0; index < vars.size(); ++index) {
constraint->AddTerm(VarIndex(vars[index]), coefficients[index]);
}
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
int64 ub) {
RowConstraint* const constraint = new RowConstraint(lb, ub);
constraint->AddTerm(VarIndex(v1), coeff1);
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
int64 ub) {
RowConstraint* const constraint = new RowConstraint(lb, ub);
constraint->AddTerm(VarIndex(v1), coeff1);
constraint->AddTerm(VarIndex(v2), coeff2);
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
IntVar* const v3,
int64 coeff3,
int64 ub) {
RowConstraint* const constraint = new RowConstraint(lb, ub);
constraint->AddTerm(VarIndex(v1), coeff1);
constraint->AddTerm(VarIndex(v2), coeff2);
constraint->AddTerm(VarIndex(v3), coeff3);
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
IntVar* const v3,
int64 coeff3,
IntVar* const v4,
int64 coeff4,
int64 ub) {
RowConstraint* const constraint = new RowConstraint(lb, ub);
constraint->AddTerm(VarIndex(v1), coeff1);
constraint->AddTerm(VarIndex(v2), coeff2);
constraint->AddTerm(VarIndex(v3), coeff3);
constraint->AddTerm(VarIndex(v4), coeff4);
return Store(constraint);
}
ConstraintRef GlobalArithmeticConstraint::MakeOrConstraint(
ConstraintRef left_ref,
ConstraintRef right_ref) {
OrConstraint* const constraint =
new OrConstraint(constraints_[left_ref.index()],
constraints_[right_ref.index()]);
return Store(constraint);
}
void GlobalArithmeticConstraint::Add(ConstraintRef ref) {
propagator_->AddConstraint(constraints_[ref.index()]);
}
int GlobalArithmeticConstraint::VarIndex(IntVar* const var) {
hash_map<IntVar*, int>::const_iterator it = var_indices_.find(var);
if (it == var_indices_.end()) {
const int new_index = var_indices_.size();
var_indices_.insert(make_pair(var, new_index));
propagator_->AddVariable(var->Min(), var->Max());
return new_index;
} else {
return it->second;
}
}
ConstraintRef GlobalArithmeticConstraint::Store(
ArithmeticConstraint* const constraint) {
const int constraint_index = constraints_.size();
constraints_.push_back(constraint);
return ConstraintRef(constraint_index);
}
} // namespace operations_research

100
examples/cpp/global_arith.h Normal file
View File

@@ -0,0 +1,100 @@
// Copyright 2010 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef EXAMPLES_GLOBAL_ARITH_H_
#define EXAMPLES_GLOBAL_ARITH_H_
namespace operations_research {
class ArithmeticPropagator;
class ArithmeticConstraint;
class ConstraintRef {
public:
ConstraintRef(int index) : index_(index) {}
int index() const { return index_; }
private:
const int index_;
};
class GlobalArithmeticConstraint : public Constraint {
public:
GlobalArithmeticConstraint(Solver* const solver);
virtual ~GlobalArithmeticConstraint();
virtual void Post();
virtual void InitialPropagate();
void Update(int var_index);
ConstraintRef MakeScalProdGreaterOrEqualConstant(
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant);
ConstraintRef MakeScalProdLessOrEqualConstant(
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant);
ConstraintRef MakeScalProdEqualConstant(const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 constant);
ConstraintRef MakeSumGreaterOrEqualConstant(const vector<IntVar*> vars,
int64 constant);
ConstraintRef MakeSumLessOrEqualConstant(const vector<IntVar*> vars,
int64 constant);
ConstraintRef MakeSumEqualConstant(const vector<IntVar*> vars,
int64 constant);
ConstraintRef MakeRowConstraint(int64 lb,
const vector<IntVar*> vars,
const vector<int64> coefficients,
int64 ub);
ConstraintRef MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
int64 ub);
ConstraintRef MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
int64 ub);
ConstraintRef MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
IntVar* const v3,
int64 coeff3,
int64 ub);
ConstraintRef MakeRowConstraint(int64 lb,
IntVar* const v1,
int64 coeff1,
IntVar* const v2,
int64 coeff2,
IntVar* const v3,
int64 coeff3,
IntVar* const v4,
int64 coeff4,
int64 ub);
ConstraintRef MakeOrConstraint(ConstraintRef left_constraint_index,
ConstraintRef right_constraint_index);
void Add(ConstraintRef constraint_index);
private:
int VarIndex(IntVar* const var);
ConstraintRef Store(ArithmeticConstraint* const constraint);
scoped_ptr<ArithmeticPropagator> propagator_;
hash_map<IntVar*, int> var_indices_;
vector<ArithmeticConstraint*> constraints_;
};
} // namespace operations_research
#endif // EXAMPLES_GLOBAL_ARITH_H_

107
examples/cpp/golomb.cc Normal file
View File

@@ -0,0 +1,107 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Golomb ruler problem
//
// find minimal ruler so that the differences between ticks are unique.
//
// First solutions:
// 0, 1
// 0, 1, 3
// 0, 1, 4, 6
// 0, 1, 4, 9, 11
// 0, 1, 4, 10, 12, 17
// 0, 1, 4, 10, 18, 23, 25
#include <cstdio>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "base/stringprintf.h"
#include "constraint_solver/constraint_solver.h"
DEFINE_bool(print, false, "If true, print the minimal solution.");
DEFINE_int32(size, 0,
"Size of the problem. If equal to 0, will test several increasing sizes.");
static const int kBestSolutions[] = {
0, 1, 3, 6, 11, 17, 25, 34, 44, 55, 72, 85,
// just for the optimistics ones, the rest:
106, 127, 151, 177, 199, 216, 246
};
static const int kKnownSolutions = 19;
namespace operations_research {
void GolombRuler(int size) {
CHECK_GE(size, 1);
Solver s("golomb");
// model
std::vector<IntVar*> ticks(size);
ticks[0] = s.MakeIntConst(0); // X(0) = 0
const int64 max = 1 + size * size * size;
for (int i = 1; i < size; ++i) {
ticks[i] = s.MakeIntVar(1, max, StringPrintf("X%02d", i));
}
std::vector<IntVar*> diffs;
for (int i = 0; i < size; ++i) {
for (int j = i + 1; j < size; ++j) {
IntVar* const diff = s.MakeDifference(ticks[j], ticks[i])->Var();
diffs.push_back(diff);
diff->SetMin(1);
}
}
s.AddConstraint(s.MakeAllDifferent(diffs));
OptimizeVar* const length = s.MakeMinimize(ticks[size-1], 1);
SolutionCollector* const collector = s.MakeLastSolutionCollector();
collector->Add(ticks);
DecisionBuilder* const db = s.MakePhase(ticks,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
s.Solve(db, collector, length); // go!
CHECK_EQ(collector->solution_count(), 1);
const int64 result = collector->Value(0, ticks[size-1]);
const int num_failures = collector->failures(0);
printf("N = %d, optimal length = %d (fails:%d)\n",
size, static_cast<int>(result), num_failures);
if (size - 1 < kKnownSolutions) {
CHECK_EQ(result, kBestSolutions[size - 1]);
}
if (FLAGS_print) {
for (int i = 0; i < size; ++i) {
const int64 tick = collector->Value(0, ticks[i]);
printf("%d ", static_cast<int>(tick));
}
printf("\n");
}
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_size != 0) {
operations_research::GolombRuler(FLAGS_size);
} else {
for (int n = 1; n < 11; ++n) {
operations_research::GolombRuler(n);
}
}
return 0;
}

80
examples/cpp/integer_programming.cc Normal file
View File

@@ -0,0 +1,80 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Integer programming example that shows how to use the API.
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "linear_solver/linear_solver.h"
namespace operations_research {
void RunIntegerProgrammingExample(
MPSolver::OptimizationProblemType optimization_problem_type) {
MPSolver solver("IntegerProgrammingExample", optimization_problem_type);
const double infinity = solver.infinity();
// x1 and x2 are integer non-negative variables.
MPVariable* const x1 = solver.MakeIntVar(0.0, infinity, "x1");
MPVariable* const x2 = solver.MakeIntVar(0.0, infinity, "x2");
// Minimize x1 + 2 * x2.
solver.SetObjectiveCoefficient(x1, 1);
solver.SetObjectiveCoefficient(x2, 2);
// 2 * x2 + 3 * x1 >= 17.
MPConstraint* const c0 = solver.MakeRowConstraint(17, infinity);
c0->SetCoefficient(x1, 3);
c0->SetCoefficient(x2, 2);
const MPSolver::ResultStatus result_status = solver.Solve();
// Check that the problem has an optimal solution.
if (result_status != MPSolver::OPTIMAL) {
LOG(FATAL) << "The problem does not have an optimal solution!";
}
LOG(INFO) << "Problem solved in " << solver.wall_time() << " milliseconds";
// The objective value of the solution.
LOG(INFO) << "Optimal objective value = " << solver.objective_value();
// The value of each variable in the solution.
LOG(INFO) << "x1 = " << x1->solution_value();
LOG(INFO) << "x2 = " << x2->solution_value();
LOG(INFO) << "Advanced usage:";
LOG(INFO) << "Problem solved in " << solver.nodes()
<< " branch-and-bound nodes";
}
void RunAllExamples() {
#if defined(USE_GLPK)
LOG(INFO) << "---- Integer programming example with GLPK ----";
RunIntegerProgrammingExample(MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING);
#endif
#if defined(USE_CBC)
LOG(INFO) << "---- Integer programming example with CBC ----";
RunIntegerProgrammingExample(MPSolver::CBC_MIXED_INTEGER_PROGRAMMING);
#endif
#if defined(USE_SCIP)
LOG(INFO) << "---- Integer programming example with SCIP ----";
RunIntegerProgrammingExample(MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING);
#endif
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::RunAllExamples();
return 0;
}

189
examples/cpp/jobshop.cc Normal file
View File

@@ -0,0 +1,189 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// This model implements a simple jobshop problem.
//
// A jobshop is a standard scheduling problem where you must schedule a
// set of jobs on a set of machines. Each job is a sequence of tasks
// (a task can only start when the preceding task finished), each of
// which occupies a single specific machine during a specific
// duration. Therefore, a job is simply given by a sequence of pairs
// (machine id, duration).
// The objective is to minimize the 'makespan', which is the duration
// between the start of the first task (across all machines) and the
// completion of the last task (across all machines).
//
// This will be modelled by sets of intervals variables (see class
// IntervalVar in constraint_solver/constraint_solver.h), one per
// task, representing the [start_time, end_time] of the task. Tasks
// in the same job will be linked by precedence constraints. Tasks on
// the same machine will be covered by Sequence constraints.
//
// Search will then be applied on the sequence constraints.
#include <cstdio>
#include <cstdlib>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "constraint_solver/constraint_solver.h"
#include "cpp/jobshop.h"
DEFINE_string(
data_file,
"",
"Required: input file description the scheduling problem to solve, "
"in our jssp format:\n"
" - the first line is \"instance <instance name>\"\n"
" - the second line is \"<number of jobs> <number of machines>\"\n"
" - then one line per job, with a single space-separated "
"list of \"<machine index> <duration>\"\n"
"note: jobs with one task are not supported");
DEFINE_int32(time_limit_in_ms, 0, "Time limit in ms, 0 means no limit.");
namespace operations_research {
void Jobshop(const JobShopData& data) {
Solver solver("jobshop");
const int machine_count = data.machine_count();
const int job_count = data.job_count();
const int horizon = data.horizon();
// ----- Creates all Intervals and vars -----
// Stores all tasks attached interval variables per job.
std::vector<std::vector<IntervalVar*> > jobs_to_tasks(job_count);
// machines_to_tasks stores the same interval variables as above, but
// grouped my machines instead of grouped by jobs.
std::vector<std::vector<IntervalVar*> > machines_to_tasks(machine_count);
// Creates all individual interval variables.
for (int job_id = 0; job_id < job_count; ++job_id) {
const std::vector<JobShopData::Task>& tasks = data.TasksOfJob(job_id);
for (int task_index = 0; task_index < tasks.size(); ++task_index) {
const JobShopData::Task& task = tasks[task_index];
CHECK_EQ(job_id, task.job_id);
const string name = StringPrintf("J%dM%dI%dD%d",
task.job_id,
task.machine_id,
task_index,
task.duration);
IntervalVar* const one_task =
solver.MakeFixedDurationIntervalVar(0,
horizon,
task.duration,
false,
name);
jobs_to_tasks[task.job_id].push_back(one_task);
machines_to_tasks[task.machine_id].push_back(one_task);
}
}
// ----- Creates model -----
// Creates precedences inside jobs.
for (int job_id = 0; job_id < job_count; ++job_id) {
const int task_count = jobs_to_tasks[job_id].size();
for (int task_index = 0; task_index < task_count - 1; ++task_index) {
IntervalVar* const t1 = jobs_to_tasks[job_id][task_index];
IntervalVar* const t2 = jobs_to_tasks[job_id][task_index + 1];
Constraint* const prec =
solver.MakeIntervalVarRelation(t2, Solver::STARTS_AFTER_END, t1);
solver.AddConstraint(prec);
}
}
// Adds disjunctive constraints on unary resources.
for (int machine_id = 0; machine_id < machine_count; ++machine_id) {
solver.AddConstraint(
solver.MakeDisjunctiveConstraint(machines_to_tasks[machine_id]));
}
// Creates sequences variables on machines. A sequence variable is a
// dedicated variable whose job is to sequence interval variables.
std::vector<SequenceVar*> all_sequences;
for (int machine_id = 0; machine_id < machine_count; ++machine_id) {
const string name = StringPrintf("Machine_%d", machine_id);
SequenceVar* const sequence =
solver.MakeSequenceVar(machines_to_tasks[machine_id], name);
all_sequences.push_back(sequence);
}
// Creates array of end_times of jobs.
std::vector<IntVar*> all_ends;
for (int job_id = 0; job_id < job_count; ++job_id) {
const int task_count = jobs_to_tasks[job_id].size();
IntervalVar* const task = jobs_to_tasks[job_id][task_count - 1];
all_ends.push_back(task->EndExpr()->Var());
}
// Objective: minimize the makespan (maximum end times of all tasks)
// of the problem.
IntVar* const objective_var = solver.MakeMax(all_ends)->Var();
OptimizeVar* const objective_monitor = solver.MakeMinimize(objective_var, 1);
// ----- Search monitors and decision builder -----
// This decision builder will rank all tasks on all machines.
DecisionBuilder* const sequence_phase =
solver.MakePhase(all_sequences, Solver::SEQUENCE_DEFAULT);
// After the ranking of tasks, the schedule is still loose and any
// task can be postponed at will. But, because the problem is now a PERT
// (http://en.wikipedia.org/wiki/Program_Evaluation_and_Review_Technique),
// we can schedule each task at its earliest start time. This is
// conveniently done by fixing the objective variable to its
// minimum value.
DecisionBuilder* const obj_phase =
solver.MakePhase(objective_var,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
// The main decision builder (ranks all tasks, then fixes the
// objective_variable).
DecisionBuilder* const main_phase =
solver.Compose(sequence_phase, obj_phase);
// Search log.
const int kLogFrequency = 1000000;
SearchMonitor* const search_log =
solver.MakeSearchLog(kLogFrequency, objective_monitor);
SearchLimit* limit = NULL;
if (FLAGS_time_limit_in_ms > 0) {
limit = solver.MakeTimeLimit(FLAGS_time_limit_in_ms);
}
// Search.
solver.Solve(main_phase, search_log, objective_monitor, limit);
}
} // namespace operations_research
static const char kUsage[] =
"Usage: see flags.\nThis program runs a simple job shop optimization "
"output besides the debug LOGs of the solver.";
int main(int argc, char **argv) {
google::SetUsageMessage(kUsage);
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_data_file.empty()) {
LOG(FATAL) << "Please supply a data file with --data_file=";
}
operations_research::JobShopData data;
data.Load(FLAGS_data_file);
operations_research::Jobshop(data);
return 0;
}

219
examples/cpp/jobshop.h Normal file
View File

@@ -0,0 +1,219 @@
// Copyright 2011 Google Inc. All Rights Reserved.
//
// This model implements a simple jobshop problem.
//
// A jobshop is a standard scheduling problem where you must schedule a
// set of jobs on a set of machines. Each job is a sequence of tasks
// (a task can only start when the preceding task finished), each of
// which occupies a single specific machine during a specific
// duration. Therefore, a job is simply given by a sequence of pairs
// (machine id, duration).
// The objective is to minimize the 'makespan', which is the duration
// between the start of the first task (across all machines) and the
// completion of the last task (across all machines).
//
// This will be modelled by sets of intervals variables (see class
// IntervalVar in constraint_solver/constraint_solver.h), one per
// task, representing the [start_time, end_time] of the task. Tasks
// in the same job will be linked by precedence constraints. Tasks on
// the same machine will be covered by Sequence constraints.
#ifndef OR_TOOLS_EXAMPLES_JOBSHOP_H_
#define OR_TOOLS_EXAMPLES_JOBSHOP_H_
#include <cstdio>
#include <cstdlib>
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/strtoint.h"
#include "base/file.h"
#include "base/filelinereader.h"
#include "base/split.h"
namespace operations_research {
// ----- JobShopData -----
// A JobShopData parses data files and stores all data internally for
// easy retrieval.
class JobShopData {
public:
// A task is the basic block of a jobshop.
struct Task {
Task(int j, int m, int d) : job_id(j), machine_id(m), duration(d) {}
int job_id;
int machine_id;
int duration;
};
enum ProblemType {
UNDEFINED,
JSSP,
TAILLARD
};
enum TaillardState {
START,
JOBS_READ,
MACHINES_READ,
SEED_READ,
JOB_ID_READ,
JOB_LENGTH_READ,
JOB_READ
};
JobShopData()
: name_(""),
machine_count_(0),
job_count_(0),
horizon_(0),
current_job_index_(0),
problem_type_(UNDEFINED),
taillard_state_(START) {}
~JobShopData() {}
// Parses a file in jssp or taillard format and loads the model. See the flag
// --data_file for a description of the format. Note that the format
// is only partially checked: bad inputs might cause undefined
// behavior.
void Load(const string& filename) {
FileLineReader reader(filename.c_str());
reader.set_line_callback(NewPermanentCallback(
this,
&JobShopData::ProcessNewLine));
reader.Reload();
if (!reader.loaded_successfully()) {
LOG(ERROR) << "Could not open jobshop file";
}
}
// The number of machines in the jobshop.
int machine_count() const { return machine_count_; }
// The number of jobs in the jobshop.
int job_count() const { return job_count_; }
// The name of the jobshop instance.
const string& name() const { return name_; }
// The horizon of the workshop (the sum of all durations), which is
// a trivial upper bound of the optimal make_span.
int horizon() const { return horizon_; }
// Returns the tasks of a job, ordered by precedence.
const std::vector<Task>& TasksOfJob(int job_id) const {
return all_tasks_[job_id];
}
private:
void ProcessNewLine(char* const line) {
// TODO(user): more robust logic to support single-task jobs.
static const char kWordDelimiters[] = " ";
std::vector<string> words;
SplitStringUsing(line, kWordDelimiters, &words);
switch (problem_type_) {
case UNDEFINED: {
if (words.size() == 2 && words[0] == "instance") {
problem_type_ = JSSP;
LOG(INFO) << "Reading jssp instance " << words[1];
name_ = words[1];
} else if (words.size() == 1 && atoi32(words[0]) > 0) {
problem_type_ = TAILLARD;
taillard_state_ = JOBS_READ;
job_count_ = atoi32(words[0]);
CHECK_GT(job_count_, 0);
all_tasks_.resize(job_count_);
}
break;
}
case JSSP: {
if (words.size() == 2) {
job_count_ = atoi32(words[0]);
machine_count_ = atoi32(words[1]);
CHECK_GT(machine_count_, 0);
CHECK_GT(job_count_, 0);
LOG(INFO) << machine_count_ << " machines and "
<< job_count_ << " jobs";
all_tasks_.resize(job_count_);
}
if (words.size() > 2 && machine_count_ != 0) {
CHECK_EQ(words.size(), machine_count_ * 2);
for (int i = 0; i < machine_count_; ++i) {
const int machine_id = atoi32(words[2 * i]);
const int duration = atoi32(words[2 * i + 1]);
AddTask(current_job_index_, machine_id, duration);
}
current_job_index_++;
}
break;
}
case TAILLARD: {
switch (taillard_state_) {
case START: {
LOG(FATAL) << "Should not be here";
break;
}
case JOBS_READ: {
CHECK_EQ(1, words.size());
machine_count_ = atoi32(words[0]);
CHECK_GT(machine_count_, 0);
taillard_state_ = MACHINES_READ;
break;
}
case MACHINES_READ: {
CHECK_EQ(1, words.size());
const int seed = atoi32(words[0]);
LOG(INFO) << "Taillard instance with " << job_count_
<< " jobs, and " << machine_count_
<< " machines, generated with a seed of " << seed;
taillard_state_ = SEED_READ;
break;
}
case SEED_READ:
case JOB_READ: {
CHECK_EQ(1, words.size());
current_job_index_ = atoi32(words[0]);
taillard_state_ = JOB_ID_READ;
break;
}
case JOB_ID_READ: {
CHECK_EQ(1, words.size());
taillard_state_ = JOB_LENGTH_READ;
break;
}
case JOB_LENGTH_READ: {
CHECK_EQ(machine_count_, words.size());
for (int i = 0; i < machine_count_; ++i) {
const int duration = atoi32(words[i]);
AddTask(current_job_index_, i, duration);
}
taillard_state_ = JOB_READ;
break;
}
}
break;
}
}
}
void AddTask(int job_id, int machine_id, int duration) {
all_tasks_[job_id].push_back(Task(job_id, machine_id, duration));
horizon_ += duration;
}
string name_;
int machine_count_;
int job_count_;
int horizon_;
std::vector<std::vector<Task> > all_tasks_;
int current_job_index_;
ProblemType problem_type_;
TaillardState taillard_state_;
};
} // namespace operations_research
#endif // OR_TOOLS_EXAMPLES_JOBSHOP_H_

460
examples/cpp/jobshop_ls.cc Normal file
View File

@@ -0,0 +1,460 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// This model implements a simple jobshop problem.
//
// A jobshop is a standard scheduling problem where you must schedule a
// set of jobs on a set of machines. Each job is a sequence of tasks
// (a task can only start when the preceding task finished), each of
// which occupies a single specific machine during a specific
// duration. Therefore, a job is simply given by a sequence of pairs
// (machine id, duration).
// The objective is to minimize the 'makespan', which is the duration
// between the start of the first task (across all machines) and the
// completion of the last task (across all machines).
//
// This will be modelled by sets of intervals variables (see class
// IntervalVar in constraint_solver/constraint_solver.h), one per
// task, representing the [start_time, end_time] of the task. Tasks
// in the same job will be linked by precedence constraints. Tasks on
// the same machine will be covered by Sequence constraints.
//
// Search will be implemented as local search on the sequence variables.
#include <cstdio>
#include <cstdlib>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/bitmap.h"
#include "constraint_solver/constraint_solver.h"
#include "constraint_solver/constraint_solveri.h"
#include "cpp/jobshop.h"
DEFINE_string(
data_file,
"",
"Required: input file description the scheduling problem to solve, "
"in our jssp format:\n"
" - the first line is \"instance <instance name>\"\n"
" - the second line is \"<number of jobs> <number of machines>\"\n"
" - then one line per job, with a single space-separated "
"list of \"<machine index> <duration>\"\n"
"note: jobs with one task are not supported");
DEFINE_int32(time_limit_in_ms, 60000, "Time limit in ms, 0 means no limit.");
DEFINE_int32(shuffle_length, 4, "Length of sub-sequences to shuffle LS.");
DEFINE_int32(sub_sequence_length, 4,
"Length of sub-sequences to relax in LNS.");
DEFINE_int32(lns_seed, 1, "Seed of the LNS random search");
DEFINE_int32(lns_limit, 30,
"Limit the size of the search tree in a LNS fragment");
namespace operations_research {
// ----- Exchange 2 intervals on a sequence variable -----
class SwapIntervals : public SequenceVarLocalSearchOperator {
public:
SwapIntervals(const SequenceVar* const* vars, int size)
: SequenceVarLocalSearchOperator(vars, size),
current_var_(-1),
current_first_(-1),
current_second_(-1) {}
virtual ~SwapIntervals() {}
virtual bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta) {
CHECK_NOTNULL(delta);
while (true) {
RevertChanges(true);
if (!Increment()) {
VLOG(1) << "finished neighborhood";
return false;
}
std::vector<int> sequence = Sequence(current_var_);
const int tmp = sequence[current_first_];
sequence[current_first_] = sequence[current_second_];
sequence[current_second_] = tmp;
SetForwardSequence(current_var_, sequence);
if (ApplyChanges(delta, deltadelta)) {
VLOG(1) << "Delta = " << delta->DebugString();
return true;
}
}
return false;
}
protected:
virtual void OnStart() {
VLOG(1) << "start neighborhood";
current_var_ = 0;
current_first_ = 0;
current_second_ = 0;
}
private:
bool Increment() {
const SequenceVar* const var = Var(current_var_);
if (++current_second_ >= var->size()) {
if (++current_first_ >= var->size() - 1) {
current_var_++;
current_first_ = 0;
}
current_second_ = current_first_ + 1;
}
return current_var_ < Size();
}
int current_var_;
int current_first_;
int current_second_;
};
// ----- Shuffle a fixed-length sub-sequence on one sequence variable -----
class ShuffleIntervals : public SequenceVarLocalSearchOperator {
public:
ShuffleIntervals(const SequenceVar* const* vars, int size, int max_length)
: SequenceVarLocalSearchOperator(vars, size),
max_length_(max_length),
current_var_(-1),
current_first_(-1),
current_index_(-1),
current_length_(-1) {}
virtual ~ShuffleIntervals() {}
virtual bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta) {
CHECK_NOTNULL(delta);
while (true) {
RevertChanges(true);
if (!Increment()) {
VLOG(1) << "finished neighborhood";
return false;
}
std::vector<int> sequence = Sequence(current_var_);
std::vector<int> sequence_backup(current_length_);
for (int i = 0; i < current_length_; ++i) {
sequence_backup[i] = sequence[i + current_first_];
}
for (int i = 0; i < current_length_; ++i) {
sequence[i + current_first_] =
sequence_backup[current_permutation_[i]];
}
SetForwardSequence(current_var_, sequence);
if (ApplyChanges(delta, deltadelta)) {
VLOG(1) << "Delta = " << delta->DebugString();
return true;
}
}
return false;
}
protected:
virtual void OnStart() {
VLOG(1) << "start neighborhood";
current_var_ = 0;
current_first_ = 0;
current_index_ = -1;
current_length_ = std::min(Var(current_var_)->size(), max_length_);
current_permutation_.resize(current_length_);
for (int i = 0; i < current_length_; ++i) {
current_permutation_[i] = i;
}
}
private:
bool Increment() {
if (!std::next_permutation(current_permutation_.begin(),
current_permutation_.end())) {
if (++current_first_ >= Var(current_var_)->size() - current_length_) {
if (++current_var_ >= Size()) {
return false;
}
current_first_ = 0;
current_length_ = std::min(Var(current_var_)->size(), max_length_);
current_permutation_.resize(current_length_);
}
current_index_ = 0;
}
return true;
}
const int max_length_;
int current_var_;
int current_first_;
int current_index_;
int current_length_;
std::vector<int> current_permutation_;
};
// ----- LNS Operator -----
class SequenceLns : public SequenceVarLocalSearchOperator {
public:
SequenceLns(const SequenceVar* const* vars,
int size,
int seed,
int max_length)
: SequenceVarLocalSearchOperator(vars, size),
random_(seed),
max_length_(max_length) {}
virtual ~SequenceLns() {}
virtual bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta) {
CHECK_NOTNULL(delta);
while (true) {
RevertChanges(true);
if (random_.Uniform(2) == 0) {
FreeTimeWindow();
} else {
FreeTwoResources();
}
if (ApplyChanges(delta, deltadelta)) {
VLOG(1) << "Delta = " << delta->DebugString();
return true;
}
}
return false;
}
private:
void FreeTimeWindow() {
for (int i = 0; i < Size(); ++i) {
std::vector<int> sequence = Sequence(i);
const int current_length =
std::min(static_cast<int>(sequence.size()), max_length_);
const int start_position =
random_.Uniform(sequence.size() - current_length);
std::vector<int> forward;
for (int j = 0; j < start_position; ++j) {
forward.push_back(sequence[j]);
}
std::vector<int> backward;
for (int j = sequence.size() - 1;
j >= start_position + current_length;
--j) {
backward.push_back(sequence[j]);
}
SetForwardSequence(i, forward);
SetBackwardSequence(i, backward);
}
}
void FreeTwoResources() {
std::vector<int> free_sequence;
SetForwardSequence(random_.Uniform(Size()), free_sequence);
SetForwardSequence(random_.Uniform(Size()), free_sequence);
}
ACMRandom random_;
const int max_length_;
};
// ----- Model and Solve -----
void JobshopLs(const JobShopData& data) {
Solver solver("jobshop");
const int machine_count = data.machine_count();
const int job_count = data.job_count();
const int horizon = data.horizon();
// ----- Creates all Intervals and vars -----
// Stores all tasks attached interval variables per job.
std::vector<std::vector<IntervalVar*> > jobs_to_tasks(job_count);
// machines_to_tasks stores the same interval variables as above, but
// grouped my machines instead of grouped by jobs.
std::vector<std::vector<IntervalVar*> > machines_to_tasks(machine_count);
// Creates all individual interval variables.
for (int job_id = 0; job_id < job_count; ++job_id) {
const std::vector<JobShopData::Task>& tasks = data.TasksOfJob(job_id);
for (int task_index = 0; task_index < tasks.size(); ++task_index) {
const JobShopData::Task& task = tasks[task_index];
CHECK_EQ(job_id, task.job_id);
const string name = StringPrintf("J%dM%dI%dD%d",
task.job_id,
task.machine_id,
task_index,
task.duration);
IntervalVar* const one_task =
solver.MakeFixedDurationIntervalVar(0,
horizon,
task.duration,
false,
name);
jobs_to_tasks[task.job_id].push_back(one_task);
machines_to_tasks[task.machine_id].push_back(one_task);
}
}
// ----- Creates model -----
// Creates precedences inside jobs.
for (int job_id = 0; job_id < job_count; ++job_id) {
const int task_count = jobs_to_tasks[job_id].size();
for (int task_index = 0; task_index < task_count - 1; ++task_index) {
IntervalVar* const t1 = jobs_to_tasks[job_id][task_index];
IntervalVar* const t2 = jobs_to_tasks[job_id][task_index + 1];
Constraint* const prec =
solver.MakeIntervalVarRelation(t2, Solver::STARTS_AFTER_END, t1);
solver.AddConstraint(prec);
}
}
// Adds disjunctive constraints on unary resources.
for (int machine_id = 0; machine_id < machine_count; ++machine_id) {
solver.AddConstraint(
solver.MakeDisjunctiveConstraint(machines_to_tasks[machine_id]));
}
// Creates sequences variables on machines. A sequence variable is a
// dedicated variable whose job is to sequence interval variables.
std::vector<SequenceVar*> all_sequences;
for (int machine_id = 0; machine_id < machine_count; ++machine_id) {
const string name = StringPrintf("Machine_%d", machine_id);
SequenceVar* const sequence =
solver.MakeSequenceVar(machines_to_tasks[machine_id], name);
all_sequences.push_back(sequence);
}
// Creates array of end_times of jobs.
std::vector<IntVar*> all_ends;
for (int job_id = 0; job_id < job_count; ++job_id) {
const int task_count = jobs_to_tasks[job_id].size();
IntervalVar* const task = jobs_to_tasks[job_id][task_count - 1];
all_ends.push_back(task->EndExpr()->Var());
}
// Objective: minimize the makespan (maximum end times of all tasks)
// of the problem.
IntVar* const objective_var = solver.MakeMax(all_ends)->Var();
// ----- Search monitors and decision builder -----
// This decision builder will rank all tasks on all machines.
DecisionBuilder* const sequence_phase =
solver.MakePhase(all_sequences, Solver::SEQUENCE_DEFAULT);
// After the ranking of tasks, the schedule is still loose and any
// task can be postponed at will. But, because the problem is now a PERT
// (http://en.wikipedia.org/wiki/Program_Evaluation_and_Review_Technique),
// we can schedule each task at its earliest start time. This is
// conveniently done by fixing the objective variable to its
// minimum value.
DecisionBuilder* const obj_phase =
solver.MakePhase(objective_var,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
Assignment* const first_solution = solver.MakeAssignment();
first_solution->Add(all_sequences);
first_solution->AddObjective(objective_var);
// Store the first solution in the 'solution' object.
DecisionBuilder* const store_db = solver.MakeStoreAssignment(first_solution);
// The main decision builder (ranks all tasks, then fixes the
// objective_variable).
DecisionBuilder* const first_solution_phase =
solver.Compose(sequence_phase, obj_phase, store_db);
LOG(INFO) << "Looking for the first solution";
const bool first_solution_found = solver.Solve(first_solution_phase);
if (first_solution_found) {
LOG(INFO) << "Solution found with makespan = "
<< first_solution->ObjectiveValue();
} else {
LOG(INFO) << "No initial solution found!";
return;
}
LOG(INFO) << "Switching to local search";
std::vector<LocalSearchOperator*> operators;
LOG(INFO) << " - use swap operator";
LocalSearchOperator* const swap_operator =
solver.RevAlloc(new SwapIntervals(all_sequences.data(),
all_sequences.size()));
operators.push_back(swap_operator);
LOG(INFO) << " - use shuffle operator with a max length of "
<< FLAGS_shuffle_length;
LocalSearchOperator* const shuffle_operator =
solver.RevAlloc(new ShuffleIntervals(all_sequences.data(),
all_sequences.size(),
FLAGS_shuffle_length));
operators.push_back(shuffle_operator);
LOG(INFO) << " - use free sub sequences of length "
<< FLAGS_sub_sequence_length << " lns operator";
LocalSearchOperator* const lns_operator =
solver.RevAlloc(new SequenceLns(all_sequences.data(),
all_sequences.size(),
FLAGS_lns_seed,
FLAGS_sub_sequence_length));
operators.push_back(lns_operator);
// Creates the local search decision builder.
LocalSearchOperator* const concat =
solver.ConcatenateOperators(operators, true);
SearchLimit* const ls_limit =
solver.MakeLimit(kint64max, FLAGS_lns_limit, kint64max, kint64max);
DecisionBuilder* const random_sequence_phase =
solver.MakePhase(all_sequences, Solver::CHOOSE_RANDOM_RANK_FORWARD);
DecisionBuilder* const ls_db =
solver.MakeSolveOnce(solver.Compose(random_sequence_phase, obj_phase),
ls_limit);
LocalSearchPhaseParameters* const parameters =
solver.MakeLocalSearchPhaseParameters(concat, ls_db);
DecisionBuilder* const final_db =
solver.MakeLocalSearchPhase(first_solution, parameters);
OptimizeVar* const objective_monitor = solver.MakeMinimize(objective_var, 1);
// Search log.
const int kLogFrequency = 1000000;
SearchMonitor* const search_log =
solver.MakeSearchLog(kLogFrequency, objective_monitor);
SearchLimit* const limit = FLAGS_time_limit_in_ms > 0 ?
solver.MakeTimeLimit(FLAGS_time_limit_in_ms) :
NULL;
// Search.
solver.Solve(final_db, search_log, objective_monitor, limit);
}
} // namespace operations_research
static const char kUsage[] =
"Usage: see flags.\nThis program runs a simple job shop optimization "
"output besides the debug LOGs of the solver.";
int main(int argc, char **argv) {
google::SetUsageMessage(kUsage);
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_data_file.empty()) {
LOG(FATAL) << "Please supply a data file with --data_file=";
}
operations_research::JobShopData data;
data.Load(FLAGS_data_file);
operations_research::JobshopLs(data);
return 0;
}

56
examples/cpp/linear_assignment_api.cc Normal file
View File

@@ -0,0 +1,56 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "graph/ebert_graph.h"
#include "graph/linear_assignment.h"
namespace operations_research {
// Test assignment on a 4x4 matrix. Example taken from
// http://www.ee.oulu.fi/~mpa/matreng/eem1_2-1.htm with kCost[0][1]
// modified so the optimum solution is unique.
void AssignmentOn4x4Matrix() {
LOG(INFO) << "Assignment on 4x4 Matrix";
const int kNumSources = 4;
const int kNumTargets = 4;
const CostValue kCost[kNumSources][kNumTargets] = {
{ 90, 76, 75, 80 },
{ 35, 85, 55, 65 },
{ 125, 95, 90, 105 },
{ 45, 110, 95, 115 }
};
const CostValue kExpectedCost =
kCost[0][3] + kCost[1][2] + kCost[2][1] + kCost[3][0];
ForwardStarGraph graph(
kNumSources + kNumTargets, kNumSources * kNumTargets);
LinearSumAssignment<ForwardStarGraph> assignment(graph, kNumSources);
for (NodeIndex source = 0; source < kNumSources; ++source) {
for (NodeIndex target = 0; target < kNumTargets; ++target) {
ArcIndex arc = graph.AddArc(source, kNumSources + target);
assignment.SetArcCost(arc, kCost[source][target]);
}
}
CHECK(assignment.ComputeAssignment());
CostValue total_cost = assignment.GetCost();
CHECK_EQ(kExpectedCost, total_cost);
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::AssignmentOn4x4Matrix();
return 0;
}

105
examples/cpp/linear_programming.cc Normal file
View File

@@ -0,0 +1,105 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Linear programming example that shows how to use the API.
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "linear_solver/linear_solver.h"
namespace operations_research {
void RunLinearProgrammingExample(
MPSolver::OptimizationProblemType optimization_problem_type) {
MPSolver solver("LinearProgrammingExample", optimization_problem_type);
const double infinity = solver.infinity();
// x1, x2 and x3 are continuous non-negative variables.
MPVariable* const x1 = solver.MakeNumVar(0.0, infinity, "x1");
MPVariable* const x2 = solver.MakeNumVar(0.0, infinity, "x2");
MPVariable* const x3 = solver.MakeNumVar(0.0, infinity, "x3");
// Maximize 10 * x1 + 6 * x2 + 4 * x3.
solver.SetObjectiveCoefficient(x1, 10);
solver.SetObjectiveCoefficient(x2, 6);
solver.SetObjectiveCoefficient(x3, 4);
solver.SetMaximization();
// x1 + x2 + x3 <= 100.
MPConstraint* const c0 = solver.MakeRowConstraint(-infinity, 100.0);
c0->SetCoefficient(x1, 1);
c0->SetCoefficient(x2, 1);
c0->SetCoefficient(x3, 1);
// 10 * x1 + 4 * x2 + 5 * x3 <= 600.
MPConstraint* const c1 = solver.MakeRowConstraint(-infinity, 600.0);
c1->SetCoefficient(x1, 10);
c1->SetCoefficient(x2, 4);
c1->SetCoefficient(x3, 5);
// 2 * x1 + 2 * x2 + 6 * x3 <= 300.
MPConstraint* const c2 = solver.MakeRowConstraint(-infinity, 300.0);
c2->SetCoefficient(x1, 2);
c2->SetCoefficient(x2, 2);
c2->SetCoefficient(x3, 6);
// TODO(user): Change example to show = and >= constraints.
LOG(INFO) << "Number of variables = " << solver.NumVariables();
LOG(INFO) << "Number of constraints = " << solver.NumConstraints();
const MPSolver::ResultStatus result_status = solver.Solve();
// Check that the problem has an optimal solution.
if (result_status != MPSolver::OPTIMAL) {
LOG(FATAL) << "The problem does not have an optimal solution!";
}
LOG(INFO) << "Problem solved in " << solver.wall_time() << " milliseconds";
// The objective value of the solution.
LOG(INFO) << "Optimal objective value = " << solver.objective_value();
// The value of each variable in the solution.
LOG(INFO) << "x1 = " << x1->solution_value();
LOG(INFO) << "x2 = " << x2->solution_value();
LOG(INFO) << "x3 = " << x3->solution_value();
LOG(INFO) << "Advanced usage:";
LOG(INFO) << "Problem solved in " << solver.iterations() << " iterations";
LOG(INFO) << "x1: reduced cost = " << x1->reduced_cost();
LOG(INFO) << "x2: reduced cost = " << x2->reduced_cost();
LOG(INFO) << "x3: reduced cost = " << x3->reduced_cost();
LOG(INFO) << "c0: dual value = " << c0->dual_value()
<< " activity = " << c0->activity();
LOG(INFO) << "c1: dual value = " << c1->dual_value()
<< " activity = " << c1->activity();
LOG(INFO) << "c2: dual value = " << c2->dual_value()
<< " activity = " << c2->activity();
}
void RunAllExamples() {
#if defined(USE_GLPK)
LOG(INFO) << "---- Linear programming example with GLPK ----";
RunLinearProgrammingExample(MPSolver::GLPK_LINEAR_PROGRAMMING);
#endif // USE_GLPK
#if defined(USE_CLP)
LOG(INFO) << "---- Linear programming example with CLP ----";
RunLinearProgrammingExample(MPSolver::CLP_LINEAR_PROGRAMMING);
#endif // USE_CLP
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::RunAllExamples();
return 0;
}

113
examples/cpp/linear_solver_protocol_buffers.cc Normal file
View File

@@ -0,0 +1,113 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include <string>
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "linear_solver/linear_solver.h"
#include "linear_solver/linear_solver.pb.h"
namespace operations_research {
void BuildLinearProgrammingMaxExample(MPSolver::OptimizationProblemType type) {
const double kObjCoef[] = {10.0, 6.0, 4.0};
const string kVarName[] = {"x1", "x2", "x3"};
const int numVars = 3;
const int kNumConstraints = 3;
const string kConstraintName[] = {"c1", "c2", "c3"};
const double kConstraintCoef1[] = {1.0, 1.0, 1.0};
const double kConstraintCoef2[] = {10.0, 4.0, 5.0};
const double kConstraintCoef3[] = {2.0, 2.0, 6.0};
const double* kConstraintCoef[] = {kConstraintCoef1,
kConstraintCoef2,
kConstraintCoef3};
const double kConstraintUb[] = {100.0, 600.0, 300.0};
const double infinity = MPSolver::infinity();
MPModelProto model_proto;
model_proto.set_name("Max_Example");
// Create variables and objective function
for (int j = 0; j < numVars; ++j) {
MPVariableProto* x = model_proto.add_variables();
x->set_id(kVarName[j]);
x->set_lb(0.0);
x->set_ub(infinity);
x->set_integer(false);
MPTermProto* obj_term = model_proto.add_objective_terms();
obj_term->set_variable_id(kVarName[j]);
obj_term->set_coefficient(kObjCoef[j]);
}
model_proto.set_maximize(true);
// Create constraints
for (int i = 0; i < kNumConstraints; ++i) {
MPConstraintProto* constraint_proto = model_proto.add_constraints();
constraint_proto->set_id(kConstraintName[i]);
constraint_proto->set_lb(-infinity);
constraint_proto->set_ub(kConstraintUb[i]);
for (int j = 0; j < numVars; ++j) {
MPTermProto* term = constraint_proto->add_terms();
term->set_variable_id(kVarName[j]);
term->set_coefficient(kConstraintCoef[i][j]);
}
}
MPModelRequest model_request;
model_request.mutable_model()->CopyFrom(model_proto);
#if defined(USE_GLPK)
if (type == MPSolver::GLPK_LINEAR_PROGRAMMING) {
model_request.set_problem_type(MPModelRequest::GLPK_LINEAR_PROGRAMMING);
}
#endif // USE_GLPK
#if defined(USE_CLP)
if (type == MPSolver::CLP_LINEAR_PROGRAMMING) {
model_request.set_problem_type(MPModelRequest::CLP_LINEAR_PROGRAMMING);
}
#endif // USE_CLP
MPSolutionResponse solution_response;
MPSolver::SolveWithProtocolBuffers(model_request, &solution_response);
// The problem has an optimal solution.
CHECK_EQ(MPSolutionResponse::OPTIMAL, solution_response.result_status());
LOG(INFO) << "objective = " << solution_response.objective_value();
const int num_non_zeros = solution_response.solution_values_size();
for (int j = 0; j < num_non_zeros; ++j) {
MPSolutionValue solution_value = solution_response.solution_values(j);
LOG(INFO) << solution_value.variable_id() << " = "
<< solution_value.value();
}
if (num_non_zeros != numVars) {
LOG(INFO) << "All other variables have zero value";
}
}
void RunAllExamples() {
#if defined(USE_GLPK)
LOG(INFO) << "----- Running Max Example with GLPK -----";
BuildLinearProgrammingMaxExample(MPSolver::GLPK_LINEAR_PROGRAMMING);
#endif // USE_GLPK
#if defined(USE_CLP)
LOG(INFO) << "----- Running Max Example with Coin LP -----";
BuildLinearProgrammingMaxExample(MPSolver::CLP_LINEAR_PROGRAMMING);
#endif // USE_CLP
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
operations_research::RunAllExamples();
return 0;
}

163
examples/cpp/magic_square.cc Normal file
View File

@@ -0,0 +1,163 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Magic square problem
//
// Solves the problem where all numbers in an nxn array have to be different
// while the sums on diagonals, rows, and columns have to be the same.
// The problem is trivial for odd orders, but not for even orders.
// We do not handle odd orders with the trivial method here.
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "constraint_solver/constraint_solver.h"
DEFINE_int32(size, 0, "Size of the magic square.");
DEFINE_bool(impact, false, "Use impact search.");
DEFINE_int32(restart, -1, "Parameter for constant restart monitor.");
DEFINE_bool(luby, false,
"Use luby restart monitor instead of constant restart monitor.");
DEFINE_bool(run_all_heuristics, false, "Run all heuristics.");
DEFINE_int32(heuristics_period, 200, "Frequency to run all heuristics.");
DEFINE_int32(choose_var_strategy, 0,
"Selection strategy for variable: 0 = max sum impact, "
"1 = max average impact, "
"2 = max individual impact.");
DEFINE_bool(select_max_impact_value, false,
"Select the value with max impact instead of min impact.");
DEFINE_double(restart_log_size, -1.0,
"Threshold for automatic restarting the search in default"
" phase.");
DEFINE_bool(verbose_impact, false, "Verbose output of impact search.");
DEFINE_bool(use_nogoods, false, "Use no goods in automatic restart.");
namespace operations_research {
void MagicSquare(int grid_size) {
Solver solver("magicsquare");
const int total_size = grid_size * grid_size;
const int sum = grid_size * (total_size + 1) / 2;
// create the variables
std::vector<IntVar*> vars;
solver.MakeIntVarArray(total_size, 1, total_size, "", &vars);
solver.AddConstraint(solver.MakeAllDifferent(vars));
// create the constraints
std::vector<IntVar*> diag1(grid_size);
std::vector<IntVar*> diag2(grid_size);
for (int n = 0; n < grid_size; ++n) {
std::vector<IntVar *> sub_set(grid_size);
for (int m = 0; m < grid_size; ++m) { // extract row indices
sub_set[m] = vars[m + n * grid_size];
}
solver.AddConstraint(solver.MakeSumEquality(sub_set, sum));
for (int m = 0; m < grid_size; ++m) {
sub_set[m] = vars[m * grid_size + n]; // extract column indices
}
solver.AddConstraint(solver.MakeSumEquality(sub_set, sum));
diag1[n] = vars[n + n * grid_size]; // extract first diagonal indices
diag2[n] = vars[(grid_size - 1 - n) + n * grid_size]; // second diagonal
}
solver.AddConstraint(solver.MakeSumEquality(diag1, sum));
solver.AddConstraint(solver.MakeSumEquality(diag2, sum));
// To break a simple symmetry: the upper right corner
// must be less than the lower left corner
solver.AddConstraint(solver.MakeLess(vars[grid_size - 1],
vars[(grid_size - 1) * grid_size]));
// TODO(user) use local search
DefaultPhaseParameters parameters;
parameters.run_all_heuristics = FLAGS_run_all_heuristics;
parameters.heuristic_period = FLAGS_heuristics_period;
parameters.restart_log_size = FLAGS_restart_log_size;
parameters.display_level = FLAGS_verbose_impact ?
DefaultPhaseParameters::VERBOSE :
DefaultPhaseParameters::NORMAL;
parameters.use_no_goods = FLAGS_use_nogoods;
switch (FLAGS_choose_var_strategy) {
case 0: {
parameters.var_selection_schema =
DefaultPhaseParameters::CHOOSE_MAX_SUM_IMPACT;
break;
}
case 1: {
parameters.var_selection_schema =
DefaultPhaseParameters::CHOOSE_MAX_AVERAGE_IMPACT;
break;
}
case 2: {
parameters.var_selection_schema =
DefaultPhaseParameters::CHOOSE_MAX_VALUE_IMPACT;
break;
}
default: {
LOG(FATAL) << "Should not be here";
}
}
parameters.value_selection_schema = FLAGS_select_max_impact_value?
DefaultPhaseParameters::SELECT_MAX_IMPACT:
DefaultPhaseParameters::SELECT_MIN_IMPACT;
DecisionBuilder* const db = FLAGS_impact?
solver.MakeDefaultPhase(vars, parameters):
solver.MakePhase(vars,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
std::vector<SearchMonitor*> monitors;
SearchMonitor* const log = solver.MakeSearchLog(100000);
monitors.push_back(log);
SearchMonitor* const restart = FLAGS_restart != -1?
(FLAGS_luby?
solver.MakeLubyRestart(FLAGS_restart):
solver.MakeConstantRestart(FLAGS_restart)):
NULL;
if (restart) {
monitors.push_back(restart);
}
solver.NewSearch(db, monitors);
if (solver.NextSolution()) {
for (int n = 0; n < grid_size; ++n) {
string output;
for (int m = 0; m < grid_size; ++m) { // extract row indices
int64 v = vars[n * grid_size + m]->Value();
StringAppendF(&output, "%3lld ", v);
}
LG << output;
}
LG << "";
} else {
LG << "No solution found!";
}
solver.EndSearch();
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_size != 0) {
operations_research::MagicSquare(FLAGS_size);
} else {
for (int n = 3; n < 6; ++n) {
operations_research::MagicSquare(n);
}
}
return 0;
}

424
examples/cpp/model_util.cc Normal file
View File

@@ -0,0 +1,424 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/macros.h"
#include "base/file.h"
#include "base/recordio.h"
#include "constraint_solver/constraint_solver.h"
#include "constraint_solver/model.pb.h"
#include "util/graph_export.h"
#include "util/string_array.h"
DEFINE_string(input, "", "Input file of the problem.");
DEFINE_string(output, "", "Output file when doing modifications.");
DEFINE_string(dot_file, "", "Exports model to dot file.");
DEFINE_string(gml_file, "", "Exports model to gml file.");
DEFINE_bool(print_proto, false, "Prints the raw model protobuf.");
DEFINE_bool(test_proto, false, "Performs various tests on the model protobuf.");
DEFINE_bool(model_stats, false, "Prints model statistics.");
DEFINE_bool(print_model, false, "Pretty print loaded model.");
DEFINE_string(rename_model, "", "Renames to the model.");
DEFINE_bool(strip_limit, false, "Strips limits from the model.");
DEFINE_bool(strip_groups, false, "Strips variable groups from the model.");
DEFINE_bool(upgrade_proto, false, "Upgrade the model to the latest version.");
DEFINE_string(insert_license, "",
"Insert content of the given file into the license file.");
DEFINE_bool(collect_variables, false,
"Shows effect of the variable collector.");
namespace operations_research {
// ----- Utilities -----
static const int kProblem = -1;
static const int kOk = 0;
// Colors
static const char kGreen1[] = "#A2CD5A";
static const char kGreen2[] = "#76EEC6";
static const char kGreen3[] = "#00CD00";
static const char kWhite[] = "#FAFAFA";
static const char kBlue[] = "#87CEFA";
static const char kYellow[] = "#FFF68F";
static const char kRed[] = "#A52A2A";
// Creates node labels.
string ExprLabel(int index) {
return StringPrintf("expr_%i", index);
}
string IntervalLabel(int index) {
return StringPrintf("interval_%i", index);
}
string SequenceLabel(int index) {
return StringPrintf("sequence_%i", index);
}
string ConstraintLabel(int index) {
return StringPrintf("ct_%i", index);
}
// Scans argument to add links in the graph.
template <class T> void ExportLinks(const CPModelProto& model,
const string& source,
const T& proto,
GraphExporter* const exporter) {
const string& arg_name = model.tags(proto.argument_index());
if (proto.has_integer_expression_index()) {
exporter->WriteLink(source,
ExprLabel(proto.integer_expression_index()),
arg_name);
}
for (int i = 0; i < proto.integer_expression_array_size(); ++i) {
exporter->WriteLink(source,
ExprLabel(proto.integer_expression_array(i)),
arg_name);
}
if (proto.has_interval_index()) {
exporter->WriteLink(source,
IntervalLabel(proto.interval_index()),
arg_name);
}
for (int i = 0; i < proto.interval_array_size(); ++i) {
exporter->WriteLink(source,
IntervalLabel(proto.interval_array(i)),
arg_name);
}
if (proto.has_sequence_index()) {
exporter->WriteLink(source,
SequenceLabel(proto.sequence_index()),
arg_name);
}
for (int i = 0; i < proto.sequence_array_size(); ++i) {
exporter->WriteLink(source,
SequenceLabel(proto.sequence_array(i)),
arg_name);
}
}
// Scans the expression protobuf to see if it corresponds to an
// integer variable with min_value == max_value.
bool GetValueIfConstant(const CPModelProto& model,
const CPIntegerExpressionProto& proto,
int64* const value) {
CHECK_NOTNULL(value);
const int expr_type = proto.type_index();
if (model.tags(expr_type) != ModelVisitor::kIntegerVariable) {
return false;
}
if (proto.arguments_size() != 2) {
return false;
}
const CPArgumentProto& arg1_proto = proto.arguments(0);
if (model.tags(arg1_proto.argument_index()) != ModelVisitor::kMinArgument) {
return false;
}
const int64 value1 = arg1_proto.integer_value();
const CPArgumentProto& arg2_proto = proto.arguments(1);
if (model.tags(arg2_proto.argument_index()) != ModelVisitor::kMaxArgument) {
return false;
}
const int64 value2 = arg2_proto.integer_value();
if (value1 == value2) {
*value = value1;
return true;
} else {
return false;
}
}
// Declares a labelled expression in the graph file.
void DeclareExpression(int index,
const CPModelProto& proto,
GraphExporter* const exporter) {
const CPIntegerExpressionProto& expr = proto.expressions(index);
const string label = ExprLabel(index);
int64 value = 0;
if (expr.has_name()) {
exporter->WriteNode(label, expr.name(), "oval", kGreen1);
} else if (GetValueIfConstant(proto, expr, &value)) {
exporter->WriteNode(label,
StringPrintf("%lld", value),
"oval",
kYellow);
} else {
const string& type = proto.tags(expr.type_index());
exporter->WriteNode(label, type, "oval", kWhite);
}
}
void DeclareInterval(int index,
const CPModelProto& proto,
GraphExporter* const exporter) {
const CPIntervalVariableProto& interval = proto.intervals(index);
const string label = IntervalLabel(index);
if (interval.has_name()) {
exporter->WriteNode(label, interval.name(), "circle", kGreen2);
} else {
const string& type = proto.tags(interval.type_index());
exporter->WriteNode(label, type, "circle", kWhite);
}
}
void DeclareSequence(int index,
const CPModelProto& proto,
GraphExporter* const exporter) {
const CPSequenceVariableProto& sequence = proto.sequences(index);
const string label = SequenceLabel(index);
if (sequence.has_name()) {
exporter->WriteNode(label, sequence.name(), "circle", kGreen3);
} else {
const string& type = proto.tags(sequence.type_index());
exporter->WriteNode(label, type, "circle", kWhite);
}
}
void DeclareConstraint(int index,
const CPModelProto& proto,
GraphExporter* const exporter) {
const CPConstraintProto& ct = proto.constraints(index);
const string& type = proto.tags(ct.type_index());
const string label = ConstraintLabel(index);
exporter->WriteNode(label, type, "rectangle", kBlue);
}
// Parses the proto and exports it to a graph file.
void ExportToGraphFile(const CPModelProto& proto,
File* const file,
GraphExporter::GraphFormat format) {
scoped_ptr<GraphExporter> exporter(
GraphExporter::MakeFileExporter(file, format));
exporter->WriteHeader(proto.model());
for (int i = 0; i < proto.expressions_size(); ++i) {
DeclareExpression(i, proto, exporter.get());
}
for (int i = 0; i < proto.intervals_size(); ++i) {
DeclareInterval(i, proto, exporter.get());
}
for (int i = 0; i < proto.sequences_size(); ++i) {
DeclareSequence(i, proto, exporter.get());
}
for (int i = 0; i < proto.constraints_size(); ++i) {
DeclareConstraint(i, proto, exporter.get());
}
const char kObjLabel[] = "obj";
if (proto.has_objective()) {
const string name = proto.objective().maximize() ? "Maximize" : "Minimize";
exporter->WriteNode(kObjLabel, name, "diamond", kRed);
}
for (int i = 0; i < proto.expressions_size(); ++i) {
const CPIntegerExpressionProto& expr = proto.expressions(i);
const string label = ExprLabel(i);
for (int j = 0; j < expr.arguments_size(); ++j) {
ExportLinks(proto, label, expr.arguments(j), exporter.get());
}
}
for (int i = 0; i < proto.intervals_size(); ++i) {
const CPIntervalVariableProto& interval = proto.intervals(i);
const string label = IntervalLabel(i);
for (int j = 0; j < interval.arguments_size(); ++j) {
ExportLinks(proto, label, interval.arguments(j), exporter.get());
}
}
for (int i = 0; i < proto.sequences_size(); ++i) {
const CPSequenceVariableProto& sequence = proto.sequences(i);
const string label = SequenceLabel(i);
for (int j = 0; j < sequence.arguments_size(); ++j) {
ExportLinks(proto, label, sequence.arguments(j), exporter.get());
}
}
for (int i = 0; i < proto.constraints_size(); ++i) {
const CPConstraintProto& ct = proto.constraints(i);
const string label = ConstraintLabel(i);
for (int j = 0; j < ct.arguments_size(); ++j) {
ExportLinks(proto, label, ct.arguments(j), exporter.get());
}
}
if (proto.has_objective()) {
const CPObjectiveProto& obj = proto.objective();
exporter->WriteLink(kObjLabel,
ExprLabel(obj.objective_index()),
ModelVisitor::kExpressionArgument);
}
exporter->WriteFooter();
}
// ----- Main Method -----
int Run() {
// ----- Load input file into protobuf -----
File::Init();
File* const file = File::Open(FLAGS_input, "r");
if (file == NULL) {
LOG(WARNING) << "Cannot open " << FLAGS_input;
return kProblem;
}
CPModelProto model_proto;
RecordReader reader(file);
if (!(reader.ReadProtocolMessage(&model_proto) && reader.Close())) {
LOG(INFO) << "No model found in " << file->CreateFileName();
return kProblem;
}
// ----- Display loaded protobuf -----
LOG(INFO) << "Read model " << model_proto.model();
if (model_proto.has_license_text()) {
LOG(INFO) << "License = " << model_proto.license_text();
}
// ----- Modifications -----
if (!FLAGS_rename_model.empty()) {
model_proto.set_model(FLAGS_rename_model);
}
if (FLAGS_strip_limit) {
model_proto.clear_search_limit();
}
if (FLAGS_strip_groups) {
model_proto.clear_variable_groups();
}
if (FLAGS_upgrade_proto) {
if (!Solver::UpgradeModel(&model_proto)) {
LOG(ERROR) << "Model upgrade failed";
return kProblem;
}
}
if (!FLAGS_insert_license.empty()) {
File* const license = File::Open(FLAGS_insert_license, "rb");
if (license == NULL) {
LOG(WARNING) << "Cannot open " << FLAGS_insert_license;
return kProblem;
}
const int size = license->Size();
char* const text = new char[size + 1];
license->Read(text, size);
text[size] = '\0';
model_proto.set_license_text(text);
license->Close();
}
// ----- Reporting -----
if (FLAGS_print_proto) {
LOG(INFO) << model_proto.DebugString();
}
if (FLAGS_test_proto ||
FLAGS_model_stats ||
FLAGS_print_model ||
FLAGS_collect_variables) {
Solver solver(model_proto.model());
std::vector<SearchMonitor*> monitors;
if (!solver.LoadModel(model_proto, &monitors)) {
LOG(INFO) << "Could not load model into the solver";
return kProblem;
}
if (FLAGS_test_proto) {
LOG(INFO) << "Model " << model_proto.model() << " loaded OK";
}
if (FLAGS_model_stats) {
ModelVisitor* const visitor = solver.MakeStatisticsModelVisitor();
solver.Accept(visitor, monitors);
}
if (FLAGS_print_model) {
ModelVisitor* const visitor = solver.MakePrintModelVisitor();
solver.Accept(visitor, monitors);
}
if (FLAGS_collect_variables) {
std::vector<IntVar*> primary_integer_variables;
std::vector<IntVar*> secondary_integer_variables;
std::vector<SequenceVar*> sequence_variables;
std::vector<IntervalVar*> interval_variables;
solver.CollectDecisionVariables(&primary_integer_variables,
&secondary_integer_variables,
&sequence_variables,
&interval_variables);
LOG(INFO) << "Primary integer variables = "
<< DebugStringVector(primary_integer_variables, ", ");
LOG(INFO) << "Secondary integer variables = "
<< DebugStringVector(secondary_integer_variables, ", ");
LOG(INFO) << "Sequence variables = "
<< DebugStringVector(sequence_variables, ", ");
LOG(INFO) << "Interval variables = "
<< DebugStringVector(interval_variables, ", ");
}
}
// ----- Output -----
if (!FLAGS_output.empty()) {
File* const output = File::Open(FLAGS_output, "wb");
if (output == NULL) {
LOG(INFO) << "Cannot open " << FLAGS_output;
return kProblem;
}
RecordWriter writer(output);
if (!(writer.WriteProtocolMessage(model_proto) && writer.Close())) {
return kProblem;
} else {
LOG(INFO) << "Model successfully written to " << FLAGS_output;
}
}
if (!FLAGS_dot_file.empty()) {
File* const dot_file = File::Open(FLAGS_dot_file, "w");
if (dot_file == NULL) {
LOG(INFO) << "Cannot open " << FLAGS_dot_file;
return kProblem;
}
ExportToGraphFile(model_proto, dot_file, GraphExporter::DOT_FORMAT);
dot_file->Close();
}
if (!FLAGS_gml_file.empty()) {
File* const gml_file = File::Open(FLAGS_gml_file, "w");
if (gml_file == NULL) {
LOG(INFO) << "Cannot open " << FLAGS_gml_file;
return kProblem;
}
ExportToGraphFile(model_proto, gml_file, GraphExporter::GML_FORMAT);
gml_file->Close();
}
return kOk;
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_input.empty()) {
LOG(FATAL) << "Filename not specified";
}
return operations_research::Run();
}

333
examples/cpp/multidim_knapsack.cc Normal file
View File

@@ -0,0 +1,333 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// This model implements a multidimensional knapsack problem.
#include <cstdio>
#include <cstdlib>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/strtoint.h"
#include "base/file.h"
#include "base/filelinereader.h"
#include "base/split.h"
#include "constraint_solver/constraint_solver.h"
DEFINE_string(
data_file,
"",
"Required: input file description the muldi-dimensional knapsack problem\n "
"to solve. It supports two file format as described in:\n"
" - http://elib.zib.de/pub/Packages/mp-testdata/ip/sac94-suite/readme\n"
" - http://hces.bus.olemiss.edu/tools.html\n");
DEFINE_int32(time_limit_in_ms, 0, "Time limit in ms, <= 0 means no limit.");
DEFINE_int32(simplex_frequency, 0, "Number of nodes explored between each"
" call to the simplex optimizer.");
DEFINE_bool(display_search_log, true, "Display search log.");
namespace operations_research {
// ----- Data -----
class MultiDimKnapsackData {
public:
MultiDimKnapsackData()
: name_(""), line_read_(0), mode_(0), num_dims_(-1), num_items_(-1),
current_bin_(0), current_item_(0), optimal_value_(0),
problem_type_(-1) {}
void Load(const string& filename) {
FileLineReader reader(filename.c_str());
reader.set_line_callback(NewPermanentCallback(
this,
&MultiDimKnapsackData::ProcessNewLine));
reader.Reload();
if (!reader.loaded_successfully()) {
LOG(ERROR) << "Could not open multi dimensional knapsack file";
}
if (optimal_value_ == 0) {
LOG(INFO) << "Successfully loaded problem " << name_ << " with "
<< items() << " items, " << dims() << " dimensions";
} else {
LOG(INFO) << "Successfully loaded problem " << name_ << " with "
<< items() << " items, " << dims()
<< " dimensions and an optimal value of " << optimal_value_;
}
}
// Number of items of the problem.
int items() const { return num_items_; }
// Number of dimensions of the problem.
int dims() const { return num_dims_; }
// Name of the problem.
const string& name() const { return name_; }
int capacity(int i) const { return dims_[i]; }
int profit(int j) const { return profit_[j]; }
int weight(int i, int j) const { return weight_[i][j]; }
int optimal_value() const { return optimal_value_; }
// Used internally.
void ProcessNewLine(char* const line) {
const char* const kWordDelimiters(" ");
std::vector<string> words;
SplitStringUsing(line, kWordDelimiters, &words);
line_read_++;
if (problem_type_ == -1) {
if (words.size() == 1) {
LOG(INFO) << "New data format";
problem_type_ = 1;
} else if (words.size() == 2) {
LOG(INFO) << "Original data format";
problem_type_ = 0;
}
}
if (problem_type_ == 0) {
// 0 = init
// 1 = size passed
// 2 = profit passed
// 3 = capacity passed
// 4 = constraint passed
// 5 = optimum passed
// 6 = name passed
switch (mode_) {
case 0: {
if (words.size() != 0) {
CHECK_EQ(2, words.size());
num_dims_ = atoi32(words[0]);
num_items_ = atoi32(words[1]);
weight_.resize(num_dims_);
mode_ = 1;
}
break;
}
case 1: {
for (int i = 0; i < words.size(); ++i) {
const int val = atoi32(words[i]);
profit_.push_back(val);
}
CHECK_LE(profit_.size(), num_items_);
if (profit_.size() == num_items_) {
mode_ = 2;
}
break;
}
case 2: {
for (int i = 0; i < words.size(); ++i) {
const int val = atoi32(words[i]);
dims_.push_back(val);
}
CHECK_LE(dims_.size(), num_dims_);
if (dims_.size() == num_dims_) {
mode_ = 3;
}
break;
}
case 3: {
for (int i = 0; i < words.size(); ++i) {
const int val = atoi32(words[i]);
weight_[current_bin_].push_back(val);
if (weight_[current_bin_].size() == num_items_) {
current_bin_++;
}
}
if (current_bin_ == num_dims_) {
mode_ = 4;
}
break;
}
case 4: {
if (words.size() != 0) {
CHECK_EQ(1, words.size());
optimal_value_ = atoi32(words[0]);
mode_ = 5;
}
break;
}
case 5: {
if (words.size() != 0) {
name_ = line;
mode_ = 6;
}
break;
}
case 6: {
break;
}
}
} else {
// 0 = init
// 1 = name passed
// 2 = size passed
// 3 = data passed
// 4 = capacity passed
switch (mode_) {
case 0: {
name_ = words[0];
mode_ = 1;
current_bin_ = -1;
break;
}
case 1: {
if (words.size() != 0) {
CHECK_EQ(2, words.size());
num_items_ = atoi32(words[0]);
num_dims_ = atoi32(words[1]);
weight_.resize(num_dims_);
mode_ = 2;
}
break;
}
case 2: {
for (int i = 0; i < words.size(); ++i) {
const int val = atoi32(words[i]);
if (current_bin_ == -1) {
profit_.push_back(val);
} else {
weight_[current_bin_].push_back(val);
CHECK_EQ(current_item_, weight_[current_bin_].size() - 1);
}
current_bin_++;
if (current_bin_ == num_dims_) {
current_bin_ = -1;
current_item_++;
}
if (current_item_ == num_items_) {
mode_ = 3;
}
}
break;
}
case 3: {
for (int i = 0; i < words.size(); ++i) {
const int val = atoi32(words[i]);
dims_.push_back(val);
}
CHECK_LE(dims_.size(), num_dims_);
if (dims_.size() == num_dims_) {
mode_ = 4;
}
break;
}
case 4:
break;
}
}
}
private:
string name_;
std::vector<int> dims_;
std::vector<int> profit_;
std::vector<std::vector<int> > weight_;
int line_read_;
int mode_;
int num_dims_;
int num_items_;
int current_bin_;
int current_item_;
int optimal_value_;
int problem_type_; // -1 = undefined, 0 = original, 1 = new format
};
int64 EvaluateItem(MultiDimKnapsackData* const data, int64 var, int64 val) {
if (val == 0) {
return 0LL;
}
const int profit = data->profit(var);
int max_weight = 0;
for (int i = 0; i < data->dims(); ++i) {
const int weight = data->weight(i, var);
if (weight > max_weight) {
max_weight = weight;
}
}
return -(profit * 100 / max_weight);
}
void SolveKnapsack(MultiDimKnapsackData* const data) {
Solver solver("MultiDim Knapsack");
std::vector<IntVar*> assign;
solver.MakeBoolVarArray(data->items(), "assign", &assign);
for (int i = 0; i < data->dims(); ++i) {
const int capacity = data->capacity(i);
std::vector<int64> coefs;
for (int j = 0; j < data->items(); ++j) {
coefs.push_back(data->weight(i, j));
}
solver.AddConstraint(
solver.MakeScalProdLessOrEqual(assign, coefs, capacity));
}
// Build objective.
std::vector<int64> profits;
for (int i = 0; i < data->items(); ++i) {
profits.push_back(data->profit(i));
}
IntVar* const objective = solver.MakeScalProd(assign, profits)->Var();
std::vector<SearchMonitor*> monitors;
OptimizeVar* const objective_monitor = solver.MakeMaximize(objective, 1);
monitors.push_back(objective_monitor);
if (FLAGS_display_search_log) {
SearchMonitor* const search_log = solver.MakeSearchLog(1000000, objective);
monitors.push_back(search_log);
}
DecisionBuilder* const db =
solver.MakePhase(assign,
NewPermanentCallback(&EvaluateItem, data),
Solver::CHOOSE_STATIC_GLOBAL_BEST);
if (FLAGS_time_limit_in_ms != 0) {
LOG(INFO) << "adding time limit of " << FLAGS_time_limit_in_ms << " ms";
SearchLimit* const limit = solver.MakeLimit(FLAGS_time_limit_in_ms,
kint64max,
kint64max,
kint64max);
monitors.push_back(limit);
}
if (FLAGS_simplex_frequency > 0) {
SearchMonitor* const simplex =
solver.MakeSimplexConstraint(FLAGS_simplex_frequency);
monitors.push_back(simplex);
}
if (solver.Solve(db, monitors)) {
LOG(INFO) << "Best solution found = " << objective_monitor->best();
}
}
} // namespace operations_research
static const char kUsage[] =
"Usage: see flags.\n"
"This program runs a multi-dimensional knapsack problem.";
int main(int argc, char **argv) {
google::SetUsageMessage(kUsage);
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_data_file.empty()) {
LOG(FATAL) << "Please supply a data file with --datafile=";
}
operations_research::MultiDimKnapsackData data;
data.Load(FLAGS_data_file);
operations_research::SolveKnapsack(&data);
return 0;
}

1053
examples/cpp/network_routing.cc Normal file
View File

File diff suppressed because it is too large Load Diff

275
examples/cpp/nqueens.cc Normal file
View File

@@ -0,0 +1,275 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// N-queens problem
//
// unique solutions: http://www.research.att.com/~njas/sequences/A000170
// distinct solutions: http://www.research.att.com/~njas/sequences/A002562
#include <cstdio>
#include <map>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "base/stringprintf.h"
#include "base/map-util.h"
#include "constraint_solver/constraint_solveri.h"
DEFINE_bool(print, false, "If true, print one of the solution.");
DEFINE_bool(print_all, false, "If true, print all the solutions.");
DEFINE_int32(nb_loops, 1,
"Number of solving loops to perform, for performance timing.");
DEFINE_int32(size, 0,
"Size of the problem. If equal to 0, will test several increasing sizes.");
DEFINE_bool(use_symmetry, false, "Use Symmetry Breaking methods");
DECLARE_bool(cp_no_solve);
static const int kNumSolutions[] = {
1, 0, 0, 2, 10, 4, 40, 92, 352, 724,
2680, 14200, 73712, 365596, 2279184
};
static const int kKnownSolutions = 15;
static const int kNumUniqueSolutions[] = {
1, 0, 0, 1, 2, 1, 6, 12, 46, 92, 341, 1787, 9233, 45752,
285053, 1846955, 11977939, 83263591, 621012754
};
static const int kKnownUniqueSolutions = 19;
namespace operations_research {
class NQueenSymmetry : public SymmetryBreaker {
public:
NQueenSymmetry(Solver* const s, const std::vector<IntVar*>& vars)
: solver_(s), vars_(vars), size_(vars.size()) {
for (int i = 0; i < size_; ++i) {
indices_[vars[i]] = i;
}
}
virtual ~NQueenSymmetry() {}
protected:
int Index(IntVar* const var) const {
return FindWithDefault(indices_, var, -1);
}
IntVar* Var(int index) const {
DCHECK_GE(index, 0);
DCHECK_LT(index, size_);
return vars_[index];
}
int size() const { return size_; }
int symmetric(int index) const { return size_ - 1 - index; }
Solver* const solver() const { return solver_; }
private:
Solver* const solver_;
const std::vector<IntVar*> vars_;
std::map<IntVar*, int> indices_;
const int size_;
};
// Symmetry vertical axis.
class SX : public NQueenSymmetry {
public:
SX(Solver* const s, const std::vector<IntVar*>& vars) : NQueenSymmetry(s, vars) {}
virtual ~SX() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(symmetric(index));
AddIntegerVariableEqualValueClause(other_var, value);
}
};
// Symmetry horizontal axis.
class SY : public NQueenSymmetry {
public:
SY(Solver* const s, const std::vector<IntVar*>& vars) : NQueenSymmetry(s, vars) {}
virtual ~SY() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
AddIntegerVariableEqualValueClause(var, symmetric(value));
}
};
// Symmetry first diagonal axis.
class SD1 : public NQueenSymmetry {
public:
SD1(Solver* const s, const std::vector<IntVar*>& vars) : NQueenSymmetry(s, vars) {}
virtual ~SD1() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(value);
AddIntegerVariableEqualValueClause(other_var, index);
}
};
// Symmetry second diagonal axis.
class SD2 : public NQueenSymmetry {
public:
SD2(Solver* const s, const std::vector<IntVar*>& vars) : NQueenSymmetry(s, vars) {}
virtual ~SD2() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(symmetric(value));
AddIntegerVariableEqualValueClause(other_var, symmetric(index));
}
};
// Rotate 1/4 turn.
class R90 : public NQueenSymmetry {
public:
R90(Solver* const s, const std::vector<IntVar*>& vars) : NQueenSymmetry(s, vars) {}
virtual ~R90() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(value);
AddIntegerVariableEqualValueClause(other_var, symmetric(index));
}
};
// Rotate 1/2 turn.
class R180 : public NQueenSymmetry {
public:
R180(Solver* const s, const std::vector<IntVar*>& vars)
: NQueenSymmetry(s, vars) {}
virtual ~R180() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(symmetric(index));
AddIntegerVariableEqualValueClause(other_var, symmetric(value));
}
};
// Rotate 3/4 turn.
class R270 : public NQueenSymmetry {
public:
R270(Solver* const s, const std::vector<IntVar*>& vars)
: NQueenSymmetry(s, vars) {}
virtual ~R270() {}
virtual void VisitSetVariableValue(IntVar* const var, int64 value) {
const int index = Index(var);
IntVar* const other_var = Var(symmetric(value));
AddIntegerVariableEqualValueClause(other_var, index);
}
};
void CheckNumberOfSolutions(int size, int num_solutions) {
if (FLAGS_use_symmetry) {
if (size - 1 < kKnownUniqueSolutions) {
CHECK_EQ(num_solutions, kNumUniqueSolutions[size - 1]);
} else if (!FLAGS_cp_no_solve) {
CHECK_GT(num_solutions, 0);
}
} else {
if (size - 1 < kKnownSolutions) {
CHECK_EQ(num_solutions, kNumSolutions[size - 1]);
} else if (!FLAGS_cp_no_solve) {
CHECK_GT(num_solutions, 0);
}
}
}
void NQueens(int size) {
CHECK_GE(size, 1);
Solver s("nqueens");
// model
std::vector<IntVar*> queens;
for (int i = 0; i < size; ++i) {
queens.push_back(s.MakeIntVar(0, size - 1, StringPrintf("queen%04d", i)));
}
s.AddConstraint(s.MakeAllDifferent(queens));
std::vector<IntVar*> vars(size);
for (int i = 0; i < size; ++i) {
vars[i] = s.MakeSum(queens[i], i)->Var();
}
s.AddConstraint(s.MakeAllDifferent(vars));
for (int i = 0; i < size; ++i) {
vars[i] = s.MakeSum(queens[i], -i)->Var();
}
s.AddConstraint(s.MakeAllDifferent(vars));
SolutionCollector* const solution_counter = s.MakeAllSolutionCollector(NULL);
SolutionCollector* const collector = s.MakeFirstSolutionCollector();
collector->Add(queens);
std::vector<SearchMonitor*> monitors;
monitors.push_back(solution_counter);
monitors.push_back(collector);
DecisionBuilder* const db = s.MakePhase(queens,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
if (FLAGS_use_symmetry) {
std::vector<SymmetryBreaker*> breakers;
NQueenSymmetry* const sx = s.RevAlloc(new SX(&s, queens));
breakers.push_back(sx);
NQueenSymmetry* const sy = s.RevAlloc(new SY(&s, queens));
breakers.push_back(sy);
NQueenSymmetry* const sd1 = s.RevAlloc(new SD1(&s, queens));
breakers.push_back(sd1);
NQueenSymmetry* const sd2 = s.RevAlloc(new SD2(&s, queens));
breakers.push_back(sd2);
NQueenSymmetry* const r90 = s.RevAlloc(new R90(&s, queens));
breakers.push_back(r90);
NQueenSymmetry* const r180 = s.RevAlloc(new R180(&s, queens));
breakers.push_back(r180);
NQueenSymmetry* const r270 = s.RevAlloc(new R270(&s, queens));
breakers.push_back(r270);
SearchMonitor* const symmetry_manager = s.MakeSymmetryManager(breakers);
monitors.push_back(symmetry_manager);
}
for (int loop = 0; loop < FLAGS_nb_loops; ++loop) {
s.Solve(db, monitors); // go!
CheckNumberOfSolutions(size, solution_counter->solution_count());
}
const int num_solutions = solution_counter->solution_count();
if (num_solutions > 0 && size < kKnownSolutions) {
int print_max = FLAGS_print_all ? num_solutions : FLAGS_print ? 1 : 0;
for (int n = 0; n < print_max; ++n) {
printf("--- solution #%d\n", n);
for (int i = 0; i < size; ++i) {
const int pos = static_cast<int>(collector->Value(n, queens[i]));
for (int k = 0; k < pos; ++k) printf(" . ");
printf("%2d ", i);
for (int k = pos + 1; k < size; ++k) printf(" . ");
printf("\n");
}
}
}
printf("========= number of solutions:%d\n", num_solutions);
printf(" number of failures: %lld\n", s.failures());
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_size != 0) {
operations_research::NQueens(FLAGS_size);
} else {
for (int n = 1; n < 12; ++n) {
operations_research::NQueens(n);
}
}
return 0;
}

252
examples/cpp/parse_dimacs_assignment.cc Normal file
View File

@@ -0,0 +1,252 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "cpp/parse_dimacs_assignment.h"
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <string>
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "graph/ebert_graph.h"
#include "graph/linear_assignment.h"
DEFINE_bool(assignment_maximize_cost, false,
"Negate costs so a max-cost assignment is found.");
DEFINE_bool(assignment_optimize_layout, true,
"Optimize graph layout for speed.");
namespace operations_research {
struct ParserState {
ParserState()
: bad(false),
expect_last_line(false),
nodes_described(false),
reason(NULL),
num_left_nodes(0) { }
bool bad;
bool expect_last_line;
bool nodes_described;
const char* reason;
NodeIndex num_left_nodes;
scoped_ptr<string> bad_line;
};
static void ParseProblemLine(const char* line,
ParserState* state,
ForwardStarGraph** graph) {
static const char* kIncorrectProblemLine =
"Incorrect assignment problem line.";
static const char* kAssignmentProblemType = "asn";
char problem_type[4];
NodeIndex num_nodes;
ArcIndex num_arcs;
if ((sscanf(line, "%*c%3s%d%d", // NOLINT
problem_type,
&num_nodes,
&num_arcs) != 3) ||
(strncmp(kAssignmentProblemType,
problem_type,
strlen(kAssignmentProblemType)) != 0)) {
state->bad = true;
state->reason = kIncorrectProblemLine;
state->bad_line.reset(new string(line));
return;
}
*graph = new ForwardStarGraph(num_nodes, num_arcs);
}
static void ParseNodeLine(const char* line,
ParserState* state) {
NodeIndex node_id;
if (sscanf(line, "%*c%d", &node_id) != 1) { // NOLINT
state->bad = true;
state->reason = "Syntax error in node desciption.";
state->bad_line.reset(new string(line));
return;
}
if (state->nodes_described) {
state->bad = true;
state->reason = "All node description must precede first arc description.";
state->bad_line.reset(new string(line));
return;
}
state->num_left_nodes = ::std::max(state->num_left_nodes, node_id);
}
static void ParseArcLine(
const char* line,
ParserState* state,
ForwardStarGraph* graph,
LinearSumAssignment<ForwardStarGraph>** assignment) {
if (graph == NULL) {
state->bad = true;
state->reason =
"Problem specification line must precede any arc specification.";
state->bad_line.reset(new string(line));
return;
}
if (!state->nodes_described) {
state->nodes_described = true;
DCHECK(*assignment == NULL);
*assignment = new LinearSumAssignment<ForwardStarGraph>(
*graph, state->num_left_nodes);
}
NodeIndex tail;
NodeIndex head;
CostValue cost;
if (sscanf(line, "%*c%d%d%lld", &tail, &head, &cost) != 3) { // NOLINT
state->bad = true;
state->reason = "Syntax error in arc descriptor.";
state->bad_line.reset(new string(line));
}
ArcIndex arc = graph->AddArc(tail - 1, head - 1);
(*assignment)->SetArcCost(arc, FLAGS_assignment_maximize_cost ? -cost : cost);
}
// Parameters out of style-guide order because this function is used
// as a callback that varies the input line.
static void ParseOneLine(
ParserState* state,
ForwardStarGraph** graph,
LinearSumAssignment<ForwardStarGraph>** assignment,
char* line) {
if (state->bad) {
return;
}
if (state->expect_last_line) {
state->bad = true;
state->reason = "Input line is too long.";
// state->bad_line was already set when we noticed the line
// didn't end with '\n'.
return;
}
size_t length = strlen(line);
// The final line might not end with newline. Any other line
// that seems not to is actually a line that was too long
// for our input buffer.
if (line[length - 1] != '\n') {
state->expect_last_line = true;
// Prepare for the worst; we might need to inform the user of
// an error on this line even though we can't detect the error
// yet.
state->bad_line.reset(new string(line));
}
switch (line[0]) {
case 'p': {
// Problem-specification line
ParseProblemLine(line, state, graph);
break;
}
case 'c': {
// Comment; do nothing.
return;
}
case 'n': {
// Node line defining a node on the left side
ParseNodeLine(line, state);
break;
}
case 'a': {
ParseArcLine(line, state, *graph, assignment);
break;
}
case '0':
case '\n':
break;
default: {
state->bad = true;
state->reason = "Unknown line type in the input.";
state->bad_line.reset(new string(line));
break;
}
}
}
void ParseFileByLines(const string& filename,
Callback1<char*>* line_parser) {
FILE* fp = fopen(filename.c_str(), "r");
const int kMaximumLineSize = 1024;
char line[kMaximumLineSize];
if (fp != NULL) {
char* result;
do {
result = fgets(line, kMaximumLineSize, fp);
if (result != NULL) {
line_parser->Run(result);
}
} while (result != NULL);
}
delete line_parser;
}
// Reads an assignment problem description from the given file in
// DIMACS format and returns a LinearSumAssignment object representing
// the problem description. For a description of the format, see
// http://lpsolve.sourceforge.net/5.5/DIMACS_asn.htm
//
// Also returns an error message (empty if no error) and a handle on
// the underlying graph representation. The error_message pointer must
// not be NULL because we insist on returning an explanatory message
// in the case of error. The graph_handle pointer must not be NULL
// because unless we pass a non-const pointer to the graph
// representation back to the caller, the caller lacks a good way to
// free the underlying graph (which isn't owned by the
// LinearAssignment instance).
LinearSumAssignment<ForwardStarGraph>* ParseDimacsAssignment(
const string& filename,
string* error_message,
ForwardStarGraph** graph_handle) {
CHECK_NOTNULL(error_message);
CHECK_NOTNULL(graph_handle);
ForwardStarGraph* graph = NULL;
LinearSumAssignment<ForwardStarGraph>* assignment = NULL;
ParserState state;
Callback1<char*>* cb =
NewPermanentCallback(ParseOneLine, &state, &graph, &assignment);
// ParseFileByLines takes ownership of cb and deletes it.
ParseFileByLines(filename, cb);
if (state.bad) {
*error_message = state.reason;
*error_message = *error_message + ": \"" + *state.bad_line + "\"";
return NULL;
}
if (graph == NULL) {
*error_message = "empty graph description";
return NULL;
}
if (FLAGS_assignment_optimize_layout) {
assignment->OptimizeGraphLayout(graph);
}
*error_message = "";
// Return a handle on the graph to the caller so the caller can free
// the graph's memory, because the LinearSumAssignment object does
// not take ownership of the graph and hence will not free it.
*graph_handle = graph;
return assignment;
}
} // namespace operations_research

38
examples/cpp/parse_dimacs_assignment.h Normal file
View File

@@ -0,0 +1,38 @@
// Copyright 2011 Google Inc. All Rights Reserved.
//
// Function for reading and parsing a file in DIMACS format:
// http://lpsolve.sourceforge.net/5.5/DIMACS_asn.htm
//
#ifndef OR_TOOLS_EXAMPLES_PARSE_DIMACS_ASSIGNMENT_H_
#define OR_TOOLS_EXAMPLES_PARSE_DIMACS_ASSIGNMENT_H_
#include <string>
#include "graph/ebert_graph.h"
namespace operations_research {
template <typename GraphType> class LinearSumAssignment;
// Reads an assignment problem description from the given file in
// DIMACS format and returns a LinearSumAssignment object representing
// the problem description. For a description of the format, see
// http://lpsolve.sourceforge.net/5.5/DIMACS_asn.htm
//
// Also returns an error message (empty if no error) and a handle on
// the underlying graph representation. The error_message pointer must
// not be NULL because we insist on returning an explanatory message
// in the case of error. The graph_handle pointer must not be NULL
// because unless we pass a non-const pointer to the graph
// representation back to the caller, the caller lacks a good way to
// free the underlying graph (which isn't owned by the
// LinearAssignment instance).
LinearSumAssignment<ForwardStarGraph>* ParseDimacsAssignment(
const string& filename,
string* error_message,
ForwardStarGraph** graph);
} // namespace operations_research
#endif // OR_TOOLS_EXAMPLES_PARSE_DIMACS_ASSIGNMENT_H_

321
examples/cpp/pdptw.cc Normal file
View File

@@ -0,0 +1,321 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Pickup and Delivery Problem with Time Windows.
// The overall objective is to minimize the length of the routes delivering
// quantities of goods between pickup and delivery locations, taking into
// account vehicle capacities and node time windows.
// Given a set of pairs of pickup and delivery nodes, find the set of routes
// visiting all the nodes, such that
// - corresponding pickup and delivery nodes are visited on the same route,
// - the pickup node is visited before the corresponding delivery node,
// - the quantity picked up at the pickup node is the same as the quantity
// delivered at the delivery node,
// - the total quantity carried by a vehicle at any time is less than its
// capacity,
// - each node must be visited within its time window (time range during which
// the node is accessible).
// The maximum number of vehicles used (i.e. the number of routes used) is
// specified in the data but can be overriden using the --pdp_force_vehicles
// flag.
//
// A further description of the problem can be found here:
// http://http://en.wikipedia.org/wiki/Vehicle_routing_problem
// http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.123.9965&rep=rep1&type=pdf.
// Reads data in the format defined by Li & Lim
// (http://www.sintef.no/static/am/opti/projects/top/vrp/format_pdp.htm).
#include <vector>
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/strtoint.h"
#include "base/file.h"
#include "base/split.h"
#include "base/mathutil.h"
#include "constraint_solver/routing.h"
DEFINE_string(pdp_file, "",
"File containing the Pickup and Delivery Problem to solve.");
DEFINE_int32(pdp_force_vehicles, 0,
"Force the number of vehicles used (maximum number of routes.");
DEFINE_bool(pdp_display_solution, false,
"Displays the solution of the Pickup and Delivery Problem.");
namespace operations_research {
// Scaling factor used to scale up distances, allowing a bit more precision
// from Euclidean distances.
const int64 kScalingFactor = 1000;
// Vector of (x,y) node coordinates, *unscaled*, in some imaginary planar,
// metric grid.
typedef std::vector<std::pair<int, int> > Coordinates;
// Returns the scaled Euclidean distance between two nodes, coords holding the
// coordinates of the nodes.
int64 Travel(const Coordinates* const coords,
RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) {
DCHECK(coords != NULL);
const int xd =
coords->at(from.value()).first - coords->at(to.value()).first;
const int yd =
coords->at(from.value()).second - coords->at(to.value()).second;
return static_cast<int64>(kScalingFactor * sqrt(1.0L * xd * xd + yd * yd));
}
// Returns the scaled service time at a given node, service_times holding the
// service times.
int64 ServiceTime(const std::vector<int64>* const service_times,
RoutingModel::NodeIndex node) {
return kScalingFactor * service_times->at(node.value());
}
// Returns the scaled (distance plus service time) between two nodes, coords
// holding the coordinates of the nodes and service_times holding the service
// times.
// The service time is the time spent to execute a delivery or a pickup.
int64 TravelPlusServiceTime(const Coordinates* const coords,
const std::vector<int64>* const service_times,
RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) {
return ServiceTime(service_times, from) + Travel(coords, from, to);
}
// Returns the demand (quantity picked up or delivered) of a node, demands
// holds the demand of each node.
int64 Demand(const std::vector<int64>* const demands,
RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) {
return demands->at(from.value());
}
// Outputs a solution to the current model in a string.
string VerboseOutput(const RoutingModel& routing,
const Assignment& assignment,
const Coordinates& coords,
const std::vector<int64>& service_times) {
string output;
for (int i = 0; i < routing.vehicles(); ++i) {
StringAppendF(&output, "Vehicle %d: ", i);
int64 index = routing.Start(i);
if (routing.IsEnd(assignment.Value(routing.NextVar(index)))) {
StringAppendF(&output, "empty");
} else {
while (!routing.IsEnd(index)) {
RoutingModel::NodeIndex real_node = routing.IndexToNode(index);
StringAppendF(&output, "%d ", real_node.value());
const IntVar* vehicle = routing.VehicleVar(index);
StringAppendF(&output, "Vehicle(%lld) ", assignment.Value(vehicle));
const IntVar* arrival = routing.CumulVar(index, "time");
StringAppendF(&output, "Time(%lld..%lld) ",
assignment.Min(arrival),
assignment.Max(arrival));
const IntVar* load = routing.CumulVar(index, "demand");
StringAppendF(&output, "Load(%lld..%lld) ",
assignment.Min(load),
assignment.Max(load));
index = assignment.Value(routing.NextVar(index));
StringAppendF(&output, "Transit(%lld) ",
TravelPlusServiceTime(&coords,
&service_times,
real_node,
routing.IndexToNode(index)));
}
StringAppendF(&output, "Route end ");
const IntVar* vehicle = routing.VehicleVar(index);
StringAppendF(&output, "Vehicle(%lld) ", assignment.Value(vehicle));
const IntVar* arrival = routing.CumulVar(index, "time");
StringAppendF(&output, "Time(%lld..%lld) ",
assignment.Min(arrival),
assignment.Max(arrival));
const IntVar* load = routing.CumulVar(index, "demand");
StringAppendF(&output, "Load(%lld..%lld) ",
assignment.Min(load),
assignment.Max(load));
}
StringAppendF(&output, "\n");
}
return output;
}
namespace {
// An inefficient but convenient method to parse a whitespace-separated list
// of integers. Returns true iff the input string was entirely valid and parsed.
bool SafeParseInt64Array(const string& str, std::vector<int64>* parsed_int) {
std::vector<string> items;
static const char kWhiteSpaces[] = " \t\n\v\f\r";
SplitStringUsing(str, kWhiteSpaces, &items);
parsed_int->assign(items.size(), 0);
for (int i = 0; i < items.size(); ++i) {
const char* item = items[i].c_str();
char *endptr = NULL;
(*parsed_int)[i] = strto64(item, &endptr, 10); // NOLINT
// The whole item should have been consumed.
if (*endptr != '\0') return false;
}
return true;
}
}
// Builds and solves a model from a file in the format defined by Li & Lim
// (http://www.sintef.no/static/am/opti/projects/top/vrp/format_pdp.htm).
bool LoadAndSolve(const string& pdp_file) {
// Load all the lines of the file in RAM (it shouldn't be too large anyway).
std::vector<string> lines;
{
const int64 kMaxInputFileSize = 1 << 30; // 1GB
File* data_file = File::OpenOrDie(pdp_file.c_str(), "r");
string contents;
data_file->ReadToString(&contents, kMaxInputFileSize);
data_file->Close();
if (contents.size() == kMaxInputFileSize) {
LOG(WARNING)
<< "Input file '" << pdp_file << "' is too large (>"
<< kMaxInputFileSize << " bytes).";
return false;
}
SplitStringUsing(contents, "\n", &lines);
}
// Reading header.
if (lines.empty()) {
LOG(WARNING) << "Empty file: " << pdp_file;
return false;
}
// Parse file header.
std::vector<int64> parsed_int;
if (!SafeParseInt64Array(lines[0], &parsed_int)
|| parsed_int.size() != 3
|| parsed_int[0] < 0
|| parsed_int[1] < 0
|| parsed_int[2] < 0) {
LOG(WARNING) << "Malformed header: " << lines[0];
return false;
}
const int num_vehicles = FLAGS_pdp_force_vehicles > 0 ?
FLAGS_pdp_force_vehicles : parsed_int[0];
const int64 capacity = parsed_int[1];
// We do not care about the 'speed' field, in third position.
// Parse order data.
std::vector<int> customer_ids;
std::vector<std::pair<int, int> > coords;
std::vector<int64> demands;
std::vector<int64> open_times;
std::vector<int64> close_times;
std::vector<int64> service_times;
std::vector<RoutingModel::NodeIndex> pickups;
std::vector<RoutingModel::NodeIndex> deliveries;
int64 horizon = 0;
for (int line_index = 1; line_index < lines.size(); ++line_index) {
if (!SafeParseInt64Array(lines[line_index], &parsed_int)
|| parsed_int.size() != 9
|| parsed_int[0] < 0
|| parsed_int[4] < 0
|| parsed_int[5] < 0
|| parsed_int[6] < 0
|| parsed_int[7] < 0
|| parsed_int[8] < 0) {
LOG(WARNING)
<< "Malformed line #" << line_index << ": " << lines[line_index];
return false;
}
const int customer_id = parsed_int[0];
const int x = parsed_int[1];
const int y = parsed_int[2];
const int delivery = parsed_int[8]; // Parse 'delivery' before 'demand'.
const int64 demand = delivery == 0 ? -parsed_int[3] : parsed_int[3];
const int64 open_time = parsed_int[4];
const int64 close_time = parsed_int[5];
const int64 service_time = parsed_int[6];
const int pickup = parsed_int[7];
customer_ids.push_back(customer_id);
coords.push_back(std::make_pair(x, y));
demands.push_back(demand);
open_times.push_back(open_time);
close_times.push_back(close_time);
service_times.push_back(service_time);
pickups.push_back(RoutingModel::NodeIndex(pickup));
deliveries.push_back(RoutingModel::NodeIndex(delivery));
horizon = std::max(horizon, close_time);
}
// Build pickup and delivery model.
const int num_nodes = customer_ids.size();
RoutingModel routing(num_nodes, num_vehicles);
routing.SetCost(NewPermanentCallback(
Travel, const_cast<const Coordinates*>(&coords)));
routing.AddDimension(
NewPermanentCallback(&Demand, const_cast<const std::vector<int64>*>(&demands)),
0, capacity, "demand");
routing.AddDimension(
NewPermanentCallback(&TravelPlusServiceTime,
const_cast<const Coordinates*>(&coords),
const_cast<const std::vector<int64>*>(&service_times)),
kScalingFactor * horizon,
kScalingFactor * horizon,
"time");
Solver* const solver = routing.solver();
for (RoutingModel::NodeIndex i(0); i < num_nodes; ++i) {
const int64 index = routing.NodeToIndex(i);
if (pickups[i.value()] == 0) {
if (deliveries[i.value()] == 0) {
routing.SetDepot(i);
} else {
const int64 delivery_index = routing.NodeToIndex(deliveries[i.value()]);
solver->AddConstraint(solver->MakeEquality(
routing.VehicleVar(index),
routing.VehicleVar(delivery_index)));
solver->AddConstraint(solver->MakeLessOrEqual(
routing.CumulVar(index, "time"),
routing.CumulVar(delivery_index, "time")));
routing.AddPickupAndDelivery(i, deliveries[i.value()]);
}
}
IntVar* const cumul = routing.CumulVar(index, "time");
cumul->SetMin(kScalingFactor * open_times[i.value()]);
cumul->SetMax(kScalingFactor * close_times[i.value()]);
}
// Adding penalty costs to allow skipping orders.
const int64 kPenalty = 10000000;
for (RoutingModel::NodeIndex order(1);
order < routing.nodes(); ++order) {
std::vector<RoutingModel::NodeIndex> orders(1, order);
routing.AddDisjunction(orders, kPenalty);
}
// Set up search parameters.
routing.set_first_solution_strategy(RoutingModel::ROUTING_ALL_UNPERFORMED);
routing.SetCommandLineOption("routing_no_lns", "true");
// Solve pickup and delivery problem.
const Assignment* assignment = routing.Solve(NULL);
if (NULL != assignment) {
LG << "Cost: " << assignment->ObjectiveValue();
LG << VerboseOutput(routing, *assignment, coords, service_times);
return true;
}
return false;
}
} // namespace operations_research
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (!operations_research::LoadAndSolve(FLAGS_pdp_file)) {
LG << "Error solving " << FLAGS_pdp_file;
}
return 0;
}

74
examples/cpp/print_dimacs_assignment.cc Normal file
View File

@@ -0,0 +1,74 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "cpp/print_dimacs_assignment.h"
#include <cstdio>
#include <string>
#include "base/logging.h"
#include "base/stringprintf.h"
#include "graph/ebert_graph.h"
#include "graph/linear_assignment.h"
namespace operations_research {
static void WriteOrDie(const char* buffer,
size_t item_size,
size_t buffer_length,
FILE* fp) {
size_t written = fwrite(buffer, item_size, buffer_length, fp);
if (written != buffer_length) {
fprintf(stderr, "Write failed.\n");
exit(1);
}
}
void PrintDimacsAssignmentProblem(
const LinearSumAssignment<ForwardStarGraph>& assignment,
const TailArrayManager<ForwardStarGraph>& tail_array_manager,
const string& output_filename) {
FILE* output = fopen(output_filename.c_str(), "w");
const ForwardStarGraph& graph(assignment.Graph());
string output_line = StringPrintf("p asn %d %d\n",
graph.num_nodes(),
graph.num_arcs());
WriteOrDie(output_line.c_str(), 1, output_line.length(),
output);
for (LinearSumAssignment<ForwardStarGraph>::BipartiteLeftNodeIterator
node_it(assignment);
node_it.Ok();
node_it.Next()) {
output_line = StringPrintf("n %d\n", node_it.Index() + 1);
WriteOrDie(output_line.c_str(), 1, output_line.length(),
output);
}
tail_array_manager.BuildTailArrayFromAdjacencyListsIfForwardGraph();
for (ForwardStarGraph::ArcIterator arc_it(assignment.Graph());
arc_it.Ok();
arc_it.Next()) {
ArcIndex arc = arc_it.Index();
output_line = StringPrintf("a %d %d %lld\n",
graph.Tail(arc) + 1,
graph.Head(arc) + 1,
assignment.ArcCost(arc));
WriteOrDie(output_line.c_str(), 1, output_line.length(),
output);
}
}
} // namespace operations_research

28
examples/cpp/print_dimacs_assignment.h Normal file
View File

@@ -0,0 +1,28 @@
// Copyright 2011 Google Inc. All Rights Reserved.
//
// Function for outputting an assignment problem in DIMACS format:
// http://lpsolve.sourceforge.net/5.5/DIMACS_asn.htm
//
#ifndef OR_TOOLS_EXAMPLES_PRINT_DIMACS_ASSIGNMENT_H_
#define OR_TOOLS_EXAMPLES_PRINT_DIMACS_ASSIGNMENT_H_
#include <string>
#include "graph/ebert_graph.h"
namespace operations_research {
template <typename GraphType> class LinearSumAssignment;
// Given a LinearSumAssigment object representing an assignment problem
// description, outputs the problem in DIMACS format in the output file.
// For a description of the format, see
// http://lpsolve.sourceforge.net/5.5/DIMACS_asn.htm
void PrintDimacsAssignmentProblem(
const LinearSumAssignment<ForwardStarGraph>& assignment,
const TailArrayManager<ForwardStarGraph>& tail_array_manager,
const string& output_filename);
} // namespace operations_research
#endif // OR_TOOLS_EXAMPLES_PRINT_DIMACS_ASSIGNMENT_H_

427
examples/cpp/sports_scheduling.cc Normal file
View File

@@ -0,0 +1,427 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// Sports scheduling problem.
//
// We want to solve the problem of scheduling of team matches in a
// double round robin tournament. Given a number of teams, we want
// each team to encounter all other teams, twice, once at home, and
// once away. Furthermore, you cannot meet the same team twice in the
// same half-season.
//
// Finally, there are constraints on the sequence of home or aways:
// - You cannot have 3 consecutive homes or three consecutive aways.
// - A break is a sequence of two homes or two aways, the overall objective
// of the optimization problem is to minimize the total number of breaks.
//
// We model this problem with three matrices of variables, each with
// num_teams rows and 2*(num_teams - 1) columns: the var [i][j]
// corresponds to the match of team #i at day #j. There are
// 2*(num_teams - 1) columns because each team meets num_teams - 1
// opponents twice.
// - The 'opponent' var [i][j] is the index of the opposing team.
// - The 'home_away' var [i][j] is a boolean: 1 for 'playing away',
// 0 for 'playing at home'.
// - The 'opponent_and_home_away' var [i][j] is the 'opponent' var [i][j] +
// num_teams * the 'home_away' var [i][j].
// This aggregated variable will be useful to state constraints of the model
// and to do search on it.
//
// We use an original approch in this model as most of the constraints will
// be pre-computed and asserted using an AllowedAssignment constraint (see
// Solver::MakeAllowedAssignment() in constraint_solver.h).
// In particular:
// - Each day, we have a perfect matching between teams
// (A meets B <=> B meets A, and A is at home <=> B is away).
// A cannot meet itself.
// - For each team, over the length of the tournament, we have constraints
// on the sequence of home-aways. We will precompute all possible sequences
// of home_aways, as well as the corresponding number of breaks for that
// team.
// - For a given team and a given day, the link between the opponent var,
// the home_away var and the aggregated var (see third matrix of variables)
// is also maintained using a AllowedAssignment constraint.
//
// Usage: run this with --helpshort for a short usage manual.
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "constraint_solver/constraint_solver.h"
#include "constraint_solver/constraint_solveri.h"
// Problem main flags.
DEFINE_int32(num_teams, 10, "Number of teams in the problem.");
// General solving parameters.
DEFINE_int32(time_limit, 20000, "Time limit in ms.");
// Search tweaking parameters. These are defined to illustrate their effect.
DEFINE_bool(run_all_heuristics,
true,
"Run all heuristics in impact search, see DefaultPhaseParameters"
" in constraint_solver/constraint_solver.h for details.");
DEFINE_int32(heuristics_period,
30,
"Frequency to run all heuristics, see DefaultPhaseParameters"
" in constraint_solver/constraint_solver.h for details.");
DEFINE_double(restart_log_size,
8.0,
"Threshold for automatic restarting the search in default phase,"
" see DefaultPhaseParameters in "
"constraint_solver/constraint_solver.h for details.");
namespace operations_research {
// ---------- Utility functions to help create the model ----------
// ----- Constraints for one day and one team -----
// Computes the tuple set that links opponents, home_away, and
// signed_opponent on a single day for a single team.
void ComputeOneDayOneTeamTuples(int num_teams, IntTupleSet* const tuples) {
for (int home_away = 0; home_away <= 1; ++home_away) {
for (int opponent = 0; opponent < num_teams; ++opponent) {
tuples->Insert3(opponent, home_away, opponent + home_away * num_teams);
}
}
}
void AddOneDayOneTeamConstraint(Solver* const solver,
IntVar* const opponent,
IntVar* const home_away,
IntVar* const signed_opponent,
const IntTupleSet& intra_day_tuples) {
std::vector<IntVar*> tmp_vars;
tmp_vars.push_back(opponent);
tmp_vars.push_back(home_away);
tmp_vars.push_back(signed_opponent);
solver->AddConstraint(
solver->MakeAllowedAssignments(tmp_vars, intra_day_tuples));
}
// ----- Constraints for one day and all teams -----
// Computes all valid combination of signed_opponent for a single
// day and all teams.
void ComputeOneDayTuples(int num_teams, IntTupleSet* const day_tuples) {
LOG(INFO) << "Compute possible opponents and locations for any day.";
Solver solver("ComputeOneDayTuples");
// We create the variables.
std::vector<IntVar*> opponents;
std::vector<IntVar*> home_aways;
std::vector<IntVar*> signed_opponents;
solver.MakeIntVarArray(num_teams, 0, num_teams - 1, "opponent_", &opponents);
solver.MakeBoolVarArray(num_teams, "home_away_", &home_aways);
solver.MakeIntVarArray(num_teams,
0,
2 * num_teams - 1,
"signed_opponent_",
&signed_opponents);
// All Diff constraint.
solver.AddConstraint(solver.MakeAllDifferent(opponents));
// Cannot play against itself
for (int i = 0; i < num_teams; ++i) {
solver.AddConstraint(solver.MakeNonEquality(opponents[i], i));
}
// Matching constraint (vars[i] == j <=> vars[j] == i).
for (int i = 0; i < num_teams; ++i) {
for (int j = 0; j < num_teams; ++j) {
if (i != j) {
solver.AddConstraint(
solver.MakeEquality(solver.MakeIsEqualCstVar(opponents[i], j),
solver.MakeIsEqualCstVar(opponents[j], i)));
}
}
}
// num_teams/2 teams are home.
solver.AddConstraint(solver.MakeSumEquality(home_aways, num_teams / 2));
// Link signed_opponents, home_away and opponents
IntTupleSet one_day_one_team_tuples(3);
ComputeOneDayOneTeamTuples(num_teams, &one_day_one_team_tuples);
for (int team_index = 0; team_index < num_teams; ++team_index) {
std::vector<IntVar*> tmp_vars;
tmp_vars.push_back(opponents[team_index]);
tmp_vars.push_back(home_aways[team_index]);
tmp_vars.push_back(signed_opponents[team_index]);
solver.AddConstraint(
solver.MakeAllowedAssignments(tmp_vars, one_day_one_team_tuples));
}
// if A meets B at home, B meets A away.
for (int first_team = 0; first_team < num_teams; ++first_team) {
IntVar* const first_home_away = home_aways[first_team];
IntVar* const second_home_away =
solver.MakeElement(home_aways, opponents[first_team])->Var();
IntVar* const reverse_second_home_away =
solver.MakeDifference(1, second_home_away)->Var();
solver.AddConstraint(
solver.MakeEquality(first_home_away, reverse_second_home_away));
}
// Search for solutions and fill day_tuples.
DecisionBuilder* const db = solver.MakePhase(signed_opponents,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
solver.NewSearch(db);
while (solver.NextSolution()) {
std::vector<int> solution;
for (int i = 0; i < num_teams; ++i) {
solution.push_back(signed_opponents[i]->Value());
}
day_tuples->Insert(solution);
}
solver.EndSearch();
LOG(INFO) << day_tuples->NumTuples()
<< " solutions to the one-day matching problem";
}
// Adds all constraints relating to one teams and the complete schedule.
void AddOneTeamConstraints(Solver* const solver,
const std::vector<IntVar*>& opponents,
const std::vector<IntVar*>& home_aways,
const std::vector<IntVar*>& signed_opponents,
const IntTupleSet& home_away_tuples,
IntVar* const break_var,
int num_teams) {
const int half_season = num_teams - 1;
// Each team meet all opponents once by half season.
for (int half = 0; half <= 1; ++half) {
for (int team_index = 0; team_index < num_teams; ++team_index) {
std::vector<IntVar*> tmp_vars;
for (int day = 0; day < half_season; ++day) {
tmp_vars.push_back(opponents[half * half_season + day]);
}
solver->AddConstraint(solver->MakeAllDifferent(tmp_vars));
}
}
// We meet each opponent once at home and once away per full season.
for (int team_index = 0; team_index < num_teams; ++team_index) {
solver->AddConstraint(solver->MakeAllDifferent(signed_opponents));
}
// Constraint per team on home_aways;
for (int i = 0; i < num_teams; ++i) {
std::vector<IntVar*> tmp_vars(home_aways);
tmp_vars.push_back(break_var);
solver->AddConstraint(
solver->MakeAllowedAssignments(tmp_vars, home_away_tuples));
}
}
// ----- Constraints for one team and all days -----
// Computes all valid tuples for home_away variables for a single team
// on the full lenght of the season.
void ComputeOneTeamHomeAwayTuples(int num_teams,
IntTupleSet* const home_away_tuples) {
LOG(INFO) << "Compute possible sequence of home and aways for any team.";
const int half_season = num_teams - 1;
const int full_season = 2 * half_season;
Solver solver("compute_home_aways");
std::vector<IntVar*> home_aways;
solver.MakeBoolVarArray(full_season, "home_away_", &home_aways);
for (int day = 0; day < full_season - 2; ++day) {
std::vector<IntVar*> tmp_vars;
tmp_vars.push_back(home_aways[day]);
tmp_vars.push_back(home_aways[day + 1]);
tmp_vars.push_back(home_aways[day + 2]);
IntVar* const partial_sum = solver.MakeSum(tmp_vars)->Var();
solver.AddConstraint(solver.MakeBetweenCt(partial_sum, 1, 2));
}
DecisionBuilder* const db = solver.MakePhase(home_aways,
Solver::CHOOSE_FIRST_UNBOUND,
Solver::ASSIGN_MIN_VALUE);
solver.NewSearch(db);
while (solver.NextSolution()) {
std::vector<int> solution;
for (int i = 0; i < full_season; ++i) {
solution.push_back(home_aways[i]->Value());
}
int breaks = 0;
for (int i = 0; i < full_season - 1; ++i) {
breaks += (solution[i] == solution[i + 1]);
}
solution.push_back(breaks);
home_away_tuples->Insert(solution);
}
solver.EndSearch();
LOG(INFO) << home_away_tuples->NumTuples()
<< " combination of home_aways for a team on the full season";
}
// ---------- Main solving method ----------
// Solves the sports scheduling problem with a given number of teams.
void SportsScheduling(int num_teams) {
const int half_season = num_teams - 1;
const int full_season = 2 * half_season;
Solver solver("Sports Scheduling");
// ----- Variables -----
// The index of the opponent of a team on a given day.
std::vector<std::vector<IntVar*> > opponents(num_teams);
// The location of the match (home or away).
std::vector<std::vector<IntVar*> > home_aways(num_teams);
// Disambiguated version of the opponent variable incorporating the
// home_away result.
std::vector<std::vector<IntVar*> > signed_opponents(num_teams);
for (int team_index = 0; team_index < num_teams; ++team_index) {
solver.MakeIntVarArray(full_season,
0,
num_teams - 1,
StringPrintf("opponent_%d_", team_index),
&opponents[team_index]);
solver.MakeBoolVarArray(full_season,
StringPrintf("home_away_%d_", team_index),
&home_aways[team_index]);
solver.MakeIntVarArray(full_season,
0,
2 * num_teams - 1,
StringPrintf("signed_opponent_%d", team_index),
&signed_opponents[team_index]);
}
// ----- Constraints -----
// Constraints on a given day.
IntTupleSet one_day_tuples(num_teams);
ComputeOneDayTuples(num_teams, &one_day_tuples);
for (int day = 0; day < full_season; ++day) {
std::vector<IntVar*> all_vars;
for (int i = 0; i < num_teams; ++i) {
all_vars.push_back(signed_opponents[i][day]);
}
solver.AddConstraint(
solver.MakeAllowedAssignments(all_vars, one_day_tuples));
}
// Links signed_opponents, home_away and opponents.
IntTupleSet one_day_one_team_tuples(3);
ComputeOneDayOneTeamTuples(num_teams, &one_day_one_team_tuples);
for (int day = 0; day < full_season; ++day) {
for (int team_index = 0; team_index < num_teams; ++team_index) {
AddOneDayOneTeamConstraint(&solver,
opponents[team_index][day],
home_aways[team_index][day],
signed_opponents[team_index][day],
one_day_one_team_tuples);
}
}
// Constraints on a team.
IntTupleSet home_away_tuples(full_season + 1);
ComputeOneTeamHomeAwayTuples(num_teams, &home_away_tuples);
std::vector<IntVar*> team_breaks;
solver.MakeIntVarArray(num_teams,
0,
full_season,
"team_break_",
&team_breaks);
for (int team = 0; team < num_teams; ++team) {
AddOneTeamConstraints(&solver,
opponents[team],
home_aways[team],
signed_opponents[team],
home_away_tuples,
team_breaks[team],
num_teams);
}
// ----- Search -----
std::vector<SearchMonitor*> monitors;
// Objective.
IntVar* const objective_var =
solver.MakeSum(team_breaks)->VarWithName("SumOfBreaks");
OptimizeVar* const objective_monitor = solver.MakeMinimize(objective_var, 1);
monitors.push_back(objective_monitor);
// Store all variables in a single array.
std::vector<IntVar*> all_signed_opponents;
for (int team_index = 0; team_index < num_teams; ++team_index) {
for (int day = 0; day < full_season; ++day) {
all_signed_opponents.push_back(signed_opponents[team_index][day]);
}
}
// Build default phase decision builder.
DefaultPhaseParameters parameters;
parameters.run_all_heuristics = FLAGS_run_all_heuristics;
parameters.heuristic_period = FLAGS_heuristics_period;
parameters.restart_log_size = FLAGS_restart_log_size;
DecisionBuilder* const db =
solver.MakeDefaultPhase(all_signed_opponents, parameters);
// Search log.
SearchMonitor* const log = solver.MakeSearchLog(1000000, objective_monitor);
monitors.push_back(log);
// Search limit.
SearchLimit* const limit = solver.MakeTimeLimit(FLAGS_time_limit);
monitors.push_back(limit);
// Solution collector.
SolutionCollector* const collector = solver.MakeLastSolutionCollector();
for (int team_index = 0; team_index < num_teams; ++team_index) {
collector->Add(opponents[team_index]);
collector->Add(home_aways[team_index]);
}
monitors.push_back(collector);
// And search.
solver.Solve(db, monitors);
// Display solution.
if (collector->solution_count() == 1) {
LOG(INFO) << "Solution found in " << solver.wall_time() << " ms, and "
<< solver.failures() << " failures.";
for (int team_index = 0; team_index < num_teams; ++team_index) {
string line;
for (int day = 0; day < full_season; ++day) {
const int opponent = collector->Value(0, opponents[team_index][day]);
const int home_away = collector->Value(0, home_aways[team_index][day]);
line += StringPrintf("%2d%s ", opponent, home_away ? "@" : " ");
}
LOG(INFO) << line;
}
}
}
} // namespace operations_research
static const char kUsage[] =
"Usage: see flags.\nThis program runs a sports scheduling problem."
"There is no output besides the debug LOGs of the solver.";
int main(int argc, char **argv) {
google::SetUsageMessage(kUsage);
google::ParseCommandLineFlags(&argc, &argv, true);
CHECK_EQ(0, FLAGS_num_teams % 2) << "The number of teams must be even";
CHECK_GE(FLAGS_num_teams, 2) << "At least 2 teams";
CHECK_LT(FLAGS_num_teams, 16) << "The model does not scale beyond 14 teams";
operations_research::SportsScheduling(FLAGS_num_teams);
return 0;
}

658
examples/cpp/strawberry_fields_with_column_generation.cc Normal file
View File

@@ -0,0 +1,658 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// Demonstration of column generation using LP toolkit.
//
// Column generation is the technique of generating columns (aka
// resource bundles aka variables) of the constraint matrix
// incrementally guided by feedback from the constraint duals
// (cost-of-resources). Frequently this lets one solve large problems
// efficiently, e.g. problems where the number of potential columns is
// exponentially large.
//
// Solves a covering problem taken from ITA Software recruiting web
// site:
//
// "Strawberries are growing in the cells of a rectangular field
// (grid). You want to build greenhouses to enclose the
// strawberries. Greenhouses are rectangular, axis-aligned with the
// field (i.e., not diagonal), and may not overlap. The cost of each
// greenhouse is $10 plus $1 per unit of area covered."
//
// Variables:
//
// for each Box (greenhouse), continuous variable b{x1,y1,x2,y2} in [0,1]
//
// Constraints:
//
// box limit:
// sum b{x1,y1,x2,y2) <= MAX_BOXES
// non-overlap (for each cell x,y):
// sum b{x1,y1,x2,y2} <= 1 (summed over containing x1<=x<=x2, y1<=y<=y2)
// coverage (for each cell x,y with a strawberry):
// sum b{x1,y1,x2,y2} = 1 (summed over containing x1<=x<=x2, y1<=y<=y2)
//
// Since the number of possible boxes is O(d^4) where d is the linear
// dimension, starts from singleton column (box) covering entire grid,
// ensuring solvability. Then iteratively the problem is solved and
// the constraint duals (aka reduced costs) used to guide the
// generation of a single new column (box), until convergence or a
// maximum number of iterations.
//
// No attempt is made to force integrality.
#include <stdio.h>
#include <string.h> // strlen
#include <map>
#include <string>
#include <utility>
#include <vector>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/macros.h"
#include "base/stringprintf.h"
#include "linear_solver/linear_solver.h"
using std::string;
DEFINE_bool(colgen_verbose, false, "print verbosely");
DEFINE_bool(colgen_complete, false, "generate all columns initially");
DEFINE_int32(colgen_max_iterations, 500, "max iterations");
DEFINE_string(colgen_solver, "clp", "solver - glpk or clp (default)");
DEFINE_int32(colgen_instance, -1, "Which instance to solve (0 - 9)");
namespace operations_research {
// ---------- Data Instances ----------
struct Instance {
int max_boxes;
int width;
int height;
const char* grid;
};
Instance kInstances[] = {
{ 4, 22, 6,
"..@@@@@..............."
"..@@@@@@........@@@..."
".....@@@@@......@@@..."
".......@@@@@@@@@@@@..."
".........@@@@@........"
".........@@@@@........"
},
{ 3, 13, 10,
"............."
"............."
"............."
"...@@@@......"
"...@@@@......"
"...@@@@......"
".......@@@..."
".......@@@..."
".......@@@..."
"............."
},
{ 4, 13, 9,
"............."
"..@.@.@......"
"...@.@.@....."
"..@.@.@......"
"..@.@.@......"
"...@.@.@....."
"....@.@......"
"..........@@@"
"..........@@@"
},
{ 4, 13, 9,
".........@..."
".........@..."
"@@@@@@@@@@..."
"..@......@..."
"..@......@..."
"..@......@..."
"..@@@@@@@@@@@"
"..@.........."
"..@.........."
},
{ 7, 25, 14,
"........................."
"..@@@@@@@@@@@@@@@@@@@@..."
"..@@@@@@@@@@@@@@@@@@@@..."
"..@@.................@..."
"..@@.................@..."
"..@@.......@@@.......@.@."
"..@@.......@@@.......@..."
"..@@...@@@@@@@@@@@@@@@..."
"..@@...@@@@@@@@@@@@@@@..."
"..@@.......@@@.......@..."
"..@@.......@@@.......@..."
"..@@.................@..."
"..@@.................@..."
"........................."
},
{ 6, 25, 16,
"........................."
"......@@@@@@@@@@@@@......"
"........................."
".....@..........@........"
".....@..........@........"
".....@......@............"
".....@......@.@@@@@@@...."
".....@......@............"
".....@......@.@@@@@@@...."
".....@......@............"
"....@@@@....@............"
"....@@@@....@............"
"..@@@@@@....@............"
"..@@@.......@............"
"..@@@...................."
"..@@@@@@@@@@@@@@@@@@@@@@@"
},
{ 5, 40, 18,
"........................................"
"........................................"
"...@@@@@@..............................."
"...@@@@@@..............................."
"...@@@@@@..............................."
"...@@@@@@.........@@@@@@@@@@............"
"...@@@@@@.........@@@@@@@@@@............"
"..................@@@@@@@@@@............"
"..................@@@@@@@@@@............"
".............@@@@@@@@@@@@@@@............"
".............@@@@@@@@@@@@@@@............"
"........@@@@@@@@@@@@...................."
"........@@@@@@@@@@@@...................."
"........@@@@@@.........................."
"........@@@@@@.........................."
"........................................"
"........................................"
"........................................"
},
{ 8, 40, 18,
"........................................"
"..@@.@.@.@.............................."
"..@@.@.@.@...............@.............."
"..@@.@.@.@............@................."
"..@@.@.@.@.............................."
"..@@.@.@.@.................@............"
"..@@.@..................@..............."
"..@@.@.................................."
"..@@.@.................................."
"..@@.@................@@@@.............."
"..@@.@..............@@@@@@@@............"
"..@@.@.................................."
"..@@.@..............@@@@@@@@............"
"..@@.@.................................."
"..@@.@................@@@@.............."
"..@@.@.................................."
"..@@.@.................................."
"........................................"
},
{ 10, 40, 19,
"@@@@@..................................."
"@@@@@..................................."
"@@@@@..................................."
"@@@@@..................................."
"@@@@@..................................."
"@@@@@...........@@@@@@@@@@@............."
"@@@@@...........@@@@@@@@@@@............."
"....................@@@@................"
"....................@@@@................"
"....................@@@@................"
"....................@@@@................"
"....................@@@@................"
"...............@@@@@@@@@@@@@@..........."
"...............@@@@@@@@@@@@@@..........."
".......@@@@@@@@@@@@@@@@@@@@@@..........."
".......@@@@@@@@@........................"
"........................................"
"........................................"
"........................................"
},
{ 10, 40, 25,
"...................@...................."
"...............@@@@@@@@@................"
"............@@@.........@@@............."
"...........@...............@............"
"..........@.................@..........."
".........@...................@.........."
".........@...................@.........."
".........@.....@@......@@....@.........."
"........@.....@@@@....@@@@....@........."
"........@.....................@........."
"........@.....................@........."
"........@..........@@.........@........."
".......@@..........@@.........@@........"
"........@.....................@........."
"........@.....................@........."
"........@......@@@@@@@@@......@........."
"........@......@@@@@@@@@......@........."
".........@...................@.........."
".........@...................@.........."
".........@...................@.........."
"..........@.................@..........."
"...........@...............@............"
"............@@@.........@@@............."
"...............@@@@@@@@@................"
"...................@...................."
}
};
const int kInstanceCount = 10;
// ---------- Box ---------
class Box {
public:
static const int kAreaCost = 1;
static const int kFixedCost = 10;
Box() {}
Box(int x_min, int x_max, int y_min, int y_max)
: x_min_(x_min), x_max_(x_max), y_min_(y_min), y_max_(y_max) {
CHECK_GE(x_max, x_min);
CHECK_GE(y_max, y_min);
}
int x_min() const { return x_min_; }
int x_max() const { return x_max_; }
int y_min() const { return y_min_; }
int y_max() const { return y_max_; }
// Lexicographic order
int Compare(const Box& box) const {
int c;
if ((c = (x_min() - box.x_min())) != 0) return c;
if ((c = (x_max() - box.x_max())) != 0) return c;
if ((c = (y_min() - box.y_min())) != 0) return c;
return y_max() - box.y_max();
}
bool Contains(int x, int y) const {
return x >= x_min() && x <= x_max() && y >= y_min() && y <= y_max();
}
int Cost() const {
return kAreaCost * (x_max() - x_min() + 1) * (y_max() - y_min() + 1)
+ kFixedCost;
}
string DebugString() const {
return StringPrintf("[%d,%dx%d,%d]c%d",
x_min(), y_min(), x_max(), y_max(), Cost());
}
private:
int x_min_;
int x_max_;
int y_min_;
int y_max_;
};
struct BoxLessThan {
bool operator()(const Box& b1, const Box& b2) const {
return b1.Compare(b2) < 0;
}
};
// ---------- Covering Problem ---------
class CoveringProblem {
public:
// Grid is a row-major string of length width*height with '@' for an
// occupied cell (strawberry) and '.' for an empty cell. Solver is
// not owned.
CoveringProblem(MPSolver* const solver, const Instance& instance)
: solver_(solver),
max_boxes_(instance.max_boxes),
width_(instance.width),
height_(instance.height),
grid_(instance.grid) {}
// Constructs initial variables and constraints. Initial column
// (box) covers entire grid, ensuring feasibility.
bool Init() {
// Check consistency.
int size = strlen(grid_);
if (size != area()) {
return false;
}
for (int i = 0; i < size; ++i) {
char c = grid_[i];
if ((c != '@') && (c != '.')) return false;
}
AddCellConstraints(); // sum for every cell is <=1 or =1
AddMaxBoxesConstraint(); // sum of box variables is <= max_boxes()
if (!FLAGS_colgen_complete) {
AddBox(Box(0, width()-1, 0, height()-1)); // grid-covering box
} else {
// Naive alternative to column generation - generate all boxes;
// works fine for smaller problems, too slow for big.
for (int y_min = 0; y_min < height(); ++y_min) {
for (int y_max = y_min; y_max < height(); ++y_max) {
for (int x_min = 0; x_min < width(); ++x_min) {
for (int x_max = x_min; x_max < width(); ++x_max) {
AddBox(Box(x_min, x_max, y_min, y_max));
}
}
}
}
}
return true;
}
int width() const { return width_; }
int height() const { return height_; }
int area() const { return width() * height(); }
int max_boxes() const { return max_boxes_; }
bool IsCellOccupied(int x, int y) const { return grid_[index(x, y)] == '@'; }
// Calculates reduced costs for each possible Box and if any is
// negative (improves cost), returns reduced cost and set target to
// the most-negative (steepest descent) one - otherwise returns 0..
//
// For a problem in standard form 'minimize c*x s.t. Ax<=b, x>=0'
// the reduced cost vector is c - transp(y) * A where y is the dual
// cost column vector.
//
// For this covering problem, in which all coefficients in A are 0
// or 1, this reduces to:
//
// reduced_cost(box) =
//
// box.Cost() - sum_{enclosed cell} cell_constraint->dual_value()
// - max_boxes_constraint_->dual_value()
//
// Since there are O(d^4) boxes, we don't also want O(d^2) sum for
// each, so pre-calculate sums of cell duals for all rectangles with
// upper-left at 0, 0, and use these to calculate the sum in
// constant time using the standard inclusion-exclusion trick.
double GetOptimalBox(Box* const target) {
// Cost change threshold for new Box
const double kCostChangeThreshold = -.01;
// Precomputes the sum of reduced costs for every upper-left
// rectangle.
std::vector<double> upper_left_sums(area());
ComputeUpperLeftSums(&upper_left_sums);
const double max_boxes_dual = max_boxes_constraint_->dual_value();
double best_reduced_cost = kCostChangeThreshold;
Box best_box;
for (int y_min = 0; y_min < height(); ++y_min) {
for (int y_max = y_min; y_max < height(); ++y_max) {
for (int x_min = 0; x_min < width(); ++x_min) {
for (int x_max = x_min; x_max < width(); ++x_max) {
Box box(x_min, x_max, y_min, y_max);
const double cell_coverage_dual = // inclusion-exclusion
+ zero_access(upper_left_sums, x_max, y_max)
- zero_access(upper_left_sums, x_max, y_min - 1)
- zero_access(upper_left_sums, x_min - 1, y_max)
+ zero_access(upper_left_sums, x_min - 1, y_min - 1);
// All coefficients for new column are 1, so no need to
// multiply constraint duals by any coefficients when
// computing the reduced cost.
const double reduced_cost = box.Cost()
- (cell_coverage_dual + max_boxes_dual);
if (reduced_cost < best_reduced_cost) {
// Even with negative reduced cost, the box may already
// exist, and even be basic (part of solution)! This
// counterintuitive situation is due to the problem's
// many redundant linear equality constraints: many
// steepest-edge pivot moves will be of zero-length.
// Ideally one would want to check the length of the
// move but that is difficult without access to the
// internals of the solver (e.g., access to B^-1 in the
// simplex algorithm).
if (boxes_.find(box) == boxes_.end()) {
best_reduced_cost = reduced_cost;
best_box = box;
}
}
}
}
}
}
if (best_reduced_cost < kCostChangeThreshold) {
if (target) {
*target = best_box;
}
return best_reduced_cost;
}
return 0;
}
// Add continuous [0,1] box variable with box.Cost() as objective
// coefficient. Add to cell constraint of all enclosed cells.
MPVariable* AddBox(const Box& box) {
CHECK(boxes_.find(box) == boxes_.end());
MPVariable* const var = solver_->MakeNumVar(0., 1., box.DebugString());
solver_->SetObjectiveCoefficient(var, box.Cost());
max_boxes_constraint_->SetCoefficient(var, 1.0);
for (int y = box.y_min(); y <= box.y_max(); ++y) {
for (int x = box.x_min(); x <= box.x_max(); ++x) {
cell(x, y)->SetCoefficient(var, 1.0);
}
}
boxes_[box] = var;
return var;
}
string PrintGrid() const {
string output = StringPrintf("width = %d, height = %d, max_boxes = %d\n",
width_,
height_,
max_boxes_);
for (int y = 0; y < height_; ++y) {
StringAppendF(&output,
"%s\n",
string(grid_ + width_ * y, width_).c_str());
}
return output;
}
// Prints covering - total cost, those variables with non-zero value,
// and graphical depiction of covering using upper case letters for
// integral coverage and lower case for coverage using combination
// of fractional boxes.
string PrintCovering() const {
static const double kTolerance = 1e-5;
string output = StringPrintf("cost = %lf\n", solver_->objective_value());
scoped_array<char> display(new char[(width_ + 1) * height_ + 1]);
for (int y = 0; y < height_; ++y) {
memcpy(display.get() + y * (width_ + 1),
grid_ + width_ * y,
width_); // Copy the original line.
display[y * (width_ + 1) + width_] = '\n';
}
display[height_ * (width_ + 1)] = '\0';
int active_box_index = 0;
for (BoxTable::const_iterator i = boxes_.begin(); i != boxes_.end(); ++i) {
const double value = i->second->solution_value();
if (value > kTolerance) {
const char box_character =
(i->second->solution_value() >= (1. - kTolerance) ? 'A' : 'a') +
active_box_index++;
StringAppendF(&output, "%c: box %s with value %lf\n",
box_character,
i->first.DebugString().c_str(),
value);
const Box& box = i->first;
for (int x = box.x_min(); x <= box.x_max(); ++x) {
for (int y = box.y_min(); y <= box.y_max(); ++y) {
display[x + y * (width_ + 1)] = box_character;
}
}
}
}
StringAppendF(&output, "%s", display.get());
return output;
}
protected:
int index(int x, int y) const { return width_ * y + x; }
MPConstraint* cell(int x, int y) { return cells_[index(x, y)]; }
const MPConstraint* cell(int x, int y) const { return cells_[index(x, y)]; }
// Adds constraints that every cell is covered at most once, exactly
// once if occupied.
void AddCellConstraints() {
cells_.resize(area());
for (int y = 0; y < height(); ++y) {
for (int x = 0; x < width(); ++x) {
cells_[index(x, y)] = solver_->MakeRowConstraint
(IsCellOccupied(x, y) ? 1. : 0., 1.);
}
}
}
// Adds constraint on maximum number of boxes used to cover.
void AddMaxBoxesConstraint() {
max_boxes_constraint_ = solver_->MakeRowConstraint(0., max_boxes());
}
// Gets 2d array element, returning 0 if out-of-bounds.
double zero_access(const std::vector<double>& array, int x, int y) const {
if (x < 0 || y < 0) {
return 0;
}
return array[index(x, y)];
}
// Precomputes the sum of reduced costs for every upper-left
// rectangle.
void ComputeUpperLeftSums(std::vector<double>* upper_left_sums) const {
for (int y = 0; y < height(); ++y) {
for (int x = 0; x < width(); ++x) {
upper_left_sums->operator[](index(x, y)) =
cell(x, y)->dual_value()
+ zero_access(*upper_left_sums, x - 1, y)
+ zero_access(*upper_left_sums, x, y - 1)
- zero_access(*upper_left_sums, x - 1, y - 1);
}
}
}
typedef std::map<Box, MPVariable*, BoxLessThan> BoxTable;
MPSolver* const solver_; // not owned
const int max_boxes_;
const int width_;
const int height_;
const char* const grid_;
std::vector<MPConstraint*> cells_;
BoxTable boxes_;
MPConstraint* max_boxes_constraint_;
};
// ---------- Main Solve Method ----------
// Solves iteratively using delayed column generation, up to maximum
// number of steps.
void SolveInstance(const Instance& instance,
MPSolver::OptimizationProblemType solver_type) {
// Prepares the solver.
MPSolver solver("ColumnGeneration", solver_type);
solver.SuppressOutput();
solver.SetMinimization();
// Construct problem.
CoveringProblem problem(&solver, instance);
CHECK(problem.Init());
LOG(INFO) << "Initial problem:\n" << problem.PrintGrid();
int step_number = 0;
while (step_number < FLAGS_colgen_max_iterations) {
if (FLAGS_colgen_verbose) {
LOG(INFO) << "Step number " << step_number;
}
// Solve with existing columns.
CHECK_EQ(MPSolver::OPTIMAL, solver.Solve());
if (FLAGS_colgen_verbose) {
LOG(INFO) << problem.PrintCovering();
}
// Find optimal new column to add, or stop if none.
Box box;
const double reduced_cost = problem.GetOptimalBox(&box);
if (reduced_cost >= 0) {
break;
}
// Add new column to problem.
if (FLAGS_colgen_verbose) {
LOG(INFO) << "Adding " << box.DebugString()
<< ", reduced_cost =" << reduced_cost;
}
problem.AddBox(box);
++step_number;
}
if (step_number >= FLAGS_colgen_max_iterations) {
// Solve one last time with all generated columns.
CHECK_EQ(MPSolver::OPTIMAL, solver.Solve());
}
LOG(INFO) << step_number << " columns added";
LOG(INFO) << "Final coverage: " << problem.PrintCovering();
}
} // namespace operations_research
int main(int argc, char** argv) {
string usage = "column_generation\n";
usage += " --colgen_verbose print verbosely\n";
usage += " --colgen_max_iterations <n> max columns to generate\n";
usage += " --colgen_complete generate all columns at start\n";
operations_research::MPSolver::OptimizationProblemType solver_type;
bool found = false;
#if defined(USE_CLP)
if (FLAGS_colgen_solver == "clp") {
solver_type = operations_research::MPSolver::CLP_LINEAR_PROGRAMMING;
found = true;
}
#endif // USE_CLP
#if defined(USE_GLPK)
if (FLAGS_colgen_solver == "glpk") {
solver_type = operations_research::MPSolver::GLPK_LINEAR_PROGRAMMING;
found = true;
}
#endif // USE_GLPK
if (!found) {
LOG(ERROR) << "Unknown solver " << FLAGS_colgen_solver;
return 1;
}
if (FLAGS_colgen_instance == -1) {
for (int i = 0; i < operations_research::kInstanceCount; ++i) {
const operations_research::Instance& instance =
operations_research::kInstances[i];
operations_research::SolveInstance(instance, solver_type);
}
} else {
CHECK_GE(FLAGS_colgen_instance, 0);
CHECK_LT(FLAGS_colgen_instance, operations_research::kInstanceCount);
const operations_research::Instance& instance =
operations_research::kInstances[FLAGS_colgen_instance];
operations_research::SolveInstance(instance, solver_type);
}
return 0;
}

160
examples/cpp/tsp.cc Normal file
View File

@@ -0,0 +1,160 @@
// Copyright 2010-2011 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Traveling Salesman Sample.
//
// This is a sample using the routing library to solve a Traveling Salesman
// Problem.
// The description of the problem can be found here:
// http://en.wikipedia.org/wiki/Travelling_salesman_problem.
// For small problems one can use the hamiltonian path library directly (cf
// graph/hamiltonian_path.h).
// The optimization engine uses local search to improve solutions, first
// solutions being generated using a cheapest addition heuristic.
// Optionally one can randomly forbid a set of random connections between nodes
// (forbidden arcs).
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/scoped_ptr.h"
#include "base/join.h"
#include "constraint_solver/routing.h"
#include "base/random.h"
using operations_research::Assignment;
using operations_research::RoutingModel;
using operations_research::ACMRandom;
using operations_research::StrCat;
using operations_research::scoped_array;
DEFINE_int32(tsp_size, 10, "Size of Traveling Salesman Problem instance.");
DEFINE_bool(tsp_use_random_matrix, true, "Use random cost matrix.");
DEFINE_int32(tsp_random_forbidden_connections, 0,
"Number of random forbidden connections.");
DEFINE_bool(tsp_use_deterministic_random_seed, false,
"Use deterministic random seeds.");
// Random seed generator.
int32 GetSeed() {
if (FLAGS_tsp_use_deterministic_random_seed) {
return ACMRandom::DeterministicSeed();
} else {
return ACMRandom::HostnamePidTimeSeed();
}
}
// Cost/distance functions.
// Sample function.
int64 MyDistance(RoutingModel::NodeIndex from, RoutingModel::NodeIndex to) {
// Put your distance code here.
return (from + to).value(); // for instance
}
// Random matrix.
class RandomMatrix {
public:
explicit RandomMatrix(int size) : size_(size) {}
void Initialize() {
matrix_.reset(new int64[size_ * size_]);
const int64 kDistanceMax = 100;
ACMRandom randomizer(GetSeed());
for (RoutingModel::NodeIndex from = RoutingModel::kFirstNode; from < size_;
++from) {
for (RoutingModel::NodeIndex to = RoutingModel::kFirstNode; to < size_;
++to) {
if (to != from) {
matrix_[MatrixIndex(from, to)] = randomizer.Uniform(kDistanceMax);
} else {
matrix_[MatrixIndex(from, to)] = 0LL;
}
}
}
}
int64 Distance(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return matrix_[MatrixIndex(from, to)];
}
private:
int64 MatrixIndex(RoutingModel::NodeIndex from,
RoutingModel::NodeIndex to) const {
return (from * size_ + to).value();
}
scoped_array<int64> matrix_;
const int size_;
};
int main(int argc, char **argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_tsp_size > 0) {
// TSP of size FLAGS_tsp_size.
// Second argument = 1 to build a single tour (it's a TSP).
// Nodes are indexed from 0 to FLAGS_tsp_size - 1, by default the start of
// the route is node 0.
RoutingModel routing(FLAGS_tsp_size, 1);
// Setting first solution heuristic (cheapest addition).
routing.SetCommandLineOption("routing_first_solution", "PathCheapestArc");
// Disabling Large Neighborhood Search, comment out to activate it.
routing.SetCommandLineOption("routing_no_lns", "true");
// Setting the cost function.
// Put a permanent callback to the distance accessor here. The callback
// has the following signature: ResultCallback2<int64, int64, int64>.
// The two arguments are the from and to node inidices.
RandomMatrix matrix(FLAGS_tsp_size);
if (FLAGS_tsp_use_random_matrix) {
matrix.Initialize();
routing.SetCost(NewPermanentCallback(&matrix, &RandomMatrix::Distance));
} else {
routing.SetCost(NewPermanentCallback(MyDistance));
}
// Forbid node connections (randomly).
ACMRandom randomizer(GetSeed());
int64 forbidden_connections = 0;
while (forbidden_connections < FLAGS_tsp_random_forbidden_connections) {
const int64 from = randomizer.Uniform(FLAGS_tsp_size - 1);
const int64 to = randomizer.Uniform(FLAGS_tsp_size - 1) + 1;
if (routing.NextVar(from)->Contains(to)) {
LG << "Forbidding connection " << from << " -> " << to;
routing.NextVar(from)->RemoveValue(to);
++forbidden_connections;
}
}
// Solve, returns a solution if any (owned by RoutingModel).
const Assignment* solution = routing.Solve();
if (solution != NULL) {
// Solution cost.
LG << "Cost " << solution->ObjectiveValue();
// Inspect solution.
// Only one route here; otherwise iterate from 0 to routing.vehicles() - 1
const int route_number = 0;
string route;
for (int64 node = routing.Start(route_number);
!routing.IsEnd(node);
node = solution->Value(routing.NextVar(node))) {
route = StrCat(route, StrCat(node, " -> "));
}
route = StrCat(route, "0");
LG << route;
} else {
LG << "No solution found.";
}
} else {
LG << "Specify an instance size greater than 0.";
}
return 0;
}