From 37c33757636b45d5949bd6c655495586eec9f3a3 Mon Sep 17 00:00:00 2001 From: "lperron@google.com" Date: Thu, 5 Apr 2012 20:24:29 +0000 Subject: [PATCH] Untested version of the sortedness constraint --- src/constraint_solver/alldiff_cst.cc | 508 +++++++++++++-------- src/constraint_solver/constraint_solver.cc | 1 + src/constraint_solver/constraint_solver.h | 6 + src/constraint_solver/io.cc | 13 + 4 files changed, 347 insertions(+), 181 deletions(-) diff --git a/src/constraint_solver/alldiff_cst.cc b/src/constraint_solver/alldiff_cst.cc index 48703e6a7c..65256a526c 100644 --- a/src/constraint_solver/alldiff_cst.cc +++ b/src/constraint_solver/alldiff_cst.cc @@ -23,6 +23,7 @@ #include "base/scoped_ptr.h" #include "constraint_solver/constraint_solver.h" #include "constraint_solver/constraint_solveri.h" +#include "util/const_ptr_array.h" #include "util/string_array.h" namespace operations_research { @@ -151,60 +152,18 @@ bool ValueAllDifferent::AllMoves() { // ---------- Bounds All Different ---------- // See http://www.cs.uwaterloo.ca/~cquimper/Papers/ijcai03_TR.pdf for details. -struct Interval { - int64 min; - int64 max; - int min_rank; - int max_rank; -}; - -// TODO(user) : use better sort, use bounding boxes of modifications to -// improve the sorting (only modified vars). - -// This method is used by the STL sort. -struct CompareIntervalMin { - bool operator()(const Interval* i1, const Interval* i2) { - return (i1->min < i2->min); - } -}; - -// This method is used by the STL sort. -struct CompareIntervalMax { - bool operator()(const Interval* i1, const Interval* i2) { - return (i1->max < i2->max); - } -}; - -void PathSet(int start, int end, int to, int* const tree) { - int k = start; - int l = start; - while (l != end) { - k = l; - l = tree[k]; - tree[k] = to; - } -} - -int PathMin(const int* const tree, int index) { - int i = index; - while (tree[i] < i) { - i = tree[i]; - } - return i; -} - -int PathMax(const int* const tree, int index) { - int i = index; - while (tree[i] > i) { - i = tree[i]; - } - return i; -} - -class BoundsAllDifferent : public BaseAllDifferent { +class RangeBipartiteMatching { public: - BoundsAllDifferent(Solver* const s, const IntVar* const * vars, int size) - : BaseAllDifferent(s, vars, size), + struct Interval { + int64 min; + int64 max; + int min_rank; + int max_rank; + }; + + RangeBipartiteMatching(Solver* const solver, int size) + : solver_(solver), + size_(size), intervals_(new Interval[size + 1]), min_sorted_(new Interval*[size]), max_sorted_(new Interval*[size]), @@ -218,7 +177,214 @@ class BoundsAllDifferent : public BaseAllDifferent { min_sorted_[i] = max_sorted_[i]; } } + + void SetRange(int index, int64 imin, int64 imax) { + intervals_[index].min = imin; + intervals_[index].max = imax; + } + + bool Propagate() { + SortArray(); + + bool modified = PropagateMin(); + modified |= PropagateMax(); + return modified; + } + + int64 Min(int index) const { + return intervals_[index].min; + } + + int64 Max(int index) const { + return intervals_[index].max; + } + + private: + // This method with sort the min_sorted_ and max_sorted_ arrays and fill + // the bounds_ array (and set the active_size_ counter). + void SortArray() { + std::sort(min_sorted_.get(), + min_sorted_.get() + size_, + CompareIntervalMin()); + std::sort(max_sorted_.get(), + max_sorted_.get() + size_, + CompareIntervalMax()); + + int64 min = min_sorted_[0]->min; + int64 max = max_sorted_[0]->max + 1; + int64 last = min - 2; + bounds_[0] = last; + + int i = 0; + int j = 0; + int nb = 0; + for (;;) { // merge min_sorted_[] and max_sorted_[] into bounds_[]. + if (i < size_ && min <= max) { // make sure min_sorted_ exhausted first. + if (min != last) { + last = min; + bounds_[++nb] = last; + } + min_sorted_[i]->min_rank = nb; + if (++i < size_) { + min = min_sorted_[i]->min; + } + } else { + if (max != last) { + last = max; + bounds_[++nb] = last; + } + max_sorted_[j]->max_rank = nb; + if (++j == size_) { + break; + } + max = max_sorted_[j]->max + 1; + } + } + active_size_ = nb; + bounds_[nb + 1] = bounds_[nb] + 2; + } + + // These two methods will actually do the new bounds computation. + bool PropagateMin() { + bool modified = false; + + for (int i = 1; i <= active_size_ + 1; ++i) { + hall_[i] = i - 1; + tree_[i] = i - 1; + diff_[i] = bounds_[i] - bounds_[i - 1]; + } + // visit intervals in increasing max order + for (int i = 0; i < size_; ++i) { + const int x = max_sorted_[i]->min_rank; + const int y = max_sorted_[i]->max_rank; + int z = PathMax(tree_.get(), x + 1); + int j = tree_[z]; + if (--diff_[z] == 0) { + tree_[z] = z + 1; + z = PathMax(tree_.get(), z + 1); + tree_[z] = j; + } + PathSet(x + 1, z, z, tree_.get()); // path compression + if (diff_[z] < bounds_[z] - bounds_[y]) { + solver_->Fail(); + } + if (hall_[x] > x) { + int w = PathMax(hall_.get(), hall_[x]); + max_sorted_[i]->min = bounds_[w]; + PathSet(x, w, w, hall_.get()); // path compression + modified = true; + } + if (diff_[z] == bounds_[z] - bounds_[y]) { + PathSet(hall_[y], j - 1, y, hall_.get()); // mark hall interval + hall_[y] = j - 1; + } + } + return modified; + } + + + bool PropagateMax() { + bool modified = false; + + for (int i = 0; i<= active_size_; i++) { + tree_[i] = i + 1; + hall_[i] = i + 1; + diff_[i] = bounds_[i + 1] - bounds_[i]; + } + // visit intervals in decreasing min order + for (int i = size_ - 1; i >= 0; --i) { + const int x = min_sorted_[i]->max_rank; + const int y = min_sorted_[i]->min_rank; + int z = PathMin(tree_.get(), x - 1); + int j = tree_[z]; + if (--diff_[z] == 0) { + tree_[z] = z - 1; + z = PathMin(tree_.get(), z - 1); + tree_[z] = j; + } + PathSet(x - 1, z, z, tree_.get()); + if (diff_[z] < bounds_[y] - bounds_[z]) { + solver_->Fail(); + // useless. Should have been caught by the PropagateMin() method. + } + if (hall_[x] < x) { + int w = PathMin(hall_.get(), hall_[x]); + min_sorted_[i]->max = bounds_[w] - 1; + PathSet(x, w, w, hall_.get()); + modified = true; + } + if (diff_[z] == bounds_[y] - bounds_[z]) { + PathSet(hall_[y], j + 1, y, hall_.get()); + hall_[y] = j + 1; + } + } + return modified; + } + + // TODO(user) : use better sort, use bounding boxes of modifications to + // improve the sorting (only modified vars). + + // This method is used by the STL sort. + struct CompareIntervalMin { + bool operator()(const Interval* i1, const Interval* i2) { + return (i1->min < i2->min); + } + }; + + // This method is used by the STL sort. + struct CompareIntervalMax { + bool operator()(const Interval* i1, const Interval* i2) { + return (i1->max < i2->max); + } + }; + + void PathSet(int start, int end, int to, int* const tree) { + int k = start; + int l = start; + while (l != end) { + k = l; + l = tree[k]; + tree[k] = to; + } + } + + int PathMin(const int* const tree, int index) { + int i = index; + while (tree[i] < i) { + i = tree[i]; + } + return i; + } + + int PathMax(const int* const tree, int index) { + int i = index; + while (tree[i] > i) { + i = tree[i]; + } + return i; + } + + Solver* const solver_; + const int size_; + scoped_array intervals_; + scoped_array min_sorted_; + scoped_array max_sorted_; + // bounds_[1..active_size_] hold set of min & max in the n intervals_ + // while bounds_[0] and bounds_[active_size_ + 1] allow sentinels. + scoped_array bounds_; + scoped_array tree_; // tree links. + scoped_array diff_; // diffs between critical capacities. + scoped_array hall_; // hall interval links. + int active_size_; +}; + +class BoundsAllDifferent : public BaseAllDifferent { + public: + BoundsAllDifferent(Solver* const s, const IntVar* const * vars, int size) + : BaseAllDifferent(s, vars, size), matching_(s, size) {} + virtual ~BoundsAllDifferent() {} + virtual void Post() { Demon* range = MakeDelayedConstraintDemon0( solver(), @@ -248,18 +414,13 @@ class BoundsAllDifferent : public BaseAllDifferent { virtual void IncrementalPropagate() { for (int i = 0; i < size_; ++i) { - intervals_[i].min = vars_[i]->Min(); - intervals_[i].max = vars_[i]->Max(); + matching_.SetRange(i, vars_[i]->Min(), vars_[i]->Max()); } - SortArray(); - bool modified = PropagateMin(); - modified |= PropagateMax(); - - if (modified) { + if (matching_.Propagate()) { for (int i = 0; i < size_; ++i) { - vars_[i]->SetRange(intervals_[i].min, intervals_[i].max); + vars_[i]->SetRange(matching_.Min(i), matching_.Max(i)); } } } @@ -288,139 +449,117 @@ class BoundsAllDifferent : public BaseAllDifferent { } private: - // This method with sort the min_sorted_ and max_sorted_ arrays and fill - // the bounds_ array (and set the active_size_ counter). - void SortArray(); - - // These two methods will actually do the new bounds computation. - bool PropagateMin(); - bool PropagateMax(); - - int64 stamp_; - scoped_array intervals_; - scoped_array min_sorted_; - scoped_array max_sorted_; - // bounds_[1..active_size_] hold set of min & max in the n intervals_ - // while bounds_[0] and bounds_[active_size_ + 1] allow sentinels. - scoped_array bounds_; - scoped_array tree_; // tree links. - scoped_array diff_; // diffs between critical capacities. - scoped_array hall_; // hall interval links. - int active_size_; + RangeBipartiteMatching matching_; }; -void BoundsAllDifferent::SortArray() { - std::sort(min_sorted_.get(), min_sorted_.get() + size_, CompareIntervalMin()); - std::sort(max_sorted_.get(), max_sorted_.get() + size_, CompareIntervalMax()); - - int64 min = min_sorted_[0]->min; - int64 max = max_sorted_[0]->max + 1; - int64 last = min - 2; - bounds_[0] = last; - - int i = 0; - int j = 0; - int nb = 0; - for (;;) { // merge min_sorted_[] and max_sorted_[] into bounds_[]. - if (i < size_ && min <= max) { // make sure min_sorted_ exhausted first. - if (min != last) { - last = min; - bounds_[++nb] = last; - } - min_sorted_[i]->min_rank = nb; - if (++i < size_) { - min = min_sorted_[i]->min; - } - } else { - if (max != last) { - last = max; - bounds_[++nb] = last; - } - max_sorted_[j]->max_rank = nb; - if (++j == size_) { - break; - } - max = max_sorted_[j]->max + 1; +class SortConstraint : public Constraint { + public: + SortConstraint(Solver* const solver, + const vector& original_vars, + const vector& sorted_vars) + : Constraint(solver), + ovars_(original_vars), + svars_(sorted_vars), + mins_(NULL), + maxs_(NULL), + size_(original_vars.size()), + matching_(solver, size_) { + if (size_ > 0) { + mins_.reset(new int64[size_]); + memset(mins_.get(), 0, size_ * sizeof(*mins_.get())); + maxs_.reset(new int64[size_]); + memset(maxs_.get(), 0, size_ * sizeof(*maxs_.get())); } } - active_size_ = nb; - bounds_[nb + 1] = bounds_[nb] + 2; -} -bool BoundsAllDifferent::PropagateMin() { - bool modified = false; + virtual ~SortConstraint() {} - for (int i = 1; i <= active_size_ + 1; ++i) { - hall_[i] = i - 1; - tree_[i] = i - 1; - diff_[i] = bounds_[i] - bounds_[i - 1]; - } - for (int i = 0; i < size_; ++i) { // visit intervals in increasing max order - const int x = max_sorted_[i]->min_rank; - const int y = max_sorted_[i]->max_rank; - int z = PathMax(tree_.get(), x + 1); - int j = tree_[z]; - if (--diff_[z] == 0) { - tree_[z] = z + 1; - z = PathMax(tree_.get(), z + 1); - tree_[z] = j; + virtual void Post() { + Demon* const demon = + solver()->MakeDelayedConstraintInitialPropagateCallback(this); + for (int i = 0; i < size_; ++i) { + ovars_[i]->WhenRange(demon); + svars_[i]->WhenRange(demon); } - PathSet(x + 1, z, z, tree_.get()); // path compression - if (diff_[z] < bounds_[z] - bounds_[y]) { + } + + virtual void InitialPropagate() { + for (int i = 0; i < size_; ++i) { + int64 vmin = 0; + int64 vmax = 0; + ovars_[i]->Range(&vmin, &vmax); + mins_[i] = vmin; + maxs_[i] = vmax; + } + // One direction, from variables to sorted variables. + std::sort(mins_.get(), mins_.get() + size_); + std::sort(maxs_.get(), maxs_.get() + size_); + for (int i = 0; i < size_; ++i) { + svars_[i]->SetRange(mins_[i], maxs_[i]); + } + // Reverse direction. + for (int i = 0; i < size_; ++i) { + int64 imin = 0; + int64 imax = 0; + FindIntersectionRange(i, &imin, &imax); + matching_.SetRange(i, imin, imax); + } + if (matching_.Propagate()) { + for (int i = 0; i < size_; ++i) { + const int64 vmin = mins_[matching_.Min(i)]; + const int64 vmax = maxs_[matching_.Max(i)]; + ovars_[i]->SetRange(vmin, vmax); + } + } + } + + virtual void Accept(ModelVisitor* const visitor) const { + visitor->BeginVisitConstraint(ModelVisitor::kSorted, this); + visitor->VisitIntegerVariableArrayArgument(ModelVisitor::kVarsArgument, + ovars_); + visitor->VisitIntegerVariableArrayArgument(ModelVisitor::kTargetArgument, + svars_); + visitor->EndVisitConstraint(ModelVisitor::kSorted, this); + } + + virtual string DebugString() const { + return StringPrintf("Sort(%s, %s)", + ovars_.DebugString().c_str(), + svars_.DebugString().c_str()); + } + + private: + void FindIntersectionRange(int index, + int64* const range_min, + int64* const range_max) const { + // Naive version. + int64 imin = 0; + while (imin < size_ && NotIntersect(index, imin)) { + imin++; + } + if (imin == size_) { solver()->Fail(); } - if (hall_[x] > x) { - int w = PathMax(hall_.get(), hall_[x]); - max_sorted_[i]->min = bounds_[w]; - PathSet(x, w, w, hall_.get()); // path compression - modified = true; - } - if (diff_[z] == bounds_[z] - bounds_[y]) { - PathSet(hall_[y], j - 1, y, hall_.get()); // mark hall interval - hall_[y] = j - 1; + int64 imax = size_ - 1; + while (imax >= 0 && NotIntersect(index, imax)) { + imax--; } + *range_min = imin; + *range_max = imax; } - return modified; -} - -bool BoundsAllDifferent::PropagateMax() { - bool modified = false; - - for (int i = 0; i<= active_size_; i++) { - tree_[i] = i + 1; - hall_[i] = i + 1; - diff_[i] = bounds_[i + 1] - bounds_[i]; + bool NotIntersect(int oindex, int sindex) const { + return ovars_[oindex]->Min() > svars_[sindex]->Max() || + ovars_[oindex]->Max() < svars_[sindex]->Min(); } - // visit intervals in decreasing min order - for (int i = size_ - 1; i >= 0; --i) { - const int x = min_sorted_[i]->max_rank; - const int y = min_sorted_[i]->min_rank; - int z = PathMin(tree_.get(), x - 1); - int j = tree_[z]; - if (--diff_[z] == 0) { - tree_[z] = z - 1; - z = PathMin(tree_.get(), z - 1); - tree_[z] = j; - } - PathSet(x - 1, z, z, tree_.get()); - if (diff_[z] < bounds_[y] - bounds_[z]) { - solver()->Fail(); - // useless. Should have been caught by the PropagateMin() method. - } - if (hall_[x] < x) { - int w = PathMin(hall_.get(), hall_[x]); - min_sorted_[i]->max = bounds_[w] - 1; - PathSet(x, w, w, hall_.get()); - modified = true; - } - if (diff_[z] == bounds_[y] - bounds_[z]) { - PathSet(hall_[y], j + 1, y, hall_.get()); - hall_[y] = j + 1; - } - } - return modified; -} + + ConstPtrArray ovars_; + ConstPtrArray svars_; + scoped_array mins_; + scoped_array maxs_; + const int size_; + RangeBipartiteMatching matching_; +}; } // namespace Constraint* Solver::MakeAllDifferent(const std::vector& vars) { @@ -451,4 +590,11 @@ Constraint* Solver::MakeAllDifferent(const IntVar* const* vars, } } } + +Constraint* Solver::MakeSorted(const std::vector& vars, + const std::vector& sorted) { + CHECK_EQ(vars.size(), sorted.size()); + return RevAlloc(new SortConstraint(this, vars, sorted)); +} + } // namespace operations_research diff --git a/src/constraint_solver/constraint_solver.cc b/src/constraint_solver/constraint_solver.cc index ed3960eca3..8d2128b56c 100644 --- a/src/constraint_solver/constraint_solver.cc +++ b/src/constraint_solver/constraint_solver.cc @@ -2632,6 +2632,7 @@ const char ModelVisitor::kScalProdGreaterOrEqual[] = const char ModelVisitor::kScalProdLessOrEqual[] = "ScalarProductLessOrEqual"; const char ModelVisitor::kSemiContinuous[] = "SemiContinuous"; const char ModelVisitor::kSequenceVariable[] = "SequenceVariable"; +const char ModelVisitor::kSorted[] = "Sorted"; const char ModelVisitor::kSquare[] = "Square"; const char ModelVisitor::kStartExpr[]= "StartExpression"; const char ModelVisitor::kSum[] = "Sum"; diff --git a/src/constraint_solver/constraint_solver.h b/src/constraint_solver/constraint_solver.h index 593970b87a..c61caed688 100644 --- a/src/constraint_solver/constraint_solver.h +++ b/src/constraint_solver/constraint_solver.h @@ -1559,6 +1559,11 @@ class Solver { Constraint* MakeAllDifferent(const IntVar* const* vars, int size, bool stronger_propagation); + // Maintains the relation between an array of variable and its + // sorted counterpart. + Constraint* MakeSorted(const std::vector& vars, + const std::vector& sorted); + // Prevent cycles, nexts variables representing the next in the chain. // Active variables indicate if the corresponding next variable is active; // this could be useful to model unperformed nodes in a routing problem. @@ -3254,6 +3259,7 @@ class ModelVisitor : public BaseObject { static const char kScalProdLessOrEqual[]; static const char kSemiContinuous[]; static const char kSequenceVariable[]; + static const char kSorted[]; static const char kSquare[]; static const char kStartExpr[]; static const char kSum[]; diff --git a/src/constraint_solver/io.cc b/src/constraint_solver/io.cc index 63e5b21615..5c208da8e5 100644 --- a/src/constraint_solver/io.cc +++ b/src/constraint_solver/io.cc @@ -2070,6 +2070,18 @@ SequenceVar* BuildSequenceVariable(CPModelLoader* const builder, return builder->solver()->MakeSequenceVar(vars, proto.name()); } +// ----- kAllDifferent ----- + +Constraint* BuildSorted(CPModelLoader* const builder, + const CPConstraintProto& proto) { + std::vector vars; + VERIFY(builder->ScanArguments(ModelVisitor::kVarsArgument, proto, &vars)); + std::vector targets; + VERIFY(builder->ScanArguments( + ModelVisitor::kTargetArgument, proto, &targets)); + return builder->solver()->MakeSorted(vars, targets); +} + // ----- kSquare ----- IntExpr* BuildSquare(CPModelLoader* const builder, @@ -2601,6 +2613,7 @@ void Solver::InitBuilders() { REGISTER(kScalProdLessOrEqual, BuildScalProdLessOrEqual); REGISTER(kSemiContinuous, BuildSemiContinuous); REGISTER(kSequenceVariable, BuildSequenceVariable); + REGISTER(kSorted, BuildSorted); REGISTER(kSquare, BuildSquare); REGISTER(kStartExpr, BuildStartExpr); REGISTER(kSum, BuildSum);