Scippy

SCIP

Solving Constraint Integer Programs

polyscip.cpp
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 polyscip.cpp
17  * @brief Implements PolySCIP solver class
18  * @author Sebastian Schenker
19  *
20  */
21 
22 
23 #include "polyscip.h"
24 
25 #include <algorithm> //std::transform, std::max, std::copy
26 #include <array>
27 #include <cmath> //std::abs, std::modf
28 #include <cstddef> //std::size_t
29 #include <fstream>
30 #include <functional> //std::plus, std::negate, std::function, std::reference_wrapper
31 #include <iomanip> //std::set_precision
32 #include <iostream>
33 #include <iterator> //std::advance, std::back_inserter
34 #include <limits>
35 #include <list>
36 #include <ostream>
37 #include <map>
38 #include <memory> //std::addressof, std::unique_ptr
39 #include <numeric> //std::inner_product
40 #include <string>
41 #include <type_traits> //std::remove_const
42 #include <utility> //std::make_pair
43 #include <vector>
44 
46 #include "scip/scip.h"
48 #include "cmd_line_args.h"
49 #include "global_functions.h"
50 #include "polyscip_types.h"
51 #include "prob_data_objectives.h"
52 #include "ReaderMOP.h"
54 
55 using std::addressof;
56 using std::array;
57 using std::begin;
58 using std::cout;
59 using std::end;
60 using std::list;
61 using std::make_pair;
62 using std::map;
63 using std::max;
64 using std::ostream;
65 using std::pair;
66 using std::reference_wrapper;
67 using std::size_t;
68 using std::string;
69 using std::vector;
70 
71 namespace polyscip {
72 
74 
75  /**
76  * Default constructor
77  * @param outcome Corresponding outcome to take two-dimensional projection of
78  * @param first First (objective) index of outcome to consider for projection
79  * @param second Second (objective) index of outcome to consider for projection
80  */
81  TwoDProj::TwoDProj(const OutcomeType& outcome, size_t first, size_t second)
82  : proj_(outcome.at(first), outcome.at(second))
83  {}
84 
85  /**
86  * Ostream operator
87  * @param os Output stream
88  * @param p Projection to write to stream
89  * @return Output stream
90  */
91  ostream& operator<<(ostream& os, const TwoDProj& p) {
92  os << "Proj = [" << p.proj_.first << ", " << p.proj_.second << "]";
93  return os;
94  }
95 
96  /**
97  * Ostream operator
98  * @param os Output stream
99  * @param nd Non-dominated projections to write to stream
100  * @return Output stream
101  */
102  ostream& operator<<(ostream& os, const NondomProjections& nd) {
103  os << "Nondominated projections: ";
104  for (const auto& p_pair : nd.nondom_projections_)
105  os << p_pair.first << " ";
106  return os;
107  }
108 
109  /**
110  * Add projection and corresponding result to non-dominated projections
111  * @param proj Projection to add
112  * @param res Corresponding result of projections
113  * @return Iterator pointing to proj
114  */
115  NondomProjections::ProjMap::iterator NondomProjections::add(TwoDProj proj, Result res) {
116  auto ret_find = nondom_projections_.find(proj);
117  if (ret_find == end(nondom_projections_)) { // key not found
118  auto ret = nondom_projections_.emplace(std::move(proj), ResultContainer{std::move(res)});
119  return ret.first;
120  }
121  else { // key found
122  nondom_projections_[proj].push_back(std::move(res));
123  return ret_find;
124  }
125  }
126 
127  /**
128  * Default constructor
129  * @param epsilon Error value for comparisons
130  * @param supported Results to take non-dominated projections
131  * @param first First (objective) index to consider for projection
132  * @param second Second (objective) index to consider for projection
133  */
135  const ResultContainer &supported,
136  size_t first,
137  size_t second)
138  : epsilon_(eps),
139  nondom_projections_([eps](const TwoDProj& lhs, const TwoDProj& rhs){
140  if (lhs.getFirst() + eps < rhs.getFirst())
141  return true;
142  else if (rhs.getFirst() + eps < lhs.getFirst())
143  return false;
144  else
145  return lhs.getSecond() < rhs.getSecond();
146  })
147  {
148  assert (first < second);
149  assert (!supported.empty());
150  for (const auto& res : supported) {
151  add(TwoDProj(res.second, first, second), res);
152  }
153 
154  auto it = begin(nondom_projections_);
155  while (std::next(it) != end(nondom_projections_)) {
156  auto next = std::next(it);
157  if (epsilonDominates(it->first, next->first)) {
158  nondom_projections_.erase(next);
159  }
160  else {
161  ++it;
162  }
163  }
164  assert (!nondom_projections_.empty());
165  current_ = begin(nondom_projections_);
166  }
167 
168  /**
169  * lhs-Projection epsilonDominates rhs-Projection if lhs.first - epsilon < rhs.first && lhs.second - epsilon < rhs.second
170  * @param lhs lhs-Projection
171  * @param rhs rhs-Projection
172  * @return true if lhs-Projection epsilon-dominated rhs-Projection; false otherwise
173  */
174  bool NondomProjections::epsilonDominates(const TwoDProj& lhs, const TwoDProj& rhs) const {
175  return lhs.getFirst() - epsilon_ < rhs.getFirst() && lhs.getSecond() - epsilon_ < rhs.getSecond();
176  }
177 
178  /**
179  * Advances current_ iterator
180  */
182  assert (current_ != std::prev(end(nondom_projections_)) && current_ != end(nondom_projections_));
183  ++current_;
184  }
185 
186  /**
187  * Incorporates a new projection and corresponding result into non-dominated projections
188  * @param proj Projection to incorporated
189  * @param res Corresponding result of projection
190  */
192  assert (current_ != std::prev(end(nondom_projections_)) && current_ != end(nondom_projections_));
193  auto it = add(proj, res);
194  if (epsilonDominates(proj, current_->first)) {
195  nondom_projections_.erase(current_);
196  current_ = it;
197  }
198 
199  while (std::next(it) != end(nondom_projections_) && epsilonDominates(proj, std::next(it)->first)) {
200  nondom_projections_.erase(std::next(it));
201  }
202  }
203 
204  /**
205  * Get outcomes corresponding to non-dominated projections
206  * @return Vector of outcomes corresponding to non-dominated projections
207  */
208  vector<OutcomeType> NondomProjections::getNondomProjOutcomes() const {
209  auto outcomes = vector<OutcomeType>{};
210  for (auto it=begin(nondom_projections_); it!=end(nondom_projections_); ++it) {
211  for (const auto& res : it->second)
212  outcomes.push_back(res.second);
213  }
214  return outcomes;
215  }
216 
217  /**
218  * Indicates that all stored projections are investigated
219  * @return true if all stored projections have been investigated; false otherwise
220  */
222  assert (current_ != end(nondom_projections_));
223  return current_ == std::prev(end(nondom_projections_));
224  }
225 
226  /**
227  * Copy constructor
228  * @param box RectangularBox to copy
229  */
230  RectangularBox::RectangularBox(const std::vector<Interval>& box)
231  : box_(begin(box), end(box))
232  {}
233 
234  /**
235  * Move constructor
236  * @param box RectangularBox to copy
237  */
238  RectangularBox::RectangularBox(std::vector<Interval>&& box)
239  : box_(begin(box), end(box))
240  {}
241 
242  /**
243  * Constructor: constructs box first_beg x ... x (first_end-1) x second x third_bex x ... x (third_end-1)
244  * @param first_beg Iterator referring to interval
245  * @param first_end Iterator referring to past-the-end interval
246  * @param second Middle interval
247  * @param third_beg Iterator referring to interval
248  * @param third_end Iterator referring to past-the-end interval
249  */
250  RectangularBox::RectangularBox(std::vector<Interval>::const_iterator first_beg,
251  std::vector<Interval>::const_iterator first_end,
252  Interval second,
253  std::vector<Interval>::const_iterator third_beg,
254  std::vector<Interval>::const_iterator third_end) {
255  std::copy(first_beg, first_end, std::back_inserter(box_));
256  box_.push_back(second);
257  std::copy(third_beg, third_end, std::back_inserter(box_));
258  }
259 
260  /**
261  * Get get number of intervals
262  * @return Dimension of rectangular box
263  */
264  size_t RectangularBox::size() const {
265  return box_.size();
266  }
267 
268  /**
269  * Get interval of box
270  * @param index Corresponding interval index
271  * @return Interval corresponding to index
272  */
274  assert (index < size());
275  return box_[index];
276  }
277 
278  /**
279  * Ostream operator
280  * @param os Output stream to write to
281  * @param box Box to write to stream
282  * @return Output stream
283  */
284  std::ostream &operator<<(std::ostream& os, const RectangularBox& box) {
285  for (auto interval : box.box_)
286  os << "[ " << interval.first << ", " << interval.second << " ) ";
287  os << "\n";
288  return os;
289  }
290 
291  /**
292  * Indicates whether given box is subset
293  * @param other Box to compare
294  * @return true if 'other' is subset; false otherwise
295  */
296  bool RectangularBox::isSupersetOf(const RectangularBox &other) const {
297  assert (other.box_.size() == this->box_.size());
298  for (size_t i=0; i<box_.size(); ++i) {
299  if (box_[i].first > other.box_[i].first || box_[i].second < other.box_[i].second)
300  return false;
301  }
302  return true;
303  }
304 
305  /**
306  * Indicates whether given box is superset
307  * @param other Box to compare
308  * @return true if 'other' is superset; false otherwise
309  */
310  bool RectangularBox::isSubsetOf(const RectangularBox &other) const {
311  assert (this->box_.size() == other.box_.size());
312  for (size_t i=0; i<box_.size(); ++i) {
313  if (box_[i].first < other.box_[i].first || box_[i].second > other.box_[i].second)
314  return false;
315  }
316  return true;
317  }
318 
319  /**
320  * Indicates whether given box is disjoint
321  * @param other Box to compare
322  * @return true if 'other' is disjoint; false otherwise
323  */
325  assert (this->box_.size() == other.box_.size());
326  for (size_t i=0; i<box_.size(); ++i) {
327  auto int_beg = std::max(box_[i].first, other.box_[i].first);
328  auto int_end = std::min(box_[i].second, other.box_[i].second);
329  if (int_beg > int_end)
330  return true;
331  }
332  return false;
333  }
334 
335  /**
336  * Indicates whether a_i + epsilon > e_i for all i
337  * @param epsilon Value to add to left interval limit
338  * @return true if a_i + epsilon <= e_i for all i; false otherwise
339  */
340  bool RectangularBox::isFeasible(double epsilon) const {
341  for (const auto& elem : box_) {
342  if (elem.first + epsilon > elem.second)
343  return false;
344  }
345  return true;
346  }
347 
348  /**
349  * Indicates whether outcome dominates entire box
350  * @param outcome Outcome to compare to
351  * @return true if given outcome dominates entire box; false otherwise
352  */
353  bool RectangularBox::isDominated(const OutcomeType& outcome) const {
354  assert (outcome.size() == box_.size());
355  for (size_t i=0; i<box_.size(); ++i) {
356  if (outcome[i] > box_[i].first) {
357  return false;
358  }
359  }
360  return true;
361  }
362 
363  /**
364  * Get interval intersection with respect to given dimension and given box
365  * @param index Interval index to take intersection
366  * @param other Box to consider intersection
367  * @return Interval intersection
368  */
369  RectangularBox::Interval RectangularBox::getIntervalIntersection(std::size_t index, const RectangularBox& other) const {
370  assert (box_.size() == other.box_.size());
371  auto int_beg = std::max(box_[index].first, other.box_[index].first);
372  auto int_end = std::min(box_[index].second, other.box_[index].second);
373  assert (int_beg <= int_end);
374  return {int_beg, int_end};
375  }
376 
377  /**
378  * Makes disjoint rectangular boxes with respect to given box
379  * @param delta Feasibility threshold
380  * @param other Box to compare to
381  * @return Container of disjoint boxes
382  */
383  vector<RectangularBox> RectangularBox::getDisjointPartsFrom(double delta, const RectangularBox &other) const {
384  auto size = this->box_.size();
385  assert (size == other.box_.size());
386  auto disjoint_partitions = vector<RectangularBox>{};
387  auto intersections = vector<Interval>{};
388  for (size_t i=0; i<size; ++i) {
389  if (box_[i].first < other.box_[i].first) { // non-empty to the left
390  auto new_box = RectangularBox(begin(intersections), end(intersections),
391  {box_[i].first, other.box_[i].first},
392  begin(box_)+(i+1), end(box_));
393  if (new_box.isFeasible(delta))
394  disjoint_partitions.push_back(std::move(new_box));
395 
396  }
397  if (other.box_[i].second < box_[i].second) { // non-empty to the right
398  auto new_box = RectangularBox(begin(intersections), end(intersections),
399  {other.box_[i].second, box_[i].second},
400  begin(box_)+(i+1), end(box_));
401  if (new_box.isFeasible(delta))
402  disjoint_partitions.push_back(std::move(new_box));
403  }
404  intersections.push_back(getIntervalIntersection(i, other));
405  }
406  return disjoint_partitions;
407  }
408 
409  /**
410  * Default constructor
411  * @param argc Argument count
412  * @param argv Argument vector
413  */
414  Polyscip::Polyscip(int argc, const char *const *argv)
415  : cmd_line_args_(argc, argv),
416  polyscip_status_(PolyscipStatus::Unsolved),
417  scip_(nullptr),
418  obj_sense_(SCIP_OBJSENSE_MINIMIZE), // default objective sense is minimization
419  no_objs_(0),
420  clock_total_(nullptr),
421  only_weight_space_phase_(false), // re-set in readProblem()
422  is_sub_prob_(false)
423  {
424  SCIPcreate(&scip_);
425  assert (scip_ != nullptr);
427  SCIPincludeObjReader(scip_, new ReaderMOP(scip_), TRUE);
428  SCIPcreateClock(scip_, addressof(clock_total_));
429 
430  if (cmd_line_args_.hasParameterFile()) {
431  if (filenameIsOkay(cmd_line_args_.getParameterFile())) {
432  SCIPreadParams(scip_, cmd_line_args_.getParameterFile().c_str());
433  }
434  else {
435  cout << "Invalid parameter settings file.\n";
436  polyscip_status_ = PolyscipStatus::Error;
437  }
438  }
439  if (cmd_line_args_.getDelta() <= 0) {
440  cout << "Please choose positive value for parameter --delta .\n";
441  polyscip_status_ = PolyscipStatus::Error;
442  }
443  if (cmd_line_args_.getEpsilon() <= 0) {
444  cout << "Please choose positive value for parameter --epsilon .\n";
445  polyscip_status_ = PolyscipStatus::Error;
446  }
447  if (cmd_line_args_.hasTimeLimit() && cmd_line_args_.getTimeLimit() <= 0) {
448  cout << "Please choose positive value for parameter --timelLimit .\n";
449  polyscip_status_ = PolyscipStatus::Error;
450  }
451  if (!filenameIsOkay(cmd_line_args_.getProblemFile())) {
452  cout << "Invalid problem file.\n";
453  polyscip_status_ = PolyscipStatus::Error;
454  }
455  }
456 
457  /**
458  * Constructor
459  * @param cmd_line_args Command line parameter object
460  * @param scip SCIP pointer
461  * @param no_objs Number of considered objective
462  * @param clock_total Clock measuring total computation time
463  */
464  Polyscip::Polyscip(const CmdLineArgs& cmd_line_args,
465  SCIP* scip,
466  std::size_t no_objs,
467  SCIP_CLOCK* clock_total)
468  : cmd_line_args_{cmd_line_args},
469  polyscip_status_{PolyscipStatus::ProblemRead},
470  scip_{scip},
471  obj_sense_{SCIP_OBJSENSE_MINIMIZE},//SCIPgetObjsense(scip_)
472  no_objs_{no_objs},
473  clock_total_{clock_total},
474  only_weight_space_phase_{false},
475  is_sub_prob_(true)
476  {}
477 
478  /**
479  * Destructor
480  */
482  if (!is_sub_prob_) {
483  SCIPfreeClock(scip_, addressof(clock_total_));
484  SCIPfree(addressof(scip_));
485  }
486  }
487 
488  /**
489  * Print PolySCIP status
490  * @param os Output stream to print to
491  */
492  void Polyscip::printStatus(std::ostream& os) const {
493  switch(polyscip_status_) {
495  os << "PolySCIP Status: ComputeProjectedNondomPointsPhase\n";
496  break;
498  os << "PolySCIP Status: Error\n";
499  break;
501  os << "PolySCIP Status: Successfully finished\n";
502  break;
504  os << "PolySCIP Status: LexOptPhase\n";
505  break;
507  os << "PolySCIP Status: ProblemRead\n";
508  break;
510  os << "PolySCIP Status: TimeLimitReached\n";
511  break;
513  os << "PolySCIP Status: Unsolved\n";
514  break;
516  os << "PolySCIP Status: WeightSpacePhase\n";
517  break;
518  }
519  }
520 
521  /**
522  * Compute non-dominated points of given problem
523  * @attention readProblem() needs to be called before
524  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
525  */
527  if (polyscip_status_ == PolyscipStatus::ProblemRead) {
528  SCIP_CALL(SCIPstartClock(scip_, clock_total_));
529 
530  auto obj_probdata = dynamic_cast<ProbDataObjectives*>(SCIPgetObjProbData(scip_));
531  auto nonzero_orig_vars = vector<vector<SCIP_VAR*>>{};
532  auto nonzero_orig_vals = vector<vector<ValueType>>{};
533  for (size_t obj=0; obj < no_objs_; ++obj) {
534  nonzero_orig_vars.push_back(obj_probdata->getNonZeroCoeffVars(obj));
535  assert (!nonzero_orig_vars.empty());
536  auto nonzero_obj_vals = vector<ValueType>{};
537  std::transform(nonzero_orig_vars.back().cbegin(),
538  nonzero_orig_vars.back().cend(),
539  std::back_inserter(nonzero_obj_vals),
540  [obj, obj_probdata](SCIP_VAR *var) { return obj_probdata->getObjCoeff(var, obj); });
541  nonzero_orig_vals.push_back(std::move(nonzero_obj_vals));
542  }
543 
544  SCIP_CALL( computeLexicographicOptResults(nonzero_orig_vars, nonzero_orig_vals) );
545 
546  if (polyscip_status_ == PolyscipStatus::LexOptPhase) {
547  if (no_objs_ > 3 || only_weight_space_phase_) {
548  SCIP_CALL(computeWeightSpaceResults());
549  }
550  else {
551  SCIP_CALL(computeNonLexicographicNondomResults(nonzero_orig_vars, nonzero_orig_vals));
552  }
553  }
554  SCIP_CALL(SCIPstopClock(scip_, clock_total_));
555  }
556  return SCIP_OKAY;
557  }
558 
559  /**
560  * Compute lexicographic optimal results
561  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
562  * @param orig_vals Container storing original non-zero objective coefficients for each objective
563  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
564  */
565  SCIP_RETCODE Polyscip::computeLexicographicOptResults(vector<vector<SCIP_VAR*>>& orig_vars,
566  vector<vector<ValueType>>& orig_vals) {
567  polyscip_status_ = PolyscipStatus::LexOptPhase;
568  for (size_t obj=0; obj<no_objs_; ++obj) {
569  if (polyscip_status_ == PolyscipStatus::LexOptPhase)
570  SCIP_CALL(computeLexicographicOptResult(obj, orig_vars, orig_vals));
571  else
572  break;
573  }
574  return SCIP_OKAY;
575  }
576 
577  /**
578  * Compute lexicographic optimal result with given objective having highest preference
579  * @param obj Objective with highest preference
580  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
581  * @param orig_vals Container storing original non-zero objective coefficients for each objective
582  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
583  */
584  SCIP_RETCODE Polyscip::computeLexicographicOptResult(size_t considered_obj,
585  vector<vector<SCIP_VAR*>>& orig_vars,
586  vector<vector<ValueType>>& orig_vals) {
587  assert (considered_obj < no_objs_);
588  auto current_obj = considered_obj;
589  auto obj_val_cons = vector<SCIP_CONS*>{};
590  auto weight = WeightType(no_objs_, 0.);
591  auto scip_status = SCIP_STATUS_UNKNOWN;
592  for (size_t counter = 0; counter<no_objs_; ++counter) {
593  weight[current_obj] = 1;
594  SCIP_CALL( setWeightedObjective(weight) );
595  SCIP_CALL( solve() );
596  scip_status = SCIPgetStatus(scip_);
597  if (scip_status == SCIP_STATUS_INFORUNBD)
598  scip_status = separateINFORUNBD(weight);
599 
600  if (scip_status == SCIP_STATUS_OPTIMAL) {
601  if (counter < no_objs_-1) {
602  auto opt_value = SCIPgetPrimalbound(scip_);
603 
605  auto cons = createObjValCons(orig_vars[current_obj],
606  orig_vals[current_obj],
607  opt_value,
608  opt_value);
609  SCIP_CALL (SCIPaddCons(scip_, cons));
610  obj_val_cons.push_back(cons);
611  }
612  }
613  else if (scip_status == SCIP_STATUS_UNBOUNDED) {
614  assert (current_obj == considered_obj);
615  SCIP_CALL( handleUnboundedStatus(true) );
616  unbd_orig_objs_.push_back(considered_obj);
617  break;
618  }
619  else if (scip_status == SCIP_STATUS_TIMELIMIT) {
620  polyscip_status_ = PolyscipStatus::TimeLimitReached;
621  break;
622  }
623  else if (scip_status == SCIP_STATUS_INFEASIBLE) {
624  assert (current_obj == 0);
625  polyscip_status_ = PolyscipStatus::Finished;
626  break;
627  }
628  else {
629  polyscip_status_ = PolyscipStatus::Error;
630  break;
631  }
632  weight[current_obj] = 0;
633  current_obj = (current_obj+1) % no_objs_;
634  }
635 
636  if (scip_status == SCIP_STATUS_OPTIMAL) {
637  auto lex_opt_result = getOptimalResult();
638  if (outcomeIsNew(lex_opt_result.second, begin(bounded_), end(bounded_))) {
639  bounded_.push_back(std::move(lex_opt_result));
640  }
641  }
642 
643  // release and delete added constraints
644  for (auto cons : obj_val_cons) {
645  SCIP_CALL( SCIPfreeTransform(scip_) );
646  SCIP_CALL( SCIPdelCons(scip_,cons) );
647  SCIP_CALL( SCIPreleaseCons(scip_,addressof(cons)) );
648  }
649  return SCIP_OKAY;
650  }
651 
652  /**
653  * Compute bounded non-dominated extreme points for objective for which unbounded ray exits
654  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
655  */
656  SCIP_RETCODE Polyscip::computeBoundedNondomResultsForUnbdObjs() {
657  auto dd = DDMethod(scip_, no_objs_, bounded_, unbounded_);
658  dd.computeVRep_Var1();
659  auto v_rep= dd.moveVRep();
660  for (auto unbd_obj : unbd_orig_objs_) {
661  for (const auto &v : v_rep) {
662  if (v->hasNonUnitAndNonZeroWeight()) {
663  auto weight = v->moveWeight();
664  assert (weight[unbd_obj] > 0.);
665  weight[unbd_obj] -= weight[unbd_obj] / 100;
666  SCIP_CALL(setWeightedObjective(weight));
667  SCIP_CALL(solve());
668  auto status = SCIPgetStatus(scip_);
669  if (status == SCIP_STATUS_OPTIMAL) {
670  auto bd_result = getOptimalResult();
671  if (outcomeIsNew(bd_result.second, begin(bounded_), end(bounded_))) {
672  bounded_.push_back(std::move(bd_result));
673  }
674  }
675  else if (status == SCIP_STATUS_TIMELIMIT) {
676  polyscip_status_ = PolyscipStatus::TimeLimitReached;
677  return SCIP_OKAY;
678  }
679  else {
680  polyscip_status_ = PolyscipStatus::Error;
681  return SCIP_OKAY;
682  }
683  }
684  }
685  }
686  return SCIP_OKAY;
687  }
688 
689  /**
690  * Compute non-dominated points which are not lexicographically optimal
691  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
692  * @param orig_vals Container storing original non-zero objective coefficients for each objective
693  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
694  */
695  SCIP_RETCODE Polyscip::computeNonLexicographicNondomResults(const vector<vector<SCIP_VAR *>> &orig_vars,
696  const vector<vector<ValueType>> &orig_vals) {
697  polyscip_status_ = PolyscipStatus::TwoProjPhase;
698  if (unboundedResultsExist()) {
699  SCIP_CALL( computeBoundedNondomResultsForUnbdObjs() );
700  }
701  // consider all (k over 2 ) combinations of considered objective functions
702  std::map<ObjPair, vector<OutcomeType>> proj_nondom_outcomes_map;
703  for (size_t obj_1=0; obj_1!=no_objs_-1; ++obj_1) {
704  for (auto obj_2=obj_1+1; obj_2!=no_objs_; ++obj_2) {
705  if (polyscip_status_ != PolyscipStatus::TwoProjPhase) {
706  return SCIP_OKAY;
707  }
708  else {
709  auto proj_nondom_outcomes = solveWeightedTchebycheff(orig_vars,
710  orig_vals,
711  obj_1, obj_2);
712 
713  proj_nondom_outcomes_map.insert({ObjPair(obj_1, obj_2), proj_nondom_outcomes});
714  }
715  }
716  }
717  if (no_objs_ == 3 && polyscip_status_ == PolyscipStatus::TwoProjPhase) {
718  auto feasible_boxes = computeFeasibleBoxes(proj_nondom_outcomes_map,
719  orig_vars,
720  orig_vals);
721  auto disjoint_boxes = computeDisjointBoxes(std::move(feasible_boxes));
722  assert (feasible_boxes.size() <= disjoint_boxes.size());
723 
724  for (const auto& box : disjoint_boxes) {
725  if (std::any_of(begin(bounded_),
726  end(bounded_),
727  [&box](const Result& res){return box.isDominated(res.second);})) {
728  continue;
729  }
730  if (cmd_line_args_.beVerbose()) {
731  cout << "Checking box = " << box << "\n";
732  }
733  auto new_res = computeNondomPointsInBox(box,
734  orig_vars,
735  orig_vals);
736  for (const auto &res : new_res) {
737  if (polyscip_status_ != PolyscipStatus::TwoProjPhase) {
738  return SCIP_OKAY;
739  }
740  if (is_sub_prob_) {
741  bounded_.push_back(res);
742  }
743  else {
744  if (!boxResultIsDominated(res.second, orig_vars, orig_vals)) {
745  bounded_.push_back(std::move(res));
746  }
747  }
748  }
749  }
750  }
751  if (polyscip_status_ == PolyscipStatus::TwoProjPhase)
752  polyscip_status_ = PolyscipStatus::Finished;
753 
754  return SCIP_OKAY;
755  }
756 
757  /**
758  * Indicates whether given outcome is globally dominated
759  * @param outcome Outcome to check for dominance
760  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
761  * @param orig_vals Container storing original non-zero objective coefficients for each objective
762  * @return true if given outcome is dominated; false otherwise
763  */
764  bool Polyscip::boxResultIsDominated(const OutcomeType& outcome,
765  const vector<vector<SCIP_VAR*>>& orig_vars,
766  const vector<vector<ValueType>>& orig_vals) {
767 
768  auto size = outcome.size();
769  assert (size == orig_vars.size());
770  assert (size == orig_vals.size());
771  auto is_dominated = false;
772 
774 
775  auto new_cons = vector<SCIP_CONS*>{};
776  // add objective value constraints
777  for (size_t i=0; i<size; ++i) {
778  auto cons = createObjValCons(orig_vars[i],
779  orig_vals[i],
780  -SCIPinfinity(scip_),
781  outcome[i]);
782  SCIP_CALL_ABORT( SCIPaddCons(scip_, cons) );
783  new_cons.push_back(cons);
784  }
785 
786  auto weight = WeightType(no_objs_, 1.);
787  setWeightedObjective(weight/*zero_weight*/);
788  solve(); // compute with zero objective
789  auto scip_status = SCIPgetStatus(scip_);
790 
791  // check solution status
792  if (scip_status == SCIP_STATUS_OPTIMAL) {
793  assert (weight.size() == outcome.size());
794  if (SCIPgetPrimalbound(scip_) + cmd_line_args_.getEpsilon() < std::inner_product(begin(weight),
795  end(weight),
796  begin(outcome),
797  0.)) {
798  is_dominated = true;
799  }
800  }
801  else if (scip_status == SCIP_STATUS_TIMELIMIT) {
802  polyscip_status_ = PolyscipStatus::TimeLimitReached;
803  }
804  else {
805  polyscip_status_ = PolyscipStatus::Error;
806  }
807 
809  // release and delete added constraints
810  for (auto cons : new_cons) {
811  SCIP_CALL_ABORT( SCIPdelCons(scip_, cons) );
812  SCIP_CALL_ABORT( SCIPreleaseCons(scip_, addressof(cons)) );
813  }
814 
815  return is_dominated;
816  }
817 
818  /**
819  * Compute locally non-dominated results in given rectangular box
820  * @param box Rectangular box
821  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
822  * @param orig_vals Container storing original non-zero objective coefficients for each objective
823  * @return Container with locally non-dominated results
824  */
825  ResultContainer Polyscip::computeNondomPointsInBox(const RectangularBox& box,
826  const vector<vector<SCIP_VAR *>>& orig_vars,
827  const vector<vector<ValueType>>& orig_vals) {
828  assert (box.size() == orig_vars.size());
829  assert (box.size() == orig_vals.size());
830  // add constraints on objective values given by box
831  auto obj_val_cons = vector<SCIP_CONS *>{};
832  if (SCIPisTransformed(scip_)) {
834  }
835  for (size_t i=0; i<box.size(); ++i) {
836  auto interval = box.getInterval(i);
837  auto new_cons = createObjValCons(orig_vars[i],
838  orig_vals[i],
839  interval.first,
840  interval.second - cmd_line_args_.getDelta());
841  SCIP_CALL_ABORT( SCIPaddCons(scip_, new_cons) );
842  obj_val_cons.push_back(new_cons);
843  }
844 
845  std::unique_ptr<Polyscip> sub_poly(new Polyscip(cmd_line_args_,
846  scip_,
847  no_objs_,
848  clock_total_) );
849  sub_poly->computeNondomPoints();
850 
851  // release and delete objective value constraints
852  if (SCIPisTransformed(scip_)) {
854  }
855  for (auto cons : obj_val_cons) {
856  SCIP_CALL_ABORT( SCIPdelCons(scip_, cons) );
857  SCIP_CALL_ABORT( SCIPreleaseCons(scip_, addressof(cons)) );
858  }
859 
860  // check computed subproblem results
861  assert (!sub_poly->unboundedResultsExist());
862 
863  auto new_nondom_res = ResultContainer{};
864  if (sub_poly->getStatus() == PolyscipStatus::Finished) {
865  if (sub_poly->numberOfBoundedResults() > 0) {
866  for (auto it = sub_poly->boundedCBegin(); it != sub_poly->boundedCEnd(); ++it) {
867  new_nondom_res.push_back(std::move(*it));
868  }
869  }
870  }
871  else if (sub_poly->getStatus() == PolyscipStatus::TimeLimitReached) {
872  polyscip_status_ = PolyscipStatus::TimeLimitReached;
873  }
874  else {
875  polyscip_status_ = PolyscipStatus::Error;
876  }
877  sub_poly.reset();
878  return new_nondom_res;
879  }
880 
881  /**
882  * Compute disjoint rectangular boxes from given feasible rectangular boxes
883  * @param feasible_boxes List of feasible boxes
884  * @return Vector of disjoint feasible rectangular boxes
885  */
886  vector<RectangularBox> Polyscip::computeDisjointBoxes(list<RectangularBox>&& feasible_boxes) const {
887  // delete redundant boxes
888  auto current = begin(feasible_boxes);
889  while (current != end(feasible_boxes)) {
890  auto increment_current = true;
891  auto it = begin(feasible_boxes);
892  while (it != end(feasible_boxes)) {
893  if (current != it) {
894  if (current->isSupersetOf(*it)) {
895  it = feasible_boxes.erase(it);
896  continue;
897  }
898  else if (current->isSubsetOf(*it)) {
899  current = feasible_boxes.erase(current);
900  increment_current = false;
901  break;
902  }
903  }
904  ++it;
905  }
906  if (increment_current)
907  ++current;
908  }
909  // compute disjoint boxes
910  auto disjoint_boxes = vector<RectangularBox>{};
911  while (!feasible_boxes.empty()) {
912  auto box_to_be_added = feasible_boxes.back();
913  feasible_boxes.pop_back();
914 
915  auto current_boxes = vector<RectangularBox>{};
916  for (const auto& elem : disjoint_boxes) {
917  assert (!box_to_be_added.isSubsetOf(elem));
918  if (box_to_be_added.isDisjointFrom(elem)) {
919  current_boxes.push_back(elem);
920  }
921  else if (box_to_be_added.isSupersetOf(elem)) {
922  continue;
923  }
924  else {
925  auto elem_disjoint = elem.getDisjointPartsFrom(cmd_line_args_.getDelta(), box_to_be_added);
926  std::move(begin(elem_disjoint), end(elem_disjoint), std::back_inserter(current_boxes));
927  }
928  }
929  disjoint_boxes.clear();
930  std::move(begin(current_boxes), end(current_boxes), std::back_inserter(disjoint_boxes));
931  disjoint_boxes.push_back(std::move(box_to_be_added));
932  }
933  return disjoint_boxes;
934  }
935 
936  /**
937  * Compute feasible rectangular boxes
938  * @param proj_nondom_outcomes Non-dominated outcomes which are non-dominated for objective pair
939  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
940  * @param orig_vals Container storing original non-zero objective coefficients for each objective
941  * @return List of feasible rectangular boxes
942  */
943  list<RectangularBox> Polyscip::computeFeasibleBoxes(const map<ObjPair, vector<OutcomeType>> &proj_nd_outcomes,
944  const vector<vector<SCIP_VAR *>> &orig_vars,
945  const vector<vector<ValueType>> &orig_vals) {
946 
947  auto& nd_outcomes_01 = proj_nd_outcomes.at(ObjPair(0,1));
948  assert (!nd_outcomes_01.empty());
949  auto& nd_outcomes_02 = proj_nd_outcomes.at(ObjPair(0,2));
950  assert (!nd_outcomes_02.empty());
951  auto& nd_outcomes_12 = proj_nd_outcomes.at(ObjPair(1,2));
952  assert (!nd_outcomes_12.empty());
953 
954  auto feasible_boxes = list<RectangularBox>{};
955  for (const auto& nd_01 : nd_outcomes_01) {
956  for (const auto& nd_02 : nd_outcomes_02) {
957  for (const auto& nd_12 : nd_outcomes_12) {
958  auto box = RectangularBox({{max(nd_01[0], nd_02[0]), nd_12[0]},
959  {max(nd_01[1], nd_12[1]), nd_02[1]},
960  {max(nd_02[2], nd_12[2]), nd_01[2]}});
961  if (box.isFeasible(cmd_line_args_.getDelta())) {
962  feasible_boxes.push_back(box);
963  }
964  }
965  }
966  }
967  return feasible_boxes;
968  }
969 
970 
971  /**
972  * Create constraint: new_var - beta_i * orig_vals \\cdot orig_vars >= - beta_i * rhs
973  * @param new_var Non-original variable
974  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
975  * @param orig_vals Container storing original non-zero objective coefficients for each objective
976  * @param rhs rhs value
977  * @param beta_i coefficient
978  * @return Pointer to corresponding SCIP constraint
979  */
980  SCIP_CONS* Polyscip::createNewVarTransformCons(SCIP_VAR *new_var,
981  const vector<SCIP_VAR *> &orig_vars,
982  const vector<ValueType> &orig_vals,
983  const ValueType &rhs,
984  const ValueType &beta_i) {
985  auto vars = vector<SCIP_VAR*>(begin(orig_vars), end(orig_vars));
986  auto vals = vector<ValueType>(orig_vals.size(), 0.);
987  std::transform(begin(orig_vals),
988  end(orig_vals),
989  begin(vals),
990  [beta_i](ValueType val){return -beta_i*val;});
991  vars.push_back(new_var);
992  vals.push_back(1.);
993 
994  SCIP_CONS* cons = nullptr;
995  // add contraint new_var - beta_i* vals \cdot vars >= - beta_i * rhs
997  addressof(cons),
998  "new_variable_transformation_constraint",
999  global::narrow_cast<int>(vars.size()),
1000  vars.data(),
1001  vals.data(),
1002  -beta_i*rhs,
1003  SCIPinfinity(scip_));
1004  assert (cons != nullptr);
1005  return cons;
1006  }
1007 
1008  /**
1009  * Create constraint: lhs <= vals \\cdot vars <= rhs
1010  * @param vars Considered variables
1011  * @param vals Considered coefficient values
1012  * @param lhs lhs value
1013  * @param rhs rhs value
1014  * @return Pointer to corresponding SCIP constraint
1015  */
1016  SCIP_CONS* Polyscip::createObjValCons(const vector<SCIP_VAR *>& vars,
1017  const vector<ValueType>& vals,
1018  const ValueType& lhs,
1019  const ValueType& rhs) {
1020  SCIP_CONS* cons = nullptr;
1021  std::remove_const<const vector<SCIP_VAR*>>::type non_const_vars(vars);
1022  std::remove_const<const vector<ValueType>>::type non_const_vals(vals);
1024  addressof(cons),
1025  "lhs <= c_i^T x <= rhs",
1026  global::narrow_cast<int>(vars.size()),
1027  non_const_vars.data(),
1028  non_const_vals.data(),
1029  lhs,
1030  rhs);
1031  assert (cons != nullptr);
1032  return cons;
1033  }
1034 
1035  /**
1036  * Computes non-dominated point which fulfills: obj_val_cons1 = obj_val_cons1_rhs and obj_val_cons2 = obj_val_cons2_rhs
1037  * @param obj_val_cons1 First constraint to consider
1038  * @param obj_val_cons2 Second constraint to consider
1039  * @param obj_val_cons1_rhs Corresponding rhs of first constraint
1040  * @param obj_val_cons2_rhs Corresponding rhs of second constraint
1041  * @param obj_1 Considered objective index corresponding to first constraint
1042  * @param obj_2 Considered objective index corresponding to second constraint
1043  * @param results Container to store computed non-dominated result
1044  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1045  */
1046  SCIP_RETCODE Polyscip::computeNondomProjResult(SCIP_CONS *cons1,
1047  SCIP_CONS *cons2,
1048  ValueType rhs_cons1,
1049  ValueType rhs_cons2,
1050  size_t obj_1,
1051  size_t obj_2,
1052  ResultContainer &results) {
1053 
1054  SCIP_CALL( SCIPchgRhsLinear(scip_, cons1, rhs_cons1) );
1055  SCIP_CALL( SCIPchgRhsLinear(scip_, cons2, rhs_cons2) );
1056  // set new objective function
1057  auto intermed_obj = WeightType(no_objs_, 0.);
1058  intermed_obj.at(obj_1) = 1.;
1059  intermed_obj.at(obj_2) = 1.;
1060  SCIP_CALL( setWeightedObjective(intermed_obj) );
1061 
1062  // solve auxiliary problem
1063  SCIP_CALL( solve() );
1064  auto scip_status = SCIPgetStatus(scip_);
1065  if (scip_status == SCIP_STATUS_OPTIMAL) {
1066  if (no_objs_ > 2) {
1067  auto intermed_result = getOptimalResult();
1068  SCIP_CALL(SCIPchgLhsLinear(scip_, cons1, intermed_result.second[obj_1]));
1069  SCIP_CALL(SCIPchgRhsLinear(scip_, cons1, intermed_result.second[obj_1]));
1070  SCIP_CALL(SCIPchgLhsLinear(scip_, cons2, intermed_result.second[obj_2]));
1071  SCIP_CALL(SCIPchgRhsLinear(scip_, cons2, intermed_result.second[obj_2]));
1072  SCIP_CALL(setWeightedObjective(WeightType(no_objs_, 1.)));
1073  SCIP_CALL(solve());
1074  if (SCIPgetStatus(scip_) == SCIP_STATUS_TIMELIMIT) {
1075  polyscip_status_ = PolyscipStatus::TimeLimitReached;
1076  }
1077  else {
1078  assert (SCIPgetStatus(scip_) == SCIP_STATUS_OPTIMAL);
1079  }
1080  }
1081 
1082  auto nondom_result = getOptimalResult();
1083  results.push_back(std::move(nondom_result));
1084  }
1085  else if (scip_status == SCIP_STATUS_TIMELIMIT) {
1086  polyscip_status_ = PolyscipStatus::TimeLimitReached;
1087  }
1088  else {
1089  cout << "unexpected SCIP status in computeNondomProjResult: " +
1090  std::to_string(int(SCIPgetStatus(scip_))) + "\n";
1091  polyscip_status_ = PolyscipStatus::Error;
1092  }
1093 
1094  // unset objective function
1095  SCIP_CALL( setWeightedObjective( WeightType(no_objs_, 0.)) );
1096  return SCIP_OKAY;
1097  }
1098 
1099 
1100  /**
1101  * Compute non-dominated points via subproblems with weighted Tchebycheff norm
1102  * @param orig_vars Container storing original problem variables with non-zero coefficients for each objective
1103  * @param orig_vals Container storing original non-zero objective coefficients for each objective
1104  * @param obj_1 Index of first considered objective
1105  * @param obj_2 Index of second considered objective
1106  * @return Container of non-dominated outcomes which are also non-dominated for projection onto obj_1 and obj_2
1107  */
1108  vector<OutcomeType> Polyscip::solveWeightedTchebycheff(const vector<vector<SCIP_VAR*>>& orig_vars,
1109  const vector<vector<ValueType>>& orig_vals,
1110  size_t obj_1,
1111  size_t obj_2){
1112 
1113  assert (orig_vars.size() == orig_vals.size());
1114  assert (orig_vals.size() == no_objs_);
1115  assert (obj_1 < no_objs_ && obj_2 < no_objs_);
1116 
1117  // change objective values of existing variabless to zero
1118  setWeightedObjective(WeightType(no_objs_, 0.));
1119 
1120  // add new variable with objective value = 1 (for transformed Tchebycheff norm objective)
1121  SCIP_VAR* z = nullptr;
1122  SCIPcreateVarBasic(scip_,
1123  addressof(z),
1124  "z",
1125  -SCIPinfinity(scip_),
1126  SCIPinfinity(scip_),
1127  1,
1129  assert (z != nullptr);
1130  SCIPaddVar(scip_, z);
1131 
1132  auto nondom_projs = NondomProjections(cmd_line_args_.getEpsilon(),
1133  bounded_,
1134  obj_1,
1135  obj_2);
1136 
1137  auto last_proj = nondom_projs.getLastProj();
1138  while (!nondom_projs.finished() && polyscip_status_ == PolyscipStatus::TwoProjPhase) {
1139  auto left_proj = nondom_projs.getLeftProj();
1140  auto right_proj = nondom_projs.getRightProj();
1141 
1142  assert (left_proj.getFirst() < right_proj.getFirst());
1143  assert (left_proj.getSecond() > last_proj.getSecond());
1144 
1145  // create constraint pred.first <= c_{objs.first} \cdot x <= succ.first
1146  auto obj_val_cons = vector<SCIP_CONS*>{};
1147  obj_val_cons.push_back( createObjValCons(orig_vars[obj_1],
1148  orig_vals[obj_1],
1149  left_proj.getFirst(),
1150  right_proj.getFirst()));
1151  // create constraint optimal_val_objs.second <= c_{objs.second} \cdot x <= pred.second
1152  obj_val_cons.push_back( createObjValCons(orig_vars[obj_2],
1153  orig_vals[obj_2],
1154  last_proj.getSecond(),
1155  left_proj.getSecond()));
1156  for (auto c : obj_val_cons) {
1157  SCIP_CALL_ABORT( SCIPaddCons(scip_, c) );
1158  }
1159 
1160  auto ref_point = std::make_pair(left_proj.getFirst() - 1., last_proj.getSecond() - 1.);
1161  // set beta = (beta_1,beta_2) s.t. pred and succ are both on the norm rectangle defined by beta
1162  auto beta_1 = 1.0;
1163  auto beta_2 = (right_proj.getFirst() - ref_point.first) / (left_proj.getSecond() - ref_point.second);
1164  // create constraint with respect to beta_1
1165  auto var_trans_cons = vector<SCIP_CONS*>{};
1166  var_trans_cons.push_back( createNewVarTransformCons(z,
1167  orig_vars[obj_1],
1168  orig_vals[obj_1],
1169  ref_point.first,
1170  beta_1));
1171  // create constraint with respect to beta_2
1172  var_trans_cons.push_back(createNewVarTransformCons(z,
1173  orig_vars[obj_2],
1174  orig_vals[obj_2],
1175  ref_point.second,
1176  beta_2));
1177  for (auto c : var_trans_cons) {
1178  SCIP_CALL_ABORT( SCIPaddCons(scip_, c) );
1179  }
1180 
1181  solve();
1182  std::unique_ptr<TwoDProj> new_proj(nullptr);
1183  auto scip_status = SCIPgetStatus(scip_);
1184  if (scip_status == SCIP_STATUS_OPTIMAL) {
1185  assert (SCIPisGE(scip_, SCIPgetPrimalbound(scip_), 0.));
1186  auto res = getOptimalResult();
1187  new_proj = global::make_unique<TwoDProj>(res.second, obj_1, obj_2);
1188  assert (new_proj);
1189  }
1190  else if (scip_status == SCIP_STATUS_TIMELIMIT) {
1191  polyscip_status_ = PolyscipStatus::TimeLimitReached;
1192  }
1193  else if (scip_status == SCIP_STATUS_INFEASIBLE) {
1194  std::cerr << "Numerical troubles between " << left_proj << " and " << right_proj << "\n";
1195  std::cerr << "Continuing with next subproblem.\n";
1196  nondom_projs.update();
1197  }
1198  else {
1199  std::cerr << "Unexpected SCIP status in solveWeightedTchebycheff: " +
1200  std::to_string(int(SCIPgetStatus(scip_))) + "\n";
1201  polyscip_status_ = PolyscipStatus::Error;
1202  }
1203 
1204  // release and delete variable transformation constraints
1206  for (auto c : var_trans_cons) {
1207  SCIP_CALL_ABORT( SCIPdelCons(scip_, c) );
1208  SCIP_CALL_ABORT( SCIPreleaseCons(scip_, addressof(c)) );
1209  }
1210 
1211  if (new_proj) {
1212  if (nondom_projs.epsilonDominates(left_proj, *new_proj) ||
1213  nondom_projs.epsilonDominates(right_proj, *new_proj)) {
1214  nondom_projs.update();
1215  }
1216  else {
1217  // temporarily delete variable 'z' from problem
1218  SCIP_Bool var_deleted = FALSE;
1219  SCIPdelVar(scip_, z, addressof(var_deleted));
1220  assert (var_deleted);
1221 
1222  computeNondomProjResult(obj_val_cons.front(), // constraint wrt obj_1
1223  obj_val_cons.back(), // constraint wrt obj2
1224  new_proj->getFirst(),
1225  new_proj->getSecond(),
1226  obj_1,
1227  obj_2,
1228  bounded_);
1229  auto nd_proj = TwoDProj(bounded_.back().second, obj_1, obj_2);
1230 
1231  nondom_projs.update(std::move(nd_proj), bounded_.back());
1232 
1233  // add variable 'z' back to problem
1234  SCIPaddVar(scip_, z);
1235 
1236  }
1237  new_proj.reset();
1238  }
1239 
1241  for (auto c : obj_val_cons) {
1242  SCIP_CALL_ABORT( SCIPdelCons(scip_, c) );
1243  SCIP_CALL_ABORT( SCIPreleaseCons(scip_, addressof(c)) );
1244  }
1245 
1246  }
1247 
1248  // clean up
1249  SCIP_Bool var_deleted = FALSE;
1250  SCIPdelVar(scip_, z, addressof(var_deleted));
1251  assert (var_deleted);
1252  SCIPreleaseVar(scip_, addressof(z));
1253 
1254  return nondom_projs.getNondomProjOutcomes();
1255 
1256  }
1257 
1258  /**
1259  * Get PolySCIP status
1260  * @return Current PolySCIP status
1261  */
1263  return polyscip_status_;
1264  }
1265 
1266  /**
1267  * Get number of bounded results
1268  * @return Number of computed bounded results
1269  */
1270  std::size_t Polyscip::numberOfBoundedResults() const {
1271  return bounded_.size();
1272  }
1273 
1274  /**
1275  * Get number of unbounded results
1276  * @return Number of computed unbounded results
1277  */
1279  return unbounded_.size();
1280  }
1281 
1282  /**
1283  * Resolve INFORUNBD SCIP status to either infeasible or unbounded
1284  * @param weight Weight yielding INFORUNBD status
1285  * @param with_presolving Indicates whether presolving should be used or not
1286  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1287  */
1288  SCIP_STATUS Polyscip::separateINFORUNBD(const WeightType& weight, bool with_presolving) {
1289  if (!with_presolving)
1291  auto zero_weight = WeightType(no_objs_, 0.);
1292  setWeightedObjective(zero_weight);
1293  solve(); // re-compute with zero objective
1294  if (!with_presolving)
1296  auto status = SCIPgetStatus(scip_);
1297  setWeightedObjective(weight); // re-set to previous objective
1298  if (status == SCIP_STATUS_INFORUNBD) {
1299  if (with_presolving)
1300  separateINFORUNBD(weight, false);
1301  else {
1302  cout << "INFORUNBD Status for problem with zero objective and no presolving.\n";
1303  polyscip_status_ = PolyscipStatus::Error;
1304  }
1305  }
1306  else if (status == SCIP_STATUS_UNBOUNDED) {
1307  cout << "UNBOUNDED Status for problem with zero objective.\n";
1308  polyscip_status_ = PolyscipStatus::Error;
1309  }
1310  else if (status == SCIP_STATUS_OPTIMAL) { // previous problem was unbounded
1311  return SCIP_STATUS_UNBOUNDED;
1312  }
1313  return status;
1314  }
1315 
1316  /**
1317  * Handle SCIP status that is neither optimal nor unbounded
1318  * @param status Current SCIP status
1319  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1320  */
1321  SCIP_RETCODE Polyscip::handleNonOptNonUnbdStatus(SCIP_STATUS status) {
1322  assert (status != SCIP_STATUS_OPTIMAL && status != SCIP_STATUS_UNBOUNDED);
1323  if (status == SCIP_STATUS_TIMELIMIT) {
1324  polyscip_status_ = PolyscipStatus::TimeLimitReached;
1325  }
1326  else if (is_sub_prob_) {
1327  assert (status == SCIP_STATUS_INFORUNBD || status == SCIP_STATUS_INFEASIBLE);
1328  polyscip_status_ = PolyscipStatus::Finished;
1329  }
1330  else {
1331  polyscip_status_ = PolyscipStatus::Error;
1332  }
1333  return SCIP_OKAY;
1334  }
1335 
1336  /**
1337  * Handle unbounded SCIP status
1338  * @param check_if_new_result Indicates whether to check if computed results is already known
1339  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1340  */
1341  SCIP_RETCODE Polyscip::handleUnboundedStatus(bool check_if_new_result) {
1342  if (!SCIPhasPrimalRay(scip_)) {
1344  if (SCIPisTransformed(scip_))
1345  SCIP_CALL( SCIPfreeTransform(scip_) );
1346  SCIP_CALL( solve() );
1348  if (SCIPgetStatus(scip_) != SCIP_STATUS_UNBOUNDED || !SCIPhasPrimalRay(scip_)) {
1349  polyscip_status_ = PolyscipStatus::Error;
1350  return SCIP_OKAY;
1351  }
1352  }
1353  auto result = getResult(false);
1354  if (!check_if_new_result || outcomeIsNew(result.second, false)) {
1355  unbounded_.push_back(std::move(result));
1356  }
1357  else if (cmd_line_args_.beVerbose()) {
1358  global::print(result.second, "Outcome: [", "]");
1359  cout << "not added to results.\n";
1360  }
1361  return SCIP_OKAY;
1362  }
1363 
1364  /**
1365  * Get computed result
1366  * @param outcome_is_bounded Indicates whether previous computation yielded unbounded status
1367  * @param primal_sol Corresponding SCIP primal solution pointer if previous computation yielded optimal status
1368  * @return Result type
1369  */
1370  Result Polyscip::getResult(bool outcome_is_bounded, SCIP_SOL *primal_sol) {
1371  SolType sol;
1372  auto outcome = OutcomeType(no_objs_,0.);
1373  auto no_vars = SCIPgetNOrigVars(scip_);
1374  auto vars = SCIPgetOrigVars(scip_);
1375  auto objs_probdata = dynamic_cast<ProbDataObjectives*>(SCIPgetObjProbData(scip_));
1376  for (auto i=0; i<no_vars; ++i) {
1377  auto var_sol_val = outcome_is_bounded ? SCIPgetSolVal(scip_, primal_sol, vars[i]) :
1378  SCIPgetPrimalRayVal(scip_, vars[i]);
1379 
1380  if (!SCIPisZero(scip_, var_sol_val)) {
1381  sol.emplace_back(SCIPvarGetName(vars[i]), var_sol_val);
1382  auto var_obj_vals = OutcomeType(no_objs_, 0.);
1383  for (size_t index=0; index!=no_objs_; ++index) {
1384  var_obj_vals[index] = objs_probdata->getObjVal(vars[i], index, var_sol_val);
1385  }
1386  std::transform(begin(outcome), end(outcome),
1387  begin(var_obj_vals),
1388  begin(outcome),
1389  std::plus<ValueType>());
1390 
1391  }
1392  }
1393  return {sol, outcome};
1394  }
1395 
1396  /**
1397  * Get bounded optimal result
1398  * @return Result type
1399  */
1400  Result Polyscip::getOptimalResult() {
1401  auto best_sol = SCIPgetBestSol(scip_);
1402  assert (best_sol != nullptr);
1403  SCIP_SOL *finite_sol{nullptr};
1404  SCIP_Bool same_obj_val{FALSE};
1405  auto retcode = SCIPcreateFiniteSolCopy(scip_, addressof(finite_sol), best_sol, addressof(same_obj_val));
1406  if (retcode != SCIP_OKAY)
1407  throw std::runtime_error("SCIPcreateFiniteSolCopy: return code != SCIP_OKAY.\n");
1408  if (!same_obj_val) {
1409  auto diff = std::fabs(SCIPgetSolOrigObj(scip_, best_sol) -
1410  SCIPgetSolOrigObj(scip_, finite_sol));
1411  if (diff > 1.0e-5) {
1412  std::cerr << "absolute value difference after calling SCIPcreateFiniteSolCopy: " << diff << "\n";
1413  SCIPfreeSol(scip_, addressof(finite_sol));
1414  throw std::runtime_error("SCIPcreateFiniteSolCopy: unacceptable difference in objective values.");
1415  }
1416  }
1417  assert (finite_sol != nullptr);
1418  auto new_result = getResult(true, finite_sol);
1419  SCIPfreeSol(scip_, addressof(finite_sol));
1420  return new_result;
1421  }
1422 
1423  /**
1424  * Indicates whether given outcome was not computed before
1425  * @param outcome Outcome to check
1426  * @param outcome_is_bounded Indicates whether given outcome is bounded or unbounded
1427  * @return true if given outcome was not computed before; false otherwise
1428  */
1429  bool Polyscip::outcomeIsNew(const OutcomeType& outcome, bool outcome_is_bounded) const {
1430  auto beg_it = outcome_is_bounded ? begin(bounded_) : begin(unbounded_);
1431  auto end_it = outcome_is_bounded ? end(bounded_) : end(unbounded_);
1432  return std::find_if(beg_it, end_it, [&outcome](const Result& res){return outcome == res.second;}) == end_it;
1433  }
1434 
1435 
1436  /**
1437  * Indicates whether given outcome is new with respect to other given results
1438  * @param outcome Outcome to check
1439  * @param beg Const_iterator to beginning of result container
1440  * @param last Const_iterator to past-the-end of result container
1441  * @return true if given outcome does not coincide with outcomes; false otherwise
1442  */
1443  bool Polyscip::outcomeIsNew(const OutcomeType& outcome,
1444  ResultContainer::const_iterator beg,
1445  ResultContainer::const_iterator last) const {
1446  auto eps = cmd_line_args_.getEpsilon();
1447  return std::none_of(beg, last, [eps, &outcome](const Result& res)
1448  {
1449  return outcomesCoincide(outcome, res.second, eps);
1450  });
1451  }
1452 
1453  /**
1454  * Indicates whether given outcomes coincide within some epsilon error
1455  * @param a First outcome to compare
1456  * @param b Second outcome to compare
1457  * @param epsilon Allowed error
1458  * @return true if outcomes coincides; false otherwise
1459  */
1460  bool Polyscip::outcomesCoincide(const OutcomeType& a, const OutcomeType& b, double epsilon) {
1461  assert (a.size() == b.size());
1462  return std::equal(begin(a), end(a), begin(b),
1463  [epsilon](ValueType v, ValueType w)
1464  {
1465  return std::fabs(v-w) < epsilon;
1466  });
1467  }
1468 
1469  /**
1470  * Solves currently considered SCIP instance
1471  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1472  */
1473  SCIP_RETCODE Polyscip::solve() {
1474  if (cmd_line_args_.hasTimeLimit()) { // set SCIP timelimit
1475  auto remaining_time = std::max(cmd_line_args_.getTimeLimit() -
1476  SCIPgetClockTime(scip_, clock_total_), 0.);
1477  SCIP_CALL(SCIPsetRealParam(scip_, "limits/time", remaining_time));
1478  }
1479  SCIP_CALL( SCIPsolve(scip_) ); // actual SCIP solver call
1480  return SCIP_OKAY;
1481  }
1482 
1483  /**
1484  * Set weighted objective: weight * (c_1,...,c_k) \\cdot x
1485  * @param weight Weight
1486  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1487  */
1488  SCIP_RETCODE Polyscip::setWeightedObjective(const WeightType& weight){
1489  if (SCIPisTransformed(scip_))
1490  SCIP_CALL( SCIPfreeTransform(scip_) );
1491  auto obj_probdata = dynamic_cast<ProbDataObjectives*>(SCIPgetObjProbData(scip_));
1492  assert (obj_probdata != nullptr);
1493  auto vars = SCIPgetOrigVars(scip_);
1494  auto no_vars = SCIPgetNOrigVars(scip_);
1495  if (weight == WeightType(weight.size(), 0.)) { // weight coincides with all zeros vector
1496  for (auto i = 0; i < no_vars; ++i) {
1497  SCIP_CALL(SCIPchgVarObj(scip_, vars[i], 0.));
1498  }
1499  }
1500  else { // weight != all zeros vector
1501  for (auto i = 0; i < no_vars; ++i) {
1502  auto val = obj_probdata->getWeightedObjVal(vars[i], weight);
1503  auto no_decimals = cmd_line_args_.roundWeightedObjCoeff();
1504  if (no_decimals) {
1505  ValueType intpart, fractpart;
1506  fractpart = std::modf(val, &intpart);
1507  fractpart = round(pow(10,no_decimals)*fractpart)/pow(10,no_decimals);
1508  val = intpart + fractpart;
1509  }
1510  SCIP_CALL(SCIPchgVarObj(scip_, vars[i], val));
1511  }
1512  }
1513  return SCIP_OKAY;
1514  }
1515 
1516  /**
1517  * Compute non-dominated extreme point results
1518  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1519  */
1520  SCIP_RETCODE Polyscip::computeWeightSpaceResults() {
1521  polyscip_status_ = PolyscipStatus::WeightSpacePhase;
1522  auto v_rep = DDMethod(scip_, no_objs_, bounded_, unbounded_);
1523  v_rep.computeVRep_Var1();
1524  weight_space_poly_ = global::make_unique<WeightSpacePolyhedron>(no_objs_,
1525  v_rep.moveVRep(),
1526  v_rep.moveHRep());
1527 
1528  while (weight_space_poly_->hasUntestedWeight() && polyscip_status_ == PolyscipStatus::WeightSpacePhase) {
1529  auto untested_weight = weight_space_poly_->getUntestedWeight();
1530  SCIP_CALL(setWeightedObjective(untested_weight));
1531  SCIP_CALL(solve());
1532  auto scip_status = SCIPgetStatus(scip_);
1533  if (scip_status == SCIP_STATUS_INFORUNBD && !is_sub_prob_) {
1534  scip_status = separateINFORUNBD(untested_weight);
1535  }
1536  if (scip_status == SCIP_STATUS_OPTIMAL) {
1537  auto res = getOptimalResult();
1538  auto weighted_outcome = std::inner_product(begin(untested_weight),
1539  end(untested_weight),
1540  begin(res.second),
1541  0.);
1542  if (weighted_outcome + cmd_line_args_.getEpsilon() <
1543  weight_space_poly_->getUntestedVertexWOV(untested_weight)) {
1544  weight_space_poly_->incorporateNewOutcome(scip_,
1545  untested_weight,
1546  res.second);
1547  bounded_.push_back(std::move(res));
1548  }
1549  else {
1550  weight_space_poly_->incorporateKnownOutcome(untested_weight);
1551  }
1552  }
1553  else if (scip_status == SCIP_STATUS_UNBOUNDED) {
1554  SCIP_CALL(handleUnboundedStatus()); //adds unbounded result to unbounded_
1555  weight_space_poly_->incorporateNewOutcome(scip_,
1556  untested_weight,
1557  unbounded_.back().second, // was added by handleUnboundedStatus()
1558  true);
1559  }
1560  else {
1561  SCIP_CALL(handleNonOptNonUnbdStatus(scip_status));
1562  }
1563  }
1564  if (polyscip_status_ == PolyscipStatus::WeightSpacePhase) {
1565  polyscip_status_ = PolyscipStatus::Finished;
1566  }
1567  return SCIP_OKAY;
1568  }
1569 
1570  /**
1571  * Print results
1572  * @param os Output stream to print to
1573  */
1574  void Polyscip::printResults(ostream &os) const {
1575  auto new_line = false;
1576  for (const auto& result : bounded_) {
1577  if (cmd_line_args_.outputOutcomes()) {
1578  outputOutcome(result.second, os);
1579  new_line = true;
1580  }
1581  if (cmd_line_args_.outputSols()) {
1582  printSol(result.first, os);
1583  new_line = true;
1584  }
1585  if (new_line) {
1586  os << "\n";
1587  }
1588  }
1589  for (const auto& result : unbounded_) {
1590  if (cmd_line_args_.outputOutcomes()) {
1591  outputOutcome(result.second, os, "Ray = ");
1592  new_line = true;
1593  }
1594  if (cmd_line_args_.outputSols()) {
1595  printSol(result.first, os);
1596  new_line = true;
1597  }
1598  if (new_line) {
1599  os << "\n";
1600  }
1601  }
1602  }
1603 
1604  /**
1605  * Print solution
1606  * @param sol Solution to print
1607  * @param os Output stream to write to
1608  */
1609  void Polyscip::printSol(const SolType& sol, ostream& os) const {
1610  for (const auto& elem : sol)
1611  os << elem.first << "=" << elem.second << " ";
1612  }
1613 
1614  /**
1615  * Print outcome
1616  * @param outcome Outcome to print
1617  * @param os Output stream to write to
1618  * @param desc Description to print before given outcome
1619  */
1620  void Polyscip::outputOutcome(const OutcomeType &outcome, std::ostream &os, const std::string desc) const {
1621  if (obj_sense_ == SCIP_OBJSENSE_MAXIMIZE) {
1622  global::print(outcome, desc + "[ ", "] ", os, true);
1623  }
1624  else {
1625  global::print(outcome, desc + "[ ", "] ", os);
1626  }
1627  }
1628 
1629  /**
1630  * Check whether file can be opened
1631  * @param filename Name of file to open
1632  * @return true if corresponding file can be opened; false otherwise
1633  */
1634  bool Polyscip::filenameIsOkay(const string& filename) {
1635  std::ifstream file(filename.c_str());
1636  return file.good();
1637  }
1638 
1639  /**
1640  * Print objective
1641  * @param obj_no Corresponding index of objective
1642  * @param nonzero_indices Indices of variables with non-zero coefficients
1643  * @param nonzero_vals Corresponding non-zero coefficient variable values
1644  * @param os Output stream to write to
1645  */
1646  void Polyscip::printObjective(size_t obj_no,
1647  const std::vector<int>& nonzero_indices,
1648  const std::vector<SCIP_Real>& nonzero_vals,
1649  ostream& os) const {
1650  assert (!nonzero_indices.empty());
1651  auto size = nonzero_indices.size();
1652  assert (size == nonzero_vals.size());
1653  auto obj = vector<SCIP_Real>(global::narrow_cast<size_t>(SCIPgetNOrigVars(scip_)), 0);
1654  for (size_t i=0; i<size; ++i)
1655  obj[nonzero_indices[i]] = nonzero_vals[i];
1656  global::print(obj, std::to_string(obj_no) + ". obj: [", "]", os);
1657  os << "\n";
1658  }
1659 
1660  /**
1661  * Indicates whether objective given by index is redundant
1662  * @param begin_nonzeros begin_nonzeros[i+1] = begin_nonzeros[i] + obj_probdata->getNumberNonzeroCoeffs(i)
1663  * @param obj_to_nonzero_indices indices of non-zero variables for each objective
1664  * @param obj_to_nonzero_values non-zero variables for each objective
1665  * @param index index of objective to check
1666  * @return true if checked objective is redundant; false otherwise
1667  */
1668  bool Polyscip::objIsRedundant(const vector<int>& begin_nonzeros,
1669  const vector< vector<int> >& obj_to_nonzero_indices,
1670  const vector< vector<SCIP_Real> >& obj_to_nonzero_values,
1671  size_t checked_obj) const {
1672  bool is_redundant = false;
1673  auto obj_probdata = dynamic_cast<ProbDataObjectives*>(SCIPgetObjProbData(scip_));
1674 
1675  assert (obj_probdata != nullptr);
1676  assert (checked_obj >= 1 && checked_obj < obj_probdata->getNoObjs());
1677 
1678  SCIP_LPI* lpi;
1679  auto retcode = SCIPlpiCreate(addressof(lpi), nullptr, "check objective redundancy", SCIP_OBJSEN_MINIMIZE);
1680  if (retcode != SCIP_OKAY)
1681  throw std::runtime_error("no SCIP_OKAY for SCIPlpiCreate\n.");
1682 
1683  auto no_cols = global::narrow_cast<int>(checked_obj);
1684  auto obj = vector<SCIP_Real>(checked_obj, 1.);
1685  auto lb = vector<SCIP_Real>(checked_obj, 0.);
1686  auto ub = vector<SCIP_Real>(checked_obj, SCIPlpiInfinity(lpi));
1687  auto no_nonzero = begin_nonzeros.at(checked_obj);
1688 
1689  auto beg = vector<int>(begin(begin_nonzeros), begin(begin_nonzeros)+checked_obj);
1690  auto ind = vector<int>{};
1691  ind.reserve(global::narrow_cast<size_t>(no_nonzero));
1692  auto val = vector<SCIP_Real>{};
1693  val.reserve(global::narrow_cast<size_t>(no_nonzero));
1694  for (size_t i=0; i<checked_obj; ++i) {
1695  ind.insert(end(ind), begin(obj_to_nonzero_indices[i]), end(obj_to_nonzero_indices[i]));
1696  val.insert(end(val), begin(obj_to_nonzero_values[i]), end(obj_to_nonzero_values[i]));
1697  }
1698 
1699  auto no_rows = SCIPgetNOrigVars(scip_);
1700  auto vars = SCIPgetOrigVars(scip_);
1701  auto lhs = vector<SCIP_Real>(global::narrow_cast<size_t>(no_rows), 0.);
1702  for (auto i=0; i<no_rows; ++i)
1703  lhs[i] = obj_probdata->getObjCoeff(vars[i], checked_obj);
1704 
1705  retcode = SCIPlpiLoadColLP(lpi,
1707  no_cols,
1708  obj.data(),
1709  lb.data(),
1710  ub.data(),
1711  nullptr,
1712  no_rows,
1713  lhs.data(),
1714  lhs.data(),
1715  nullptr,
1716  no_nonzero,
1717  beg.data(),
1718  ind.data(),
1719  val.data());
1720 
1721  if (retcode != SCIP_OKAY)
1722  throw std::runtime_error("no SCIP_OKAY for SCIPlpiLoadColLP\n");
1723 
1724  retcode = SCIPlpiSolvePrimal(lpi);
1725  if (retcode != SCIP_OKAY)
1726  throw std::runtime_error("no SCIP_OKAY for SCIPlpiSolvePrimal\n");
1727 
1728  if (SCIPlpiIsPrimalFeasible(lpi)) {
1729  is_redundant = true;
1730  }
1731  else {
1732  assert (SCIPlpiIsPrimalInfeasible(lpi));
1733  }
1734 
1735  retcode = SCIPlpiFree(addressof(lpi));
1736  if (retcode != SCIP_OKAY)
1737  throw std::runtime_error("no SCIP_OKAY for SCIPlpiFree\n");
1738 
1739  return is_redundant;
1740  }
1741 
1742  /**
1743  * Read multi-objective problem file
1744  * @return SCIP_OKAY if everything worked; otherwise a suitable error code is passed
1745  */
1747  if (polyscip_status_ != PolyscipStatus::Unsolved) {
1748  return SCIP_OKAY;
1749  }
1750  auto filename = cmd_line_args_.getProblemFile();
1751  SCIP_CALL( SCIPreadProb(scip_, filename.c_str(), "mop") );
1752  auto obj_probdata = dynamic_cast<ProbDataObjectives*>(SCIPgetObjProbData(scip_));
1753  assert (obj_probdata != nullptr);
1754  no_objs_ = obj_probdata->getNoObjs();
1755 
1756  if (cmd_line_args_.onlyExtremal() || SCIPgetNOrigContVars(scip_) == SCIPgetNOrigVars(scip_)) {
1757  only_weight_space_phase_ = true;
1758  }
1759 
1760  auto vars = SCIPgetOrigVars(scip_);
1761  auto begin_nonzeros = vector<int>(no_objs_, 0);
1762  for (size_t i = 0; i < no_objs_ - 1; ++i)
1763  begin_nonzeros[i + 1] = global::narrow_cast<int>(
1764  begin_nonzeros[i] + obj_probdata->getNumberNonzeroCoeffs(i));
1765 
1766  auto obj_to_nonzero_inds = vector<vector<int> >{};
1767  auto obj_to_nonzero_vals = vector<vector<SCIP_Real> >{};
1768  for (size_t obj_ind = 0; obj_ind < no_objs_; ++obj_ind) {
1769  auto nonzero_vars = obj_probdata->getNonZeroCoeffVars(obj_ind);
1770  auto size = nonzero_vars.size();
1771  if (size == 0) {
1772  cout << obj_ind << ". objective is zero objective!\n";
1773  polyscip_status_ = PolyscipStatus::Error;
1774  return SCIP_OKAY;
1775  }
1776  auto nonzero_inds = vector<int>(size, 0);
1777  std::transform(begin(nonzero_vars),
1778  end(nonzero_vars),
1779  begin(nonzero_inds),
1780  [](SCIP_VAR *var) { return SCIPvarGetProbindex(var); });
1781  std::sort(begin(nonzero_inds), end(nonzero_inds));
1782 
1783  auto nonzero_vals = vector<SCIP_Real>(size, 0.);
1784  std::transform(begin(nonzero_inds),
1785  end(nonzero_inds),
1786  begin(nonzero_vals),
1787  [&](int var_ind) { return obj_probdata->getObjCoeff(vars[var_ind], obj_ind); });
1788 
1789 
1790  if (cmd_line_args_.beVerbose())
1791  printObjective(obj_ind, nonzero_inds, nonzero_vals);
1792 
1793  obj_to_nonzero_inds.push_back(std::move(nonzero_inds)); // nonzero_inds invalid from now on
1794  obj_to_nonzero_vals.push_back(std::move(nonzero_vals)); // nonzero_vals invalid from now on
1795 
1796  if (obj_ind > 0 && objIsRedundant(begin_nonzeros, // first objective is always non-redundant
1797  obj_to_nonzero_inds,
1798  obj_to_nonzero_vals,
1799  obj_ind)) {
1800  cout << obj_ind << ". objective is non-negative linear combination of previous objectives!\n";
1801  cout << "Only problems with non-redundant objectives will be solved.\n";
1802  polyscip_status_ = PolyscipStatus::Error;
1803  return SCIP_OKAY;
1804  }
1805  }
1806 
1807  if (SCIPgetObjsense(scip_) == SCIP_OBJSENSE_MAXIMIZE) {
1808  obj_sense_ = SCIP_OBJSENSE_MAXIMIZE;
1809  // internally we treat problem as min problem and negate objective coefficients
1811  obj_probdata->negateAllCoeffs();
1812  }
1813  if (cmd_line_args_.beVerbose()) {
1814  cout << "Objective sense: ";
1815  if (obj_sense_ == SCIP_OBJSENSE_MAXIMIZE)
1816  cout << "MAXIMIZE\n";
1817  else
1818  cout << "MINIMIZE\n";
1819  cout << "Number of objectives: " << no_objs_ << "\n";
1820  }
1821  polyscip_status_ = PolyscipStatus::ProblemRead;
1822  return SCIP_OKAY;
1823  }
1824 
1825  /**
1826  * Write results to file named 'solutions_name-of-problem-file.txt'
1827  */
1829  auto prob_file = cmd_line_args_.getProblemFile();
1830  size_t prefix = prob_file.find_last_of("/"), //separate path/ and filename.mop
1831  suffix = prob_file.find_last_of("."), //separate filename and .mop
1832  start_ind = (prefix == string::npos) ? 0 : prefix + 1,
1833  end_ind = (suffix != string::npos) ? suffix : string::npos;
1834  string file_name = "solutions_" +
1835  prob_file.substr(start_ind, end_ind - start_ind) + ".txt";
1836  auto write_path = cmd_line_args_.getWritePath();
1837  if (write_path.back() != '/')
1838  write_path.push_back('/');
1839  std::ofstream solfs(write_path + file_name);
1840  if (solfs.is_open()) {
1841  printResults(solfs);
1842  solfs.close();
1843  cout << "#Solution file " << file_name
1844  << " written to: " << write_path << "\n";
1845  }
1846  else
1847  std::cerr << "ERROR writing solution file\n.";
1848  }
1849 
1850  /**
1851  * Checks whether results corresponding to given iterator is dominated or equal to other given elements
1852  * @param it Const_iterator corresponding to result
1853  * @param beg_it Const_iterator to beginning of result container
1854  * @param end_it Const_iterator to past-the-end of result container
1855  * @return true if result given by it is dominated or equal to other given results; false otherwise
1856  */
1857  bool Polyscip::isDominatedOrEqual(ResultContainer::const_iterator it,
1858  ResultContainer::const_iterator beg_it,
1859  ResultContainer::const_iterator end_it) const {
1860  for (auto curr = beg_it; curr != end_it; ++curr) {
1861  if (it == curr)
1862  continue;
1863  else if (std::equal(begin(curr->second),
1864  end(curr->second),
1865  begin(it->second),
1866  std::less_equal<ValueType>())) {
1867  outputOutcome(it->second, std::cout, "outcome ");
1868  outputOutcome(curr->second, std::cout, " dominated by ");
1869  return true;
1870  }
1871  }
1872  return false;
1873  }
1874 
1875 
1876  /**
1877  * Indicates whether dominated results were computed
1878  * @return true if dominated bounded results were computed
1879  */
1881  for (auto cur=begin(bounded_); cur!=end(bounded_); ++cur) {
1882  if (isDominatedOrEqual(cur, begin(bounded_), end(bounded_)))
1883  return true;
1884  }
1885  return false;
1886  }
1887 
1888 
1889 }
bool outputOutcomes() const
Definition: cmd_line_args.h:71
SCIP_RETCODE SCIPlpiFree(SCIP_LPI **lpi)
Definition: lpi_clp.cpp:634
NondomProjections(double epsilon, const ResultContainer &supported, std::size_t first, std::size_t second)
Definition: polyscip.cpp:134
std::vector< ValueType > WeightType
Type for weight vectors.
SCIP_RETCODE SCIPcreateConsBasicLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs)
void writeResultsToFile() const
Definition: polyscip.cpp:1828
bool epsilonDominates(const TwoDProj &lhs, const TwoDProj &rhs) const
Definition: polyscip.cpp:174
std::pair< ValueType, ValueType > Interval
Interval I = [a,b)
Definition: polyscip.h:186
SCIP_Real SCIPgetPrimalbound(SCIP *scip)
Class storing multiple objectives of given problem instance.
bool dominatedPointsFound() const
Definition: polyscip.cpp:1880
SCIP_RETCODE SCIPlpiSolvePrimal(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:1796
SCIP_RETCODE SCIPdelCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2841
int SCIPgetNOrigVars(SCIP *scip)
Definition: scip_prob.c:2430
SCIP_Real ValueType
Type for computed values.
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
bool hasParameterFile() const
Definition: cmd_line_args.h:89
bool beVerbose() const
Definition: cmd_line_args.h:53
Class for .mop file reader.
Definition: ReaderMOP.h:38
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1245
SCIP_RETCODE SCIPstopClock(SCIP *scip, SCIP_CLOCK *clck)
Definition: scip_timing.c:169
#define FALSE
Definition: def.h:87
bool isFeasible(double epsilon) const
Definition: polyscip.cpp:340
std::vector< OutcomeType > getNondomProjOutcomes() const
Definition: polyscip.cpp:208
SCIP_Real SCIPinfinity(SCIP *scip)
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
Status if an error occured.
SCIP_RETCODE SCIPsetPresolving(SCIP *scip, SCIP_PARAMSETTING paramsetting, SCIP_Bool quiet)
Definition: scip_param.c:923
int SCIPvarGetProbindex(SCIP_VAR *var)
Definition: var.c:17600
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition: scip_var.c:185
SCIP_VAR ** SCIPgetOrigVars(SCIP *scip)
Definition: scip_prob.c:2403
ValueType getFirst() const
Definition: polyscip.h:64
General types used for PolySCIP.
SCIP_RETCODE SCIPfreeClock(SCIP *scip, SCIP_CLOCK **clck)
Definition: scip_timing.c:118
Target narrow_cast(Source v)
doubledescription::DoubleDescriptionMethod DDMethod
abbreviation
Definition: polyscip.cpp:73
std::vector< ValueType > OutcomeType
Type for points, rays in objective space.
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:283
SCIP_Bool SCIPisTransformed(SCIP *scip)
Definition: scip_general.c:566
SCIP_RETCODE SCIPsetRealParam(SCIP *scip, const char *name, SCIP_Real value)
Definition: scip_param.c:594
SCIP_RETCODE SCIPlpiCreate(SCIP_LPI **lpi, SCIP_MESSAGEHDLR *messagehdlr, const char *name, SCIP_OBJSEN objsen)
Definition: lpi_clp.cpp:522
std::pair< SolType, OutcomeType > Result
A result comprises of a solution/ray in feasible space and a corresponding outcome in objective space...
real eps
std::size_t numberOfBoundedResults() const
Definition: polyscip.cpp:1270
std::size_t numberofUnboundedResults() const
Definition: polyscip.cpp:1278
SCIP_VAR * w
Definition: circlepacking.c:58
Status while results of weight space polyhedron are computed.
SCIP_RETCODE SCIPdelVar(SCIP *scip, SCIP_VAR *var, SCIP_Bool *deleted)
Definition: scip_prob.c:1788
C++ wrapper for default SCIP plugins.
SCIP_RETCODE SCIPsetObjsense(SCIP *scip, SCIP_OBJSENSE objsense)
Definition: scip_prob.c:1240
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2613
SCIP_RETCODE SCIPcreateClock(SCIP *scip, SCIP_CLOCK **clck)
Definition: scip_timing.c:67
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2768
bool isSupersetOf(const RectangularBox &other) const
Definition: polyscip.cpp:296
Status after lexicographic optimal results were computed.
SCIP_RETCODE SCIPreadProb(SCIP *scip, const char *filename, const char *extension)
Definition: scip_prob.c:329
Status if problem was solved successfully.
RectangularBox(const std::vector< Interval > &box)
Definition: polyscip.cpp:230
SCIP_STATUS SCIPgetStatus(SCIP *scip)
Definition: scip_general.c:474
SCIP_Bool SCIPlpiIsPrimalInfeasible(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:2488
SCIP_RETCODE readProblem()
Definition: polyscip.cpp:1746
void print(const Container &container, const std::string &prefix="", const std::string &suffix="", std::ostream &os=std::cout, bool negate=false, int prec=6)
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17251
bool isSubsetOf(const RectangularBox &other) const
Definition: polyscip.cpp:310
#define SCIP_CALL(x)
Definition: def.h:384
Class representing the 1-skeleton of the weight space polyhedron.
friend std::ostream & operator<<(std::ostream &os, const RectangularBox &box)
Definition: polyscip.cpp:284
TimeLimitType getTimeLimit() const
friend std::ostream & operator<<(std::ostream &os, const TwoDProj &proj)
SCIP_Bool SCIPhasPrimalRay(SCIP *scip)
Definition: scip_sol.c:3535
SCIP_RETCODE SCIPchgVarObj(SCIP *scip, SCIP_VAR *var, SCIP_Real newobj)
Definition: scip_var.c:4510
bool hasTimeLimit() const
Definition: cmd_line_args.h:83
SCIP_RETCODE SCIPfreeTransform(SCIP *scip)
Definition: scip_solve.c:3432
#define SCIP_Bool
Definition: def.h:84
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
Command line arguments for PolySCIP.
Definition: cmd_line_args.h:34
.mop file format reader
SCIP_Real SCIPlpiInfinity(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:3905
std::vector< Result > ResultContainer
Container type for results.
enum SCIP_Status SCIP_STATUS
Definition: type_stat.h:58
Status after problem instance was read successfully.
SCIP_Real SCIPgetClockTime(SCIP *scip, SCIP_CLOCK *clck)
Definition: scip_timing.c:310
Global helper functions.
double getEpsilon() const
double getDelta() const
Status while computing 2-dimensional non-dominated projection results.
Interval getInterval(std::size_t index) const
Definition: polyscip.cpp:273
SCIP_RETCODE SCIPfreeSol(SCIP *scip, SCIP_SOL **sol)
Definition: scip_sol.c:976
SCIP_RETCODE SCIPincludeObjReader(SCIP *scip, scip::ObjReader *objreader, SCIP_Bool deleteobject)
Definition: objreader.cpp:146
bool outputSols() const
Definition: cmd_line_args.h:77
int roundWeightedObjCoeff() const
Definition: cmd_line_args.h:95
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1435
void printResults(std::ostream &os=std::cout) const
Definition: polyscip.cpp:1574
SCIP_RETCODE SCIPchgRhsLinear(SCIP *scip, SCIP_CONS *cons, SCIP_Real rhs)
Class representing a two-dimensional projection of an outcome.
Definition: polyscip.h:47
Initial status after calling public constructor.
int SCIPgetNOrigContVars(SCIP *scip)
Definition: scip_prob.c:2538
scip::ObjProbData * SCIPgetObjProbData(SCIP *scip)
Stores coefficients and basic methods for objectives of given multi-objective problem.
Double description method for transforming a polyhedron given via its h-representation into its v-rep...
PolySCIP solver class.
SCIP_VAR ** b
Definition: circlepacking.c:56
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2304
Status if given time limit was reached.
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1666
SCIP_RETCODE SCIPlpiLoadColLP(SCIP_LPI *lpi, SCIP_OBJSEN objsen, int ncols, const SCIP_Real *obj, const SCIP_Real *lb, const SCIP_Real *ub, char **colnames, int nrows, const SCIP_Real *lhs, const SCIP_Real *rhs, char **rownames, int nnonz, const int *beg, const int *ind, const SCIP_Real *val)
Definition: lpi_clp.cpp:668
SCIP_RETCODE computeNondomPoints()
Definition: polyscip.cpp:526
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1110
SCIP_RETCODE SCIPreadParams(SCIP *scip, const char *filename)
Definition: scip_param.c:742
std::string getProblemFile() const
SCIP_VAR * a
Definition: circlepacking.c:57
std::pair< std::size_t, std::size_t > ObjPair
Pair of objectives indices.
Definition: polyscip.h:312
Class representing non-dominated projections.
Definition: polyscip.h:89
SCIP_Bool SCIPlpiIsPrimalFeasible(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:2507
bool isDominated(const OutcomeType &outcome) const
Definition: polyscip.cpp:353
ValueType getSecond() const
Definition: polyscip.h:70
bool isDisjointFrom(const RectangularBox &other) const
Definition: polyscip.cpp:324
TwoDProj getLastProj() const
Definition: polyscip.h:162
Polyscip(int argc, const char *const *argv)
Definition: polyscip.cpp:414
A rectangular box R = [a_1,e_1) x ... x [a_k,e_k) is a k-ary Cartesian product of half-open intervals...
Definition: polyscip.h:184
SCIP_OBJSENSE SCIPgetObjsense(SCIP *scip)
Definition: scip_prob.c:1223
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
std::string getParameterFile() const
PolyscipStatus getStatus() const
Definition: polyscip.cpp:1262
SCIP_RETCODE SCIPchgLhsLinear(SCIP *scip, SCIP_CONS *cons, SCIP_Real lhs)
void printStatus(std::ostream &os=std::cout) const
Definition: polyscip.cpp:492
bool onlyExtremal() const
Definition: cmd_line_args.h:59
std::vector< std::pair< std::string, ValueType > > SolType
Type for solutions in feasible space.
#define SCIP_CALL_ABORT(x)
Definition: def.h:363
std::string getWritePath() const
SCIP_RETCODE SCIPstartClock(SCIP *scip, SCIP_CLOCK *clck)
Definition: scip_timing.c:152
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1352
std::size_t size() const
Definition: polyscip.cpp:264
PolySCIP command line arguments.
Algorithm for transforming h-representation to v-representation.
SCIP_Real SCIPgetPrimalRayVal(SCIP *scip, SCIP_VAR *var)
Definition: scip_sol.c:3553
std::vector< RectangularBox > getDisjointPartsFrom(double delta, const RectangularBox &other) const
Definition: polyscip.cpp:383
SCIP callable library.
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:315
SCIP_RETCODE SCIPcreateFiniteSolCopy(SCIP *scip, SCIP_SOL **sol, SCIP_SOL *sourcesol, SCIP_Bool *success)
Definition: scip_sol.c:840
TwoDProj(const OutcomeType &outcome, std::size_t first, std::size_t second)
Definition: polyscip.cpp:81