42 using std::make_shared;
48 namespace doubledescription {
60 min_infeas_ind_(false, 0)
62 if (shouldNormalize(kNormalizingThreshold))
63 normalize(kNormalizingThreshold);
65 setSlacksAndMinInfeasInd(scip, h_rep);
79 std::size_t index_of_ineq,
81 : min_infeas_ind_(false, 0)
83 auto m_coeff = plus.getSlack(index_of_ineq);
84 auto p_coeff = minus.getSlack(index_of_ineq);
88 std::transform(minus.weight_.cbegin(), minus.weight_.cend(),
89 plus.weight_.cbegin(), std::back_inserter(weight_),
91 return m_coeff * m_val - p_coeff * p_val;
93 wov_ = m_coeff*minus.wov_ - p_coeff*plus.wov_;
95 if (shouldNormalize(kNormalizingThreshold))
96 normalize(kNormalizingThreshold);
98 setSlacksAndMinInfeasInd(scip, h_rep);
107 return std::tie(weight_, wov_) == std::tie(rhs.weight_, rhs.wov_);
124 ValueType V_RepT::getSlack(std::size_t index)
const {
125 assert (inds_to_slacks_.count(index) != 0);
126 return inds_to_slacks_.at(index);
134 assert (min_infeas_ind_.first);
135 return min_infeas_ind_.second;
144 return zero_slacks_.test(index);
152 return (std::count_if(weight_.cbegin(), weight_.cend(), [](
const ValueType& val){
return val != 0.;}) > 1);
161 void V_RepT::setSlacksAndMinInfeasInd(
SCIP*
scip,
const H_RepC& h_rep) {
162 for (
size_t i=0; i<h_rep.size(); ++i) {
163 auto result = std::inner_product(weight_.cbegin(), weight_.cend(),
164 h_rep[i].first.cbegin(), -(h_rep[i].second * wov_));
166 inds_to_slacks_.emplace(i, result);
168 inds_to_slacks_[i] = 0.;
169 zero_slacks_.set(i, 1);
171 else if (
SCIPisNegative(scip, result) && !min_infeas_ind_.first) {
172 min_infeas_ind_.first =
true;
173 min_infeas_ind_.second = i;
176 if (!min_infeas_ind_.first) {
177 min_infeas_ind_.first =
true;
178 min_infeas_ind_.second = h_rep.size();
187 bool V_RepT::shouldNormalize(
double threshold)
const {
190 else if (std::any_of(weight_.cbegin(), weight_.cend(), [threshold](
ValueType w){
return w > threshold;}))
200 void V_RepT::normalize(
double normalizing_val) {
201 std::transform(begin(weight_), end(weight_), begin(weight_),
202 [normalizing_val](
const ValueType &val) {
return val / normalizing_val; });
203 wov_ /= normalizing_val;
214 os <<
" Coeff = " << wov_ <<
"\n";
215 os <<
"MinInfeasIndex = " << min_infeas_ind_.second <<
"\n";
216 if (withIncidentFacets) {
219 if (zero_slacks_[i]) {
221 os <<
" " << h_rep[i].second <<
"\n";
235 return (common_zero_inds & zero_slacks_).count() >= common_zero_inds.count();
251 outcome_dimension_(no_all_objs),
252 current_hrep_index_(no_all_objs),
253 adj_pairs_(no_all_objs+bounded_results.size()+unbounded_results.size(),
AdjPairContainer{})
256 assert (!bounded_results.empty());
258 for (
size_t i=0; i<outcome_dimension_; ++i) {
259 auto unit_vec =
OutcomeType(outcome_dimension_, 0.);
261 h_rep_.emplace_back(std::move(unit_vec), 0.);
264 for (
const auto &bd : bounded_results) {
265 h_rep_.emplace_back(bd.second, 1.);
267 for (
const auto &unbd : unbounded_results) {
268 h_rep_.emplace_back(unbd.second, 0.);
278 std::tuple<bool, DoubleDescriptionMethod::VarOrder, size_t> DoubleDescriptionMethod::minInfeasCondition(
const V_RepT& r1,
const V_RepT& r2)
const {
282 assert (current_hrep_index_ < minInfeasInd1 && minInfeasInd1 < h_rep_.size());
283 return std::make_tuple(
true, VarOrder::change_var_order, minInfeasInd1);
285 else if (minInfeasInd2 < minInfeasInd1 && !r1.
isZeroSlackIndex(minInfeasInd2)) {
286 assert (current_hrep_index_ < minInfeasInd2 && minInfeasInd2 < h_rep_.size());
287 return std::make_tuple(
true, VarOrder::keep_var_order, minInfeasInd2);
290 return std::make_tuple(
false, VarOrder::keep_var_order, 0);
303 void DoubleDescriptionMethod::conditionalStoreEdge(
const V_RepT& plus,
308 bool with_adjacency_test) {
310 for (
auto index=i+1; index<k; ++index) {
314 if (!with_adjacency_test || rayPairIsAdjacent(plus, minus, v_rep)) {
315 adj_pairs_.at(k).push_back({std::cref(plus), std::cref(minus)});
325 for (
const auto& v : v_rep_)
326 v->print(os, withIncidentFacets, h_rep_);
333 auto current_v_rep = computeInitialVRep();
334 ++current_hrep_index_;
335 while (current_hrep_index_ < h_rep_.size()) {
336 current_v_rep = extendVRep(std::move(current_v_rep));
337 ++current_hrep_index_;
339 v_rep_ = current_v_rep;
350 void DoubleDescriptionMethod::applyInfeasCondition(
const V_RepT& r1,
354 bool with_adjacency_test
356 auto infeasTuple = minInfeasCondition(r1, r2);
357 if (std::get<0>(infeasTuple)) {
358 if (std::get<1>(infeasTuple) == VarOrder::keep_var_order){
359 conditionalStoreEdge(r1, r2, std::get<2>(infeasTuple), index, v_rep, with_adjacency_test);
362 conditionalStoreEdge(r2, r1, std::get<2>(infeasTuple), index, v_rep, with_adjacency_test);
371 auto current_v_rep = computeInitialVRep();
373 for (
auto r1_it=begin(current_v_rep); r1_it!=std::prev(end(current_v_rep)); ++r1_it) {
374 for (
auto r2_it=std::next(r1_it); r2_it!=end(current_v_rep); ++r2_it) {
375 const V_RepT& r1 = r1_it->operator*();
376 const V_RepT& r2 = r2_it->operator*();
377 applyInfeasCondition(r1, r2, current_v_rep, current_hrep_index_,
true);
381 ++current_hrep_index_;
382 while (current_hrep_index_ < h_rep_.size()) {
383 current_v_rep = extendVRep_Var1(std::move(current_v_rep));
384 ++current_hrep_index_;
386 v_rep_ = current_v_rep;
395 V_RepC DoubleDescriptionMethod::extendVRep_Var1(
V_RepC&& current_v_rep) {
397 auto extended_v_rep =
V_RepC{};
398 const H_RepT&
h = h_rep_[current_hrep_index_];
400 for (
const auto& v : current_v_rep) {
401 auto result = std::inner_product(v->weight_.cbegin(), v->weight_.cend(),
402 h.first.cbegin(), -(v->wov_ * h.second));
404 extended_v_rep.push_back(v);
407 extended_v_rep.push_back(v);
411 for (
auto p : adj_pairs_[current_hrep_index_]) {
412 auto new_ray = make_shared<V_RepT>(scip_, p.first.get(), p.second.get(), current_hrep_index_, h_rep_);
413 extended_v_rep.push_back(new_ray);
414 applyInfeasCondition(*new_ray, p.first.get(), extended_v_rep, current_hrep_index_,
false);
417 for (
auto r1_it=extended_v_rep.cbegin(); r1_it!=std::prev(extended_v_rep.cend()); ++r1_it) {
418 for (
auto r2_it=std::next(r1_it); r2_it!=extended_v_rep.cend(); ++r2_it) {
419 const V_RepT& r1 = r1_it->operator*();
420 const V_RepT& r2 = r2_it->operator*();
422 applyInfeasCondition(r1, r2, extended_v_rep, current_hrep_index_,
true);
427 return extended_v_rep;
435 V_RepC DoubleDescriptionMethod::extendVRep(
V_RepC&& cur_v_rep) {
436 auto extended_v_rep =
V_RepC{};
439 const H_RepT& hrep = h_rep_[current_hrep_index_];
441 for (
const auto& v : cur_v_rep) {
442 auto result = std::inner_product(v->weight_.cbegin(), v->weight_.cend(),
443 hrep.first.cbegin(), -(v->wov_ * hrep.second));
448 extended_v_rep.push_back(std::move(v));
456 auto adj_pairs = computeAdjacentPairs(plus, minus, cur_v_rep);
457 for (
const auto& p : adj_pairs)
458 extended_v_rep.push_back(make_shared<V_RepT>(scip_, p.first.get(), p.second.get(), current_hrep_index_, h_rep_));
460 for (
const auto& p : plus)
461 extended_v_rep.push_back(std::move(p));
463 return extended_v_rep;
476 const V_RepC& current_v_rep)
const {
478 for (
const auto& p_ptr : plus) {
479 for (
const auto& m_ptr : minus) {
480 assert (*p_ptr != *m_ptr);
481 if (rayPairIsAdjacent(*p_ptr, *m_ptr, current_v_rep))
482 adj_pairs.push_back({std::cref(*p_ptr), std::cref(*m_ptr)});
494 std::bitset<V_RepT::kMaxInitialHrepSize> DoubleDescriptionMethod::getCommonZeroSlackIndices(
const V_RepT &v,
496 return v.zero_slacks_ & w.zero_slacks_;
507 bool DoubleDescriptionMethod::rayPairIsAdjacent(
const V_RepT& ray1,
509 const V_RepC& v_rep)
const {
510 auto common_zero_inds = getCommonZeroSlackIndices(ray1, ray2);
512 for (
const auto& ray_sptr : v_rep) {
513 if (ray1 == *ray_sptr || ray2 == *ray_sptr) {
516 else if (ray_sptr->hasZeroIndsSuperSet(common_zero_inds)) {
518 if (!isMultiple(*ray_sptr, ray1) && !isMultiple(*ray_sptr, ray2))
533 bool DoubleDescriptionMethod::isMultiple(
const V_RepT& v,
const V_RepT& w)
const {
534 assert (v.weight_.size() == w.weight_.size());
535 auto scip_ptr = scip_;
536 if (
SCIPisEQ(scip_, v.wov_, w.wov_)) {
537 auto mismatch_pair = std::mismatch(begin(v.weight_),
541 {return SCIPisEQ(scip_ptr, v_val, w_val);});
542 if (mismatch_pair.first == end(v.weight_)) {
547 auto v_val = *mismatch_pair.first;
548 auto w_val = *mismatch_pair.second;
553 auto multiple = w_val / v_val;
554 return weightIsMultiple(scip_ptr, multiple, v, w);
564 assert (!std::all_of(begin(v.weight_), end(v.weight_), [scip_ptr](
ValueType val){return SCIPisZero(scip_ptr,val);}));
568 assert (!std::all_of(begin(w.weight_), end(w.weight_), [scip_ptr](
ValueType val){return SCIPisZero(scip_ptr,val);}));
572 auto multiple = w.wov_ / v.wov_;
573 return weightIsMultiple(scip_ptr, multiple, v, w);
586 bool DoubleDescriptionMethod::weightIsMultiple(
SCIP* scip,
double v_multiple,
const V_RepT& v,
const V_RepT& w)
const {
587 auto mismatch_weight = std::mismatch(begin(v.weight_),
591 {return SCIPisEQ(scip, v_multiple * v_val, w_val);});
592 return mismatch_weight.first == end(v.weight_);
605 V_RepC DoubleDescriptionMethod::computeInitialVRep()
const {
607 auto init_v_rep =
V_RepC{};
608 init_v_rep.push_back(make_shared<V_RepT>(scip_,
WeightType(outcome_dimension_, 0.), -1., h_rep_));
609 for (
size_t i=0; i<current_hrep_index_; ++i) {
610 auto unit_vec =
WeightType(outcome_dimension_, 0.);
612 auto wov = h_rep_[current_hrep_index_].first.at(i);
613 init_v_rep.push_back(make_shared<V_RepT>(scip_, std::move(unit_vec), std::move(wov), h_rep_));
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Some preprocessor directives.
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
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.
bool hasNonUnitAndNonZeroWeight() const
std::size_t getMinInfeasIndex() const
void printVRep(std::ostream &os=std::cout, bool withIncidentFacets=false) const
std::vector< std::shared_ptr< V_RepT > > V_RepC
Container for v-representations.
bool operator==(const V_RepT &rhs) const
std::vector< H_RepT > H_RepC
Container for h-representations.
General types used for PolySCIP.
std::vector< ValueType > OutcomeType
Type for points, rays in objective space.
bool isZeroSlackIndex(std::size_t index) const
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
void print(const Container &container, const std::string &prefix="", const std::string &suffix="", std::ostream &os=std::cout, bool negate=false, int prec=6)
C++ wrapper classes for SCIP.
SCIPInterval fabs(const SCIPInterval &x)
DoubleDescriptionMethod(SCIP *scip, std::size_t no_all_obj, const ResultContainer &bounded, const ResultContainer &unbounded)
bool operator!=(const V_RepT &rhs) const
std::vector< Result > ResultContainer
Container type for results.
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
std::pair< OutcomeType, ValueType > H_RepT
Type for element of h-representation.
Double description method for transforming a polyhedron given via its h-representation into its v-rep...