1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program PolySCIP */ 4 /* */ 5 /* Copyright (C) 2012-2021 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 45 #include "double_description_method.h" 46 #include "scip/scip.h" 47 #include "objscip/objscipdefplugins.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" 53 #include "weight_space_polyhedron.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 73 using DDMethod = doubledescription::DoubleDescriptionMethod; ///< abbreviation 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 */ TwoDProj(const OutcomeType & outcome,size_t first,size_t second)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 */ operator <<(ostream & os,const TwoDProj & p)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 */ operator <<(ostream & os,const NondomProjections & nd)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 */ add(TwoDProj proj,Result res)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 */ NondomProjections(double eps,const ResultContainer & supported,size_t first,size_t second)134 NondomProjections::NondomProjections(double eps, 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 */ epsilonDominates(const TwoDProj & lhs,const TwoDProj & rhs) const174 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 */ update()181 void NondomProjections::update() { 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 */ update(TwoDProj proj,Result res)191 void NondomProjections::update(TwoDProj proj, Result res) { 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 */ getNondomProjOutcomes() const208 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 */ finished() const221 bool NondomProjections::finished() const { 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 */ RectangularBox(const std::vector<Interval> & box)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 */ RectangularBox(std::vector<Interval> && box)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 */ RectangularBox(std::vector<Interval>::const_iterator first_beg,std::vector<Interval>::const_iterator first_end,Interval second,std::vector<Interval>::const_iterator third_beg,std::vector<Interval>::const_iterator third_end)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 */ size() const264 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 */ getInterval(size_t index) const273 RectangularBox::Interval RectangularBox::getInterval(size_t index) const { 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 */ operator <<(std::ostream & os,const RectangularBox & box)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 */ isSupersetOf(const RectangularBox & other) const296 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 */ isSubsetOf(const RectangularBox & other) const310 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 */ isDisjointFrom(const RectangularBox & other) const324 bool RectangularBox::isDisjointFrom(const RectangularBox &other) const { 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 */ isFeasible(double epsilon) const340 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 */ isDominated(const OutcomeType & outcome) const353 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 */ getIntervalIntersection(std::size_t index,const RectangularBox & other) const369 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 */ getDisjointPartsFrom(double delta,const RectangularBox & other) const383 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 */ Polyscip(int argc,const char * const * argv)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); 426 SCIPincludeDefaultPlugins(scip_); 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 */ Polyscip(const CmdLineArgs & cmd_line_args,SCIP * scip,std::size_t no_objs,SCIP_CLOCK * clock_total)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 */ ~Polyscip()481 Polyscip::~Polyscip() { 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 */ printStatus(std::ostream & os) const492 void Polyscip::printStatus(std::ostream& os) const { 493 switch(polyscip_status_) { 494 case PolyscipStatus::TwoProjPhase: 495 os << "PolySCIP Status: ComputeProjectedNondomPointsPhase\n"; 496 break; 497 case PolyscipStatus::Error: 498 os << "PolySCIP Status: Error\n"; 499 break; 500 case PolyscipStatus::Finished: 501 os << "PolySCIP Status: Successfully finished\n"; 502 break; 503 case PolyscipStatus::LexOptPhase: 504 os << "PolySCIP Status: LexOptPhase\n"; 505 break; 506 case PolyscipStatus::ProblemRead: 507 os << "PolySCIP Status: ProblemRead\n"; 508 break; 509 case PolyscipStatus::TimeLimitReached: 510 os << "PolySCIP Status: TimeLimitReached\n"; 511 break; 512 case PolyscipStatus::Unsolved: 513 os << "PolySCIP Status: Unsolved\n"; 514 break; 515 case PolyscipStatus::WeightSpacePhase: 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 */ computeNondomPoints()526 SCIP_RETCODE Polyscip::computeNondomPoints() { 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 */ computeLexicographicOptResults(vector<vector<SCIP_VAR * >> & orig_vars,vector<vector<ValueType>> & orig_vals)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 */ computeLexicographicOptResult(size_t considered_obj,vector<vector<SCIP_VAR * >> & orig_vars,vector<vector<ValueType>> & orig_vals)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 604 SCIP_CALL(SCIPfreeTransform(scip_)); 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 */ computeBoundedNondomResultsForUnbdObjs()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 */ computeNonLexicographicNondomResults(const vector<vector<SCIP_VAR * >> & orig_vars,const vector<vector<ValueType>> & orig_vals)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 */ boxResultIsDominated(const OutcomeType & outcome,const vector<vector<SCIP_VAR * >> & orig_vars,const vector<vector<ValueType>> & orig_vals)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 773 SCIP_CALL_ABORT( SCIPfreeTransform(scip_) ); 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 808 SCIP_CALL_ABORT( SCIPfreeTransform(scip_) ); 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 */ computeNondomPointsInBox(const RectangularBox & box,const vector<vector<SCIP_VAR * >> & orig_vars,const vector<vector<ValueType>> & orig_vals)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_)) { 833 SCIP_CALL_ABORT( SCIPfreeTransform(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_)) { 853 SCIP_CALL_ABORT( SCIPfreeTransform(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 */ computeDisjointBoxes(list<RectangularBox> && feasible_boxes) const886 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 */ computeFeasibleBoxes(const map<ObjPair,vector<OutcomeType>> & proj_nd_outcomes,const vector<vector<SCIP_VAR * >> & orig_vars,const vector<vector<ValueType>> & orig_vals)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 */ createNewVarTransformCons(SCIP_VAR * new_var,const vector<SCIP_VAR * > & orig_vars,const vector<ValueType> & orig_vals,const ValueType & rhs,const ValueType & beta_i)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 996 SCIPcreateConsBasicLinear(scip_, 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 */ createObjValCons(const vector<SCIP_VAR * > & vars,const vector<ValueType> & vals,const ValueType & lhs,const ValueType & rhs)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); 1023 SCIPcreateConsBasicLinear(scip_, 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 */ computeNondomProjResult(SCIP_CONS * cons1,SCIP_CONS * cons2,ValueType rhs_cons1,ValueType rhs_cons2,size_t obj_1,size_t obj_2,ResultContainer & results)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 */ solveWeightedTchebycheff(const vector<vector<SCIP_VAR * >> & orig_vars,const vector<vector<ValueType>> & orig_vals,size_t obj_1,size_t obj_2)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, 1128 SCIP_VARTYPE_CONTINUOUS); 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 1205 SCIP_CALL_ABORT( SCIPfreeTransform(scip_) ); 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 1240 SCIP_CALL_ABORT( SCIPfreeTransform(scip_) ); 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 */ getStatus() const1262 Polyscip::PolyscipStatus Polyscip::getStatus() const { 1263 return polyscip_status_; 1264 } 1265 1266 /** 1267 * Get number of bounded results 1268 * @return Number of computed bounded results 1269 */ numberOfBoundedResults() const1270 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 */ numberofUnboundedResults() const1278 std::size_t Polyscip::numberofUnboundedResults() const { 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 */ separateINFORUNBD(const WeightType & weight,bool with_presolving)1288 SCIP_STATUS Polyscip::separateINFORUNBD(const WeightType& weight, bool with_presolving) { 1289 if (!with_presolving) 1290 SCIPsetPresolving(scip_, SCIP_PARAMSETTING_OFF, TRUE); 1291 auto zero_weight = WeightType(no_objs_, 0.); 1292 setWeightedObjective(zero_weight); 1293 solve(); // re-compute with zero objective 1294 if (!with_presolving) 1295 SCIPsetPresolving(scip_, SCIP_PARAMSETTING_DEFAULT, TRUE); 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 */ handleNonOptNonUnbdStatus(SCIP_STATUS status)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 */ handleUnboundedStatus(bool check_if_new_result)1341 SCIP_RETCODE Polyscip::handleUnboundedStatus(bool check_if_new_result) { 1342 if (!SCIPhasPrimalRay(scip_)) { 1343 SCIP_CALL( SCIPsetPresolving(scip_, SCIP_PARAMSETTING_OFF, TRUE) ); 1344 if (SCIPisTransformed(scip_)) 1345 SCIP_CALL( SCIPfreeTransform(scip_) ); 1346 SCIP_CALL( solve() ); 1347 SCIP_CALL( SCIPsetPresolving(scip_, SCIP_PARAMSETTING_DEFAULT, TRUE) ); 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 */ getResult(bool outcome_is_bounded,SCIP_SOL * primal_sol)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 */ getOptimalResult()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 */ outcomeIsNew(const OutcomeType & outcome,bool outcome_is_bounded) const1429 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 */ outcomeIsNew(const OutcomeType & outcome,ResultContainer::const_iterator beg,ResultContainer::const_iterator last) const1443 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 */ outcomesCoincide(const OutcomeType & a,const OutcomeType & b,double epsilon)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 */ solve()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 */ setWeightedObjective(const WeightType & weight)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 */ computeWeightSpaceResults()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 */ printResults(ostream & os) const1574 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 */ printSol(const SolType & sol,ostream & os) const1609 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 */ outputOutcome(const OutcomeType & outcome,std::ostream & os,const std::string desc) const1620 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 */ filenameIsOkay(const string & filename)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 */ printObjective(size_t obj_no,const std::vector<int> & nonzero_indices,const std::vector<SCIP_Real> & nonzero_vals,ostream & os) const1646 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 */ objIsRedundant(const vector<int> & begin_nonzeros,const vector<vector<int>> & obj_to_nonzero_indices,const vector<vector<SCIP_Real>> & obj_to_nonzero_values,size_t checked_obj) const1668 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, 1706 SCIP_OBJSEN_MINIMIZE, 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 */ readProblem()1746 SCIP_RETCODE Polyscip::readProblem() { 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 1810 SCIPsetObjsense(scip_, SCIP_OBJSENSE_MINIMIZE); 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 */ writeResultsToFile() const1828 void Polyscip::writeResultsToFile() const { 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 */ isDominatedOrEqual(ResultContainer::const_iterator it,ResultContainer::const_iterator beg_it,ResultContainer::const_iterator end_it) const1857 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 */ dominatedPointsFound() const1880 bool Polyscip::dominatedPointsFound() const { 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 } 1890