Scippy

SCIP

Solving Constraint Integer Programs

double_description_method.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program PolySCIP */
4 /* */
5 /* Copyright (C) 2012-2022 Konrad-Zuse-Zentrum */
6 /* fuer Informationstechnik Berlin */
7 /* */
8 /* PolySCIP is distributed under the terms of the ZIB Academic License. */
9 /* */
10 /* You should have received a copy of the ZIB Academic License */
11 /* along with PolySCIP; see the file LICENCE. */
12 /* */
13 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
14 
15 /**
16  * @file double_description_method.h
17  * @brief Double description method for transforming a polyhedron given
18  * via its h-representation into its v-representation
19  * @author Sebastian Schenker
20  *
21  * The underlying algorithm is described in the paper: "Double description method revisited" by
22  * Komei Fukuda, Alain Prodon
23  */
24 
25 #ifndef POLYSCIP_SRC_DOUBLE_DESCRIPTION_H_INCLUDED
26 #define POLYSCIP_SRC_DOUBLE_DESCRIPTION_H_INCLUDED
27 
28 #include <algorithm>
29 #include <bitset>
30 #include <cstddef>
31 #include <functional> // std::reference_wrapper
32 #include <iostream>
33 #include <memory>
34 #include <unordered_map>
35 #include <ostream>
36 #include <tuple> // std::tie
37 #include <utility>
38 #include <vector>
39 
40 #include "global_functions.h"
41 #include "objscip/objscip.h"
42 #include "PolySCIPConfig.h" // defines MAX_NO_OBJS
43 #include "polyscip_types.h"
44 #include "weight_space_facet.h"
45 
46 namespace polyscip {
47 
48  namespace doubledescription {
49 
50  using H_RepT = std::pair<OutcomeType, ValueType>; ///< Type for element of h-representation
51  using H_RepC = std::vector<H_RepT>; ///< Container for h-representations
52 
53  /**
54  * @class V_RepT
55  * @brief Class for element of v-representation
56  */
57  class V_RepT {
58  public:
59 
60  using SlackMap = std::unordered_map<std::size_t, ValueType>; ///< Maps indices to slacks
61  constexpr static std::size_t kMaxInitialHrepSize = 2*POLYSCIP_MAX_NO_OBJS; ///< Maximal initial size of h-representation
63 
64  /**
65  * Constructor
66  * @param scip SCIP pointer
67  * @param weight Corresponding weight for the constructed vertex
68  * @param wov Corresponding weighted objective value for the constructed vertex
69  * @param h_rep Current h-representation
70  */
71  explicit V_RepT(SCIP* scip,
72  WeightType&& weight,
73  ValueType&& wov,
74  const H_RepC& h_rep);
75 
76  /**
77  * Constructor
78  * @param scip SCIP pointer
79  * @param plus
80  * @param minus
81  * @param index_of_ineq
82  * @param h_rep
83  */
84  explicit V_RepT(SCIP* scip,
85  const V_RepT& plus,
86  const V_RepT& minus,
87  std::size_t index_of_ineq,
88  const H_RepC& h_rep);
89 
90  /**
91  * Equal operator
92  * @param rhs v-representation object to compare with
93  * @return true if objects are considered equal; otherwise false
94  */
95  inline bool operator==(const V_RepT& rhs) const;
96 
97  /**
98  * Not equal operator
99  * @param rhs v-representation object to compare with
100  * @return true if objects are considered not equal; otherwise false
101  */
102  inline bool operator!=(const V_RepT& rhs) const;
103 
104  /**
105  * Print function
106  * @param os Output stream to print to
107  * @param withIncidentFacets Boolean indicating whether incident facets should also be printed
108  * @param h_rep
109  */
110  void print(std::ostream& os,
111  bool withIncidentFacets,
112  const H_RepC& h_rep) const;
113 
114  /**
115  * Indicates whether corresponding weight is neither unit vector nor zero vector
116  * @return true if weight is neither unit vector nor zero vector; false otherwise
117  */
118  bool hasNonUnitAndNonZeroWeight() const;
119 
120  /**
121  * Indicates whether zero indices are subset of given zero indices
122  * @param common_zero_inds Zero indices
123  * @return true if zero indices are subset of given common_zero_inds
124  */
125  bool hasZeroIndsSuperSet(const std::bitset<kMaxInitialHrepSize>& common_zero_inds) const;
126 
127  /**
128  * Get minimal infeasibility index
129  * @return minimal infeasibility index
130  */
131  std::size_t getMinInfeasIndex() const;
132 
133  /**
134  * Indicates whether given index is zero index
135  * @param index Index to check
136  * @return true if index is zero index; false otherwise
137  */
138  bool isZeroSlackIndex(std::size_t index) const;
139 
140  /**
141  * Move weight
142  * @return Rvalue reference to corresponding weight
143  */
144  WeightType&& moveWeight() {return std::move(weight_);};
145 
146  /**
147  * Get weighted objective value
148  * @return weighted objective value
149  */
150  ValueType getWov() const {return wov_;};
151 
152  private:
153  constexpr static double kNormalizingThreshold = 1e+5; ///< Normalizing threshold
154 
155  /**
156  * Indicates whether corresponding weight should be normalized
157  * @param threshold Normalizing threshold
158  * @return true if weight should be normalized; false otherwise
159  */
160  bool shouldNormalize(double threshold) const;
161 
162  /**
163  * Normalizing weight function
164  * @param normalizing_val Value by which weights are normalized
165  */
166  void normalize(double normalizing_val);
167 
168  /**
169  * Get slack value
170  * @param index Corresponding index to get slack for
171  * @return Slack value
172  */
173  ValueType getSlack(std::size_t index) const;
174 
175  /**
176  * Set slack and minimal infeasibility index
177  * @param scip SCIP pointer
178  * @param h_rep h-representation
179  */
180  void setSlacksAndMinInfeasInd(SCIP* scip,
181  const H_RepC& h_rep);
182 
183  WeightType weight_; ///< Corresponding weight of v-representation
184  ValueType wov_; ///< Corresponding weighted objective value of v-representation
185  std::pair<bool, std::size_t> min_infeas_ind_; ///< Minimal infeasibility index
186  SlackMap inds_to_slacks_; ///< Maps index to slack
187  std::bitset<kMaxInitialHrepSize> zero_slacks_; ///< Indices which have zero slack
188  };
189 
190  using V_RepC = std::vector<std::shared_ptr<V_RepT>>; ///< Container for v-representations
191 
192  /**
193  * @class DoubleDescriptionMethod
194  * @brief Algorithm for transforming h-representation to v-representation
195  * @details Details of the underlying algorithm can be found in "Double Description Method Revisited" by K. Fukuda and A. Prodon
196  */
198  public:
199  using AdjPair = std::pair<std::reference_wrapper<const V_RepT>, std::reference_wrapper<const V_RepT>>; ///< Pair of adjacent v-representation elements
200  using AdjPairContainer = std::vector<AdjPair>; ///< Container for adjacent pairs
201 
202  /**
203  * Constructor
204  * @param scip SCIP pointer
205  * @param no_all_obj Number of objectives
206  * @param bounded Bounded results
207  * @param unbounded Unbounded results
208  */
210  std::size_t no_all_obj,
211  const ResultContainer &bounded,
212  const ResultContainer &unbounded);
213 
214  /**
215  * Standard double description algorithm
216  */
217  void computeVRep();
218 
219  /**
220  * Variation 1 of double description algorithm
221  */
222  void computeVRep_Var1();
223 
224  /**
225  * Print function
226  * @param os Output stream to write to
227  * @param withIncidentFacets Indicates whether corresponding facets should be printed
228  */
229  void printVRep(std::ostream &os = std::cout,
230  bool withIncidentFacets = false) const;
231 
232  /**
233  * Get size of v-representation
234  * @return Number of elements in v-representation
235  */
236  std::size_t size() const { return v_rep_.size(); };
237 
238  /**
239  * Copy entire h-representation
240  * @return Copy of container storing entire h-representation
241  */
242  H_RepC getHRep() {return h_rep_;};
243 
244  /**
245  * Copy entire v-representation
246  * @return Copy of container storing entire v-representation
247  */
248  V_RepC getVRep() {return v_rep_;};
249 
250  /**
251  * Make entire h-representation movable
252  * @return Rvalue reference to container storing entire h-representation
253  */
254  H_RepC&& moveHRep() {return std::move(h_rep_);};
255 
256  /**
257  * Make entire v-representation movable
258  * @return Rvalue reference to container storing entire v-representation
259  */
260  V_RepC&& moveVRep() {return std::move(v_rep_);};
261 
262  private:
263  using SlackContainer = std::vector<std::vector<std::size_t>>; ///< Container for slacks
264 
265  /**
266  * Variable orders
267  */
268  enum class VarOrder {
269  keep_var_order, ///< Keep variable order
270  change_var_order ///< Exchange variable order
271  };
272 
273  /**
274  * Computes initial v-representation of following h-representation:
275  * 1) bounded \\cdot (w_1,...,w_k) -a >= 0
276  * 2) w_1 >= 0
277  * ...
278  * k+1) w_k >= 0
279  * where (w_1,...,w_k) is a weight vector
280  * @return Initial v-representation
281  */
282  V_RepC computeInitialVRep() const;
283 
284  /**
285  * Computes the intersection of zero indices for two given v-representation elements
286  * @param v First element
287  * @param w Second element
288  * @return bitset with ones indicating common zero slacks
289  */
290  std::bitset<V_RepT::kMaxInitialHrepSize> getCommonZeroSlackIndices(const V_RepT &v,
291  const V_RepT &w) const;
292 
293  /**
294  * Check for minimal infeasibility condition
295  * @param r1 First element
296  * @param r2 Second element
297  * @return Tuple
298  */
299  std::tuple<bool, VarOrder, std::size_t> minInfeasCondition(const V_RepT& r1,
300  const V_RepT& r2) const;
301 
302  /**
303  * Applies infeasibility condition
304  * @param r1
305  * @param r2
306  * @param current_v_rep
307  * @param index
308  * @param with_adjacency_test
309  */
310  void applyInfeasCondition(const V_RepT& r1,
311  const V_RepT& r2,
312  const V_RepC& current_v_rep,
313  std::size_t index,
314  bool with_adjacency_test);
315 
316 
317  /**
318  * Detailed description of algorithm can be found on page 13 in "Double Description Method Revisited"
319  * @param r1
320  * @param r2
321  * @param k
322  * @param i
323  * @param v_rep
324  * @param with_adjacency_test
325  */
326  void conditionalStoreEdge(const V_RepT& r1,
327  const V_RepT& r2,
328  std::size_t k,
329  std::size_t i,
330  const V_RepC& v_rep,
331  bool with_adjacency_test);
332 
333 
334  /**
335  * Indicates whether given elements are adjacent
336  * @param ray1
337  * @param ray2
338  * @param v_rep
339  * @return True if ray1 and ray2 are adjancent; false otherwise
340  */
341  bool rayPairIsAdjacent(const V_RepT& ray1,
342  const V_RepT& ray2,
343  const V_RepC& v_rep) const;
344 
345  /**
346  * Indicates whether first element is multiple of second element
347  * @param v First v-representation element
348  * @param w Second v-representation element
349  * @return true if v is multiple of w; false otherwise
350  */
351  bool isMultiple(const V_RepT& v, const V_RepT& w) const;
352 
353  /**
354  * Indicates whether weight of first element is a multiple (with value v_multiple) of weight of second element
355  * @param scip SCIP pointer
356  * @param v_multiple Value
357  * @param v First element
358  * @param w Second element
359  * @return true if weight of first element is a v_multiple of weight of second element; false otherwise
360  */
361  bool weightIsMultiple(SCIP* scip,
362  double v_multiple,
363  const V_RepT& v,
364  const V_RepT& w) const;
365 
366  /**
367  * Compute adjacent pairs of elements in v-representation
368  * @param plus Elements in plus partition
369  * @param minus Elements in minus partition
370  * @param current_v_rep Current v-representation
371  * @return Adjacent elements
372  */
373  AdjPairContainer computeAdjacentPairs(const V_RepC& plus,
374  const V_RepC& minus,
375  const V_RepC& current_v_rep) const;
376 
377  /**
378  * Incorporate next hyperplane into v-representation; standard algorithm
379  * @param cur_v_rep Current v-representation
380  * @return Extended v-representation
381  */
382  V_RepC extendVRep(V_RepC&& cur_v_rep);
383 
384  /**
385  * Incorporate next hyperplane into v-representation; Variation1 algorithm
386  * @param current_v_rep Current v-representation
387  * @return Extended v-representation
388  */
389  V_RepC extendVRep_Var1(V_RepC&& current_v_rep);
390 
391 
392  SCIP *scip_; ///< SCIP pointer
393  std::size_t outcome_dimension_; ///< Dimension of outcome in objective space
394  std::size_t current_hrep_index_; ///< hyperplane currently considered
395  H_RepC h_rep_; ///< h-representation
396  V_RepC v_rep_; ///< v-representation
397  std::vector< AdjPairContainer > adj_pairs_; ///< Adjacent pairs
398  };
399  }
400 }
401 #endif //POLYSCIP_SRC_DOUBLE_DESCRIPTION_H_INCLUDED
Some preprocessor directives.
bool hasZeroIndsSuperSet(const std::bitset< kMaxInitialHrepSize > &common_zero_inds) const
Class for element of v-representation.
std::vector< ValueType > WeightType
Type for weight vectors.
std::vector< AdjPair > AdjPairContainer
Container for adjacent pairs.
V_RepT(SCIP *scip, WeightType &&weight, ValueType &&wov, const H_RepC &h_rep)
SCIP_Real ValueType
Type for computed values.
void print(std::ostream &os, bool withIncidentFacets, const H_RepC &h_rep) const
static constexpr std::size_t kMaxInitialHrepSize
Maximal initial size of h-representation.
std::vector< std::shared_ptr< V_RepT > > V_RepC
Container for v-representations.
std::vector< H_RepT > H_RepC
Container for h-representations.
General types used for PolySCIP.
SCIP_VAR * w
Definition: circlepacking.c:58
bool isZeroSlackIndex(std::size_t index) const
C++ wrapper classes for SCIP.
#define POLYSCIP_MAX_NO_OBJS
Maximal number of objectives considered by PolySCIP.
std::unordered_map< std::size_t, ValueType > SlackMap
Maps indices to slacks.
std::pair< std::reference_wrapper< const V_RepT >, std::reference_wrapper< const V_RepT > > AdjPair
Pair of adjacent v-representation elements.
std::vector< Result > ResultContainer
Container type for results.
Global helper functions.
std::pair< OutcomeType, ValueType > H_RepT
Type for element of h-representation.
Class representing a facet of the weight space polyhedron.
Algorithm for transforming h-representation to v-representation.