this repo has no description
dekker.one
1/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2/*
3 * Main authors:
4 * Samuel Gagnon <samuel.gagnon92@gmail.com>
5 *
6 * Copyright:
7 * Samuel Gagnon, 2018
8 *
9 * This file is part of Gecode, the generic constraint
10 * development environment:
11 * http://www.gecode.org
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining
14 * a copy of this software and associated documentation files (the
15 * "Software"), to deal in the Software without restriction, including
16 * without limitation the rights to use, copy, modify, merge, publish,
17 * distribute, sublicense, and/or sell copies of the Software, and to
18 * permit persons to whom the Software is furnished to do so, subject to
19 * the following conditions:
20 *
21 * The above copyright notice and this permission notice shall be
22 * included in all copies or substantial portions of the Software.
23 *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
27 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
28 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
29 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
30 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31 *
32 */
33
34#ifdef GECODE_HAS_CBS
35
36#include <limits>
37#include <algorithm>
38
39namespace Gecode { namespace Int { namespace Distinct {
40
41 /**
42 * \brief Minc and Brégman factors
43 *
44 * Factors precomputed for every value in the domain of x. Thoses factors are
45 * used to compute the Minc and Brégman upper bound for the permanent in
46 * counting base search
47 */
48 const int MAX_MINC_FACTORS = 400;
49 extern const double mincFactors[MAX_MINC_FACTORS];
50
51 forceinline double
52 getMincFactor(int domSize) {
53 return mincFactors[domSize - 1];
54 }
55
56 /**
57 * \brief Liang and Bai factors
58 *
59 * Factors precomputed for every index and domain size in x. Thoses factors
60 * are used to compute the Liang and Bai upper bound for the permanent in
61 * counting base search
62 */
63 const int WIDTH_LIANG_BAI_FACTORS = 400;
64 extern const double liangBaiFactors[WIDTH_LIANG_BAI_FACTORS * WIDTH_LIANG_BAI_FACTORS];
65
66 forceinline double
67 getLiangBaiFactor(int index, int domSize) {
68 return liangBaiFactors[index*WIDTH_LIANG_BAI_FACTORS + domSize - 1];
69 }
70
71 /**
72 * \brief Mapping of each value to its permanent update
73 *
74 */
75 class ValToUpdate {
76 private:
77 /// Minimum value of the union of all variable domains in the propagator
78 const int minVal;
79 /// Minc and Brégman estimation update for each value
80 double* mincUpdate;
81 /// Liang and Bai estimation update for each value
82 double* liangUpdate;
83 public:
84 template<class View>
85 ValToUpdate(const ViewArray<View>& x,
86 int minDomVal, int maxDomVal, Region& r);
87 /**
88 * Gives the update we have to apply to the Minc and Brégman
89 * estimation of the permanent if we fix a variable of cardinalty
90 * \a varSize to the value \a val.
91 */
92 double getMincUpdate(int val, unsigned int varSize) const;
93 /**
94 * Gives the update we have to apply to the Liang and Bai
95 * estimation of the
96 * permanent if we fix a variable of cardinalty \a varSize
97 * to the value "val".
98 */
99 double getLiangUpdate(int val, unsigned int idx, unsigned int varSize) const;
100 };
101
102 template<class View>
103 forceinline
104 ValToUpdate::ValToUpdate(const ViewArray<View>& x,
105 int minDomVal, int maxDomVal, Region& r)
106 : minVal(minDomVal) {
107 unsigned int width = maxDomVal - minDomVal + 1;
108 mincUpdate = r.alloc<double>(width);
109 std::fill(mincUpdate, mincUpdate + width, 1);
110 liangUpdate = r.alloc<double>(width);
111 std::fill(liangUpdate, liangUpdate + width, 1);
112
113 for (int i=0; i<x.size(); i++) {
114 if (x[i].assigned()) continue;
115 size_t s = x[i].size();
116 for (ViewValues<View> val(x[i]); val(); ++val) {
117 int idx = val.val() - minVal;
118 mincUpdate[idx] *= getMincFactor(s-1) / getMincFactor(s);
119 liangUpdate[idx] *= getLiangBaiFactor(i, s-1) / getLiangBaiFactor(i, s);
120 }
121 }
122 }
123
124 forceinline double
125 ValToUpdate::getMincUpdate(int val, unsigned int varSize) const {
126 return mincUpdate[val-minVal] / getMincFactor(varSize-1);
127 }
128
129 forceinline double
130 ValToUpdate::getLiangUpdate(int val, unsigned int idx,
131 unsigned int varSize) const {
132 return liangUpdate[val-minVal] / getLiangBaiFactor(idx, varSize-1);
133 }
134
135
136 template<class View>
137 void cbsdistinct(Space&, unsigned int prop_id, const ViewArray<View>& x,
138 Propagator::SendMarginal send) {
139 // Computation of Minc and Brégman and Liang and Bai upper bounds for
140 // the permanent of the whole constraint
141 struct UB {
142 double minc;
143 double liangBai;
144 };
145
146 UB ub{1,1};
147 for (int i=0; i<x.size(); i++) {
148 unsigned int s = x[i].size();
149 if ((s >= MAX_MINC_FACTORS) || (s >= WIDTH_LIANG_BAI_FACTORS))
150 throw Gecode::Exception("Int::Distinct::cbsdistinct",
151 "Variable cardinality too big for using counting-based"
152 "search with distinct constraints");
153 ub.minc *= getMincFactor(s);
154 ub.liangBai *= getLiangBaiFactor(i, s);
155 }
156
157 // Minimum and maximum value of the union of all variable domains
158 int minVal = std::numeric_limits<int>::max();
159 int maxVal = std::numeric_limits<int>::min();
160 for (const auto& v : x) {
161 if (v.assigned()) continue;
162 minVal = std::min(v.min(), minVal);
163 maxVal = std::max(v.max(), maxVal);
164 }
165
166 // For each possible value, we compute the update we have to apply to the
167 // permanent of the whole constraint to get the new solution count
168 Region r;
169 ValToUpdate valToUpdate(x, minVal, maxVal, r);
170
171 // Preallocated memory for holding solution counts for all values of a
172 // variable during computation
173 double* solCounts = r.alloc<double>(maxVal - minVal + 1);
174
175 for (int i=0; i<x.size(); i++) {
176 if (x[i].assigned()) continue;
177
178 // Normalization constant for keeping densities values between 0 and 1
179 double normalization = 0;
180 // We calculate the density for every possible value assignment
181 for (ViewValues<View> val(x[i]); val(); ++val) {
182 UB localUB = ub;
183 int v = val.val();
184 unsigned int s = x[i].size();
185
186 // We update both upper bounds according to the assigned value, yielding
187 // two new estimations for the upper bound
188 localUB.minc *= valToUpdate.getMincUpdate(v, s);
189 localUB.liangBai *= valToUpdate.getLiangUpdate(v, i, s);
190
191 // We take the lower upper bound as our estimation for the permanent
192 double lowerUB = std::min(localUB.minc, ::sqrt(localUB.liangBai));
193 solCounts[val.val() - minVal] = lowerUB;
194 normalization += lowerUB;
195 }
196
197 // Because we approximate the permanent of each value for the variable, we
198 // assign densities in a separate loop where we normalize solution densities.
199 for (ViewValues<View> val(x[i]); val(); ++val) {
200 // In practice, send is going to be a function provided by a brancher.
201 // Thus, the brancher will receive each computed solution densities via
202 // this call. For more details, please see Section 4 of the dissertation
203 // "Improvement and Integration of Counting-Based Search Heuristics in
204 // Constraint Programming" by Samuel Gagnon.
205 send(prop_id,
206 x[i].id(),
207 x[i].baseval(val.val()),
208 solCounts[val.val() - minVal] / normalization);
209 }
210 }
211 }
212
213 template<class View>
214 void cbssize(const ViewArray<View>& x, Propagator::InDecision in,
215 unsigned int& size, unsigned int& size_b) {
216 size = 0;
217 size_b = 0;
218 for (const auto& v : x) {
219 if (!v.assigned()) {
220 size += v.size();
221 if (in(v.id()))
222 size_b += v.size();
223 }
224 }
225 }
226
227}}}
228
229#endif
230
231// STATISTICS: int-prop
232