OR-Tools  8.2
zero_half_cuts.cc
Go to the documentation of this file.
1 // Copyright 2010-2018 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13 
15 
16 namespace operations_research {
17 namespace sat {
18 
19 void ZeroHalfCutHelper::Reset(int size) {
20  rows_.clear();
21  shifted_lp_values_.clear();
22  bound_parity_.clear();
23  col_to_rows_.clear();
24  col_to_rows_.resize(size);
25  tmp_marked_.resize(size);
26 }
27 
29  const std::vector<double>& lp_values,
30  const std::vector<IntegerValue>& lower_bounds,
31  const std::vector<IntegerValue>& upper_bounds) {
32  Reset(lp_values.size());
33 
34  // Shift all variables to their closest bound.
35  lp_values_ = lp_values;
36  for (int i = 0; i < lp_values.size(); ++i) {
37  const double lb_dist = lp_values[i] - ToDouble(lower_bounds[i]);
38  const double ub_dist = ToDouble(upper_bounds[i]) - lp_values_[i];
39  if (lb_dist < ub_dist) {
40  shifted_lp_values_.push_back(lb_dist);
41  bound_parity_.push_back(lower_bounds[i].value() & 1);
42  } else {
43  shifted_lp_values_.push_back(ub_dist);
44  bound_parity_.push_back(upper_bounds[i].value() & 1);
45  }
46  }
47 }
48 
50  // No point pushing an all zero row with a zero rhs.
51  if (binary_row.cols.empty() && !binary_row.rhs_parity) return;
52  for (const int col : binary_row.cols) {
53  col_to_rows_[col].push_back(rows_.size());
54  }
55  rows_.push_back(binary_row);
56 }
57 
59  const glop::RowIndex row,
60  const std::vector<std::pair<glop::ColIndex, IntegerValue>>& terms,
61  IntegerValue lb, IntegerValue ub) {
62  if (terms.size() > kMaxInputConstraintSize) return;
63 
64  double activity = 0.0;
65  IntegerValue magnitude(0);
66  CombinationOfRows binary_row;
67  int rhs_adjust = 0;
68  for (const auto& term : terms) {
69  const int col = term.first.value();
70  activity += ToDouble(term.second) * lp_values_[col];
71  magnitude = std::max(magnitude, IntTypeAbs(term.second));
72 
73  // Only consider odd coefficient.
74  if ((term.second.value() & 1) == 0) continue;
75 
76  // Ignore column in the binary matrix if its lp value is almost zero.
77  if (shifted_lp_values_[col] > 1e-2) {
78  binary_row.cols.push_back(col);
79  }
80 
81  // Because we work on the shifted variable, the rhs needs to be updated.
82  rhs_adjust ^= bound_parity_[col];
83  }
84 
85  // We ignore constraint with large coefficient, since there is little chance
86  // to cancel them and because of that the efficacity of a generated cut will
87  // be limited.
88  if (magnitude > kMaxInputConstraintMagnitude) return;
89 
90  // TODO(user): experiment with the best value. probably only tight rows are
91  // best? and we could use the basis status rather than recomputing the
92  // activity for that.
93  //
94  // TODO(user): Avoid adding duplicates and just randomly pick one. Note
95  // that we should also remove duplicate in a generic way.
96  const double tighteness_threshold = 1e-2;
97  if (ToDouble(ub) - activity < tighteness_threshold) {
98  binary_row.multipliers = {{row, IntegerValue(1)}};
99  binary_row.slack = ToDouble(ub) - activity;
100  binary_row.rhs_parity = (ub.value() & 1) ^ rhs_adjust;
101  AddBinaryRow(binary_row);
102  }
103  if (activity - ToDouble(lb) < tighteness_threshold) {
104  binary_row.multipliers = {{row, IntegerValue(-1)}};
105  binary_row.slack = activity - ToDouble(lb);
106  binary_row.rhs_parity = (lb.value() & 1) ^ rhs_adjust;
107  AddBinaryRow(binary_row);
108  }
109 }
110 
112  std::function<bool(int)> extra_condition, const std::vector<int>& a,
113  std::vector<int>* b) {
114  for (const int v : *b) tmp_marked_[v] = true;
115  for (const int v : a) {
116  if (tmp_marked_[v]) {
117  tmp_marked_[v] = false;
118  } else {
119  tmp_marked_[v] = true;
120 
121  // TODO(user): optim by doing that at the end?
122  b->push_back(v);
123  }
124  }
125 
126  // Remove position that are not marked, and clear tmp_marked_.
127  int new_size = 0;
128  for (const int v : *b) {
129  if (tmp_marked_[v]) {
130  if (extra_condition(v)) {
131  (*b)[new_size++] = v;
132  }
133  tmp_marked_[v] = false;
134  }
135  }
136  b->resize(new_size);
137 }
138 
139 void ZeroHalfCutHelper::ProcessSingletonColumns() {
140  for (const int singleton_col : singleton_cols_) {
141  if (col_to_rows_[singleton_col].empty()) continue;
142  CHECK_EQ(col_to_rows_[singleton_col].size(), 1);
143  const int row = col_to_rows_[singleton_col][0];
144  int new_size = 0;
145  auto& mutable_cols = rows_[row].cols;
146  for (const int col : mutable_cols) {
147  if (col == singleton_col) continue;
148  mutable_cols[new_size++] = col;
149  }
150  CHECK_LT(new_size, mutable_cols.size());
151  mutable_cols.resize(new_size);
152  col_to_rows_[singleton_col].clear();
153  rows_[row].slack += shifted_lp_values_[singleton_col];
154  }
155  singleton_cols_.clear();
156 }
157 
158 // This is basically one step of a Gaussian elimination with the given pivot.
160  int eliminated_row) {
161  CHECK_LE(rows_[eliminated_row].slack, 1e-6);
162  CHECK(!rows_[eliminated_row].cols.empty());
163 
164  // First update the row representation of the matrix.
165  tmp_marked_.resize(std::max(col_to_rows_.size(), rows_.size()));
166  DCHECK(std::all_of(tmp_marked_.begin(), tmp_marked_.end(),
167  [](bool b) { return !b; }));
168  int new_size = 0;
169  for (const int other_row : col_to_rows_[eliminated_col]) {
170  if (other_row == eliminated_row) continue;
171  col_to_rows_[eliminated_col][new_size++] = other_row;
172 
173  SymmetricDifference([](int i) { return true; }, rows_[eliminated_row].cols,
174  &rows_[other_row].cols);
175 
176  // Update slack & parity.
177  rows_[other_row].rhs_parity ^= rows_[eliminated_row].rhs_parity;
178  rows_[other_row].slack += rows_[eliminated_row].slack;
179 
180  // Update the multipliers the same way.
181  {
182  auto& mutable_multipliers = rows_[other_row].multipliers;
183  mutable_multipliers.insert(mutable_multipliers.end(),
184  rows_[eliminated_row].multipliers.begin(),
185  rows_[eliminated_row].multipliers.end());
186  std::sort(mutable_multipliers.begin(), mutable_multipliers.end());
187  int new_size = 0;
188  for (const auto& entry : mutable_multipliers) {
189  if (new_size > 0 && entry == mutable_multipliers[new_size - 1]) {
190  // Cancel both.
191  --new_size;
192  } else {
193  mutable_multipliers[new_size++] = entry;
194  }
195  }
196  mutable_multipliers.resize(new_size);
197  }
198  }
199  col_to_rows_[eliminated_col].resize(new_size);
200 
201  // Then update the col representation of the matrix.
202  //
203  // Note that we remove from the col-wise representation any rows with a large
204  // slack.
205  {
206  int new_size = 0;
207  for (const int other_col : rows_[eliminated_row].cols) {
208  if (other_col == eliminated_col) continue;
209  const int old_size = col_to_rows_[other_col].size();
210  rows_[eliminated_row].cols[new_size++] = other_col;
212  [this](int i) { return rows_[i].slack < kSlackThreshold; },
213  col_to_rows_[eliminated_col], &col_to_rows_[other_col]);
214  if (old_size != 1 && col_to_rows_[other_col].size() == 1) {
215  singleton_cols_.push_back(other_col);
216  }
217  }
218  rows_[eliminated_row].cols.resize(new_size);
219  }
220 
221  // Clear col.
222  col_to_rows_[eliminated_col].clear();
223  rows_[eliminated_row].slack += shifted_lp_values_[eliminated_col];
224 }
225 
226 std::vector<std::vector<std::pair<glop::RowIndex, IntegerValue>>>
228  std::vector<std::vector<std::pair<glop::RowIndex, IntegerValue>>> result;
229 
230  // Initialize singleton_cols_.
231  singleton_cols_.clear();
232  for (int col = 0; col < col_to_rows_.size(); ++col) {
233  if (col_to_rows_[col].size() == 1) singleton_cols_.push_back(col);
234  }
235 
236  // Process rows by increasing size, but randomize if same size.
237  std::vector<int> to_process;
238  for (int row = 0; row < rows_.size(); ++row) to_process.push_back(row);
239  std::shuffle(to_process.begin(), to_process.end(), *random);
240  std::stable_sort(to_process.begin(), to_process.end(), [this](int a, int b) {
241  return rows_[a].cols.size() < rows_[b].cols.size();
242  });
243 
244  for (const int row : to_process) {
245  ProcessSingletonColumns();
246 
247  if (rows_[row].cols.empty()) continue;
248  if (rows_[row].slack > 1e-6) continue;
249  if (rows_[row].multipliers.size() > kMaxAggregationSize) continue;
250 
251  // Heuristic: eliminate the variable with highest shifted lp value.
252  int eliminated_col = -1;
253  double max_lp_value = 0.0;
254  for (const int col : rows_[row].cols) {
255  if (shifted_lp_values_[col] > max_lp_value) {
256  max_lp_value = shifted_lp_values_[col];
257  eliminated_col = col;
258  }
259  }
260  if (eliminated_col == -1) continue;
261 
262  EliminateVarUsingRow(eliminated_col, row);
263  }
264 
265  // As an heuristic, we just try to add zero rows with an odd rhs and a low
266  // enough slack.
267  for (const auto& row : rows_) {
268  if (row.cols.empty() && row.rhs_parity && row.slack < kSlackThreshold) {
269  result.push_back(row.multipliers);
270  }
271  }
272  VLOG(1) << "#candidates: " << result.size() << " / " << rows_.size();
273  return result;
274 }
275 
276 } // namespace sat
277 } // namespace operations_research
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define CHECK_LT(val1, val2)
Definition: base/logging.h:700
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define DCHECK(condition)
Definition: base/logging.h:884
#define CHECK_LE(val1, val2)
Definition: base/logging.h:699
#define VLOG(verboselevel)
Definition: base/logging.h:978
void AddBinaryRow(const CombinationOfRows &binary_row)
void SymmetricDifference(std::function< bool(int)> extra_condition, const std::vector< int > &a, std::vector< int > *b)
void ProcessVariables(const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds)
std::vector< std::vector< std::pair< glop::RowIndex, IntegerValue > > > InterestingCandidates(ModelRandomGenerator *random)
void AddOneConstraint(glop::RowIndex, const std::vector< std::pair< glop::ColIndex, IntegerValue >> &terms, IntegerValue lb, IntegerValue ub)
int64 value
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
double ToDouble(IntegerValue value)
Definition: integer.h:69
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::vector< double > lower_bounds
std::vector< double > upper_bounds
std::vector< std::pair< glop::RowIndex, IntegerValue > > multipliers