41 #include <type_traits> 66 using std::reference_wrapper;
82 : proj_(outcome.at(first), outcome.at(second))
92 os <<
"Proj = [" << p.proj_.first <<
", " << p.proj_.second <<
"]";
103 os <<
"Nondominated projections: ";
104 for (
const auto& p_pair : nd.nondom_projections_)
105 os << p_pair.first <<
" ";
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_)) {
118 auto ret = nondom_projections_.emplace(std::move(proj),
ResultContainer{std::move(res)});
122 nondom_projections_[proj].push_back(std::move(res));
140 if (lhs.getFirst() + eps < rhs.getFirst())
142 else if (rhs.getFirst() + eps < lhs.getFirst())
145 return lhs.getSecond() < rhs.getSecond();
148 assert (first < second);
149 assert (!supported.empty());
150 for (
const auto& res : supported) {
151 add(
TwoDProj(res.second, first, second), res);
154 auto it = begin(nondom_projections_);
155 while (std::next(it) != end(nondom_projections_)) {
156 auto next = std::next(it);
158 nondom_projections_.erase(next);
164 assert (!nondom_projections_.empty());
165 current_ = begin(nondom_projections_);
182 assert (current_ != std::prev(end(nondom_projections_)) && current_ != end(nondom_projections_));
192 assert (current_ != std::prev(end(nondom_projections_)) && current_ != end(nondom_projections_));
193 auto it = add(proj, res);
195 nondom_projections_.erase(current_);
199 while (std::next(it) != end(nondom_projections_) &&
epsilonDominates(proj, std::next(it)->first)) {
200 nondom_projections_.erase(std::next(it));
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);
222 assert (current_ != end(nondom_projections_));
223 return current_ == std::prev(end(nondom_projections_));
231 : box_(begin(box), end(box))
239 : box_(begin(box), end(box))
251 std::vector<Interval>::const_iterator first_end,
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_));
274 assert (index <
size());
285 for (
auto interval : box.box_)
286 os <<
"[ " << interval.first <<
", " << interval.second <<
" ) ";
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)
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)
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)
341 for (
const auto& elem : box_) {
342 if (elem.first + epsilon > elem.second)
354 assert (outcome.size() == box_.size());
355 for (
size_t i=0; i<box_.size(); ++i) {
356 if (outcome[i] > box_[i].first) {
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};
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) {
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));
397 if (other.box_[i].second < box_[i].second) {
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));
404 intersections.push_back(getIntervalIntersection(i, other));
406 return disjoint_partitions;
415 : cmd_line_args_(argc, argv),
420 clock_total_(nullptr),
421 only_weight_space_phase_(false),
425 assert (scip_ !=
nullptr);
435 cout <<
"Invalid parameter settings file.\n";
439 if (cmd_line_args_.
getDelta() <= 0) {
440 cout <<
"Please choose positive value for parameter --delta .\n";
444 cout <<
"Please choose positive value for parameter --epsilon .\n";
448 cout <<
"Please choose positive value for parameter --timelLimit .\n";
452 cout <<
"Invalid problem file.\n";
468 : cmd_line_args_{cmd_line_args},
473 clock_total_{clock_total},
474 only_weight_space_phase_{
false},
493 switch(polyscip_status_) {
495 os <<
"PolySCIP Status: ComputeProjectedNondomPointsPhase\n";
498 os <<
"PolySCIP Status: Error\n";
501 os <<
"PolySCIP Status: Successfully finished\n";
504 os <<
"PolySCIP Status: LexOptPhase\n";
507 os <<
"PolySCIP Status: ProblemRead\n";
510 os <<
"PolySCIP Status: TimeLimitReached\n";
513 os <<
"PolySCIP Status: Unsolved\n";
516 os <<
"PolySCIP Status: WeightSpacePhase\n";
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));
544 SCIP_CALL( computeLexicographicOptResults(nonzero_orig_vars, nonzero_orig_vals) );
547 if (no_objs_ > 3 || only_weight_space_phase_) {
551 SCIP_CALL(computeNonLexicographicNondomResults(nonzero_orig_vars, nonzero_orig_vals));
565 SCIP_RETCODE Polyscip::computeLexicographicOptResults(vector<vector<SCIP_VAR*>>& orig_vars,
566 vector<vector<ValueType>>& orig_vals) {
568 for (
size_t obj=0; obj<no_objs_; ++obj) {
570 SCIP_CALL(computeLexicographicOptResult(obj, orig_vars, 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*>{};
592 for (
size_t counter = 0; counter<no_objs_; ++counter) {
593 weight[current_obj] = 1;
594 SCIP_CALL( setWeightedObjective(weight) );
598 scip_status = separateINFORUNBD(weight);
601 if (counter < no_objs_-1) {
605 auto cons = createObjValCons(orig_vars[current_obj],
606 orig_vals[current_obj],
610 obj_val_cons.push_back(cons);
614 assert (current_obj == considered_obj);
615 SCIP_CALL( handleUnboundedStatus(
true) );
616 unbd_orig_objs_.push_back(considered_obj);
624 assert (current_obj == 0);
632 weight[current_obj] = 0;
633 current_obj = (current_obj+1) % no_objs_;
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));
644 for (
auto cons : obj_val_cons) {
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;
670 auto bd_result = getOptimalResult();
671 if (outcomeIsNew(bd_result.second, begin(bounded_), end(bounded_))) {
672 bounded_.push_back(std::move(bd_result));
695 SCIP_RETCODE Polyscip::computeNonLexicographicNondomResults(
const vector<vector<SCIP_VAR *>> &orig_vars,
696 const vector<vector<ValueType>> &orig_vals) {
698 if (unboundedResultsExist()) {
699 SCIP_CALL( computeBoundedNondomResultsForUnbdObjs() );
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) {
709 auto proj_nondom_outcomes = solveWeightedTchebycheff(orig_vars,
713 proj_nondom_outcomes_map.insert({
ObjPair(obj_1, obj_2), proj_nondom_outcomes});
718 auto feasible_boxes = computeFeasibleBoxes(proj_nondom_outcomes_map,
721 auto disjoint_boxes = computeDisjointBoxes(std::move(feasible_boxes));
722 assert (feasible_boxes.size() <= disjoint_boxes.size());
724 for (
const auto& box : disjoint_boxes) {
725 if (std::any_of(begin(bounded_),
727 [&box](
const Result& res){
return box.isDominated(res.second);})) {
731 cout <<
"Checking box = " << box <<
"\n";
733 auto new_res = computeNondomPointsInBox(box,
736 for (
const auto &res : new_res) {
741 bounded_.push_back(res);
744 if (!boxResultIsDominated(res.second, orig_vars, orig_vals)) {
745 bounded_.push_back(std::move(res));
764 bool Polyscip::boxResultIsDominated(
const OutcomeType& outcome,
765 const vector<vector<SCIP_VAR*>>& orig_vars,
766 const vector<vector<ValueType>>& orig_vals) {
768 auto size = outcome.size();
769 assert (size == orig_vars.size());
770 assert (size == orig_vals.size());
771 auto is_dominated =
false;
775 auto new_cons = vector<SCIP_CONS*>{};
777 for (
size_t i=0; i<size; ++i) {
778 auto cons = createObjValCons(orig_vars[i],
783 new_cons.push_back(cons);
787 setWeightedObjective(weight);
793 assert (weight.size() == outcome.size());
810 for (
auto cons : new_cons) {
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());
831 auto obj_val_cons = vector<SCIP_CONS *>{};
835 for (
size_t i=0; i<box.
size(); ++i) {
837 auto new_cons = createObjValCons(orig_vars[i],
840 interval.second - cmd_line_args_.
getDelta());
842 obj_val_cons.push_back(new_cons);
845 std::unique_ptr<Polyscip> sub_poly(
new Polyscip(cmd_line_args_,
849 sub_poly->computeNondomPoints();
855 for (
auto cons : obj_val_cons) {
861 assert (!sub_poly->unboundedResultsExist());
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));
878 return new_nondom_res;
886 vector<RectangularBox> Polyscip::computeDisjointBoxes(list<RectangularBox>&& feasible_boxes)
const {
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)) {
894 if (current->isSupersetOf(*it)) {
895 it = feasible_boxes.erase(it);
898 else if (current->isSubsetOf(*it)) {
899 current = feasible_boxes.erase(current);
900 increment_current =
false;
906 if (increment_current)
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();
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);
921 else if (box_to_be_added.isSupersetOf(elem)) {
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));
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));
933 return disjoint_boxes;
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) {
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());
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) {
959 {max(nd_01[1], nd_12[1]), nd_02[1]},
960 {max(nd_02[2], nd_12[2]), nd_01[2]}});
962 feasible_boxes.push_back(box);
967 return feasible_boxes;
981 const vector<SCIP_VAR *> &orig_vars,
982 const vector<ValueType> &orig_vals,
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),
990 [beta_i](
ValueType val){
return -beta_i*val;});
991 vars.push_back(new_var);
998 "new_variable_transformation_constraint",
999 global::narrow_cast<int>(vars.size()),
1004 assert (cons !=
nullptr);
1016 SCIP_CONS* Polyscip::createObjValCons(
const vector<SCIP_VAR *>& vars,
1017 const vector<ValueType>& vals,
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);
1025 "lhs <= c_i^T x <= rhs",
1026 global::narrow_cast<int>(vars.size()),
1027 non_const_vars.data(),
1028 non_const_vals.data(),
1031 assert (cons !=
nullptr);
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) );
1067 auto intermed_result = getOptimalResult();
1082 auto nondom_result = getOptimalResult();
1083 results.push_back(std::move(nondom_result));
1089 cout <<
"unexpected SCIP status in computeNondomProjResult: " +
1108 vector<OutcomeType> Polyscip::solveWeightedTchebycheff(
const vector<vector<SCIP_VAR*>>& orig_vars,
1109 const vector<vector<ValueType>>& orig_vals,
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_);
1118 setWeightedObjective(
WeightType(no_objs_, 0.));
1129 assert (z !=
nullptr);
1139 auto left_proj = nondom_projs.getLeftProj();
1140 auto right_proj = nondom_projs.getRightProj();
1142 assert (left_proj.getFirst() < right_proj.getFirst());
1143 assert (left_proj.getSecond() > last_proj.getSecond());
1146 auto obj_val_cons = vector<SCIP_CONS*>{};
1147 obj_val_cons.push_back( createObjValCons(orig_vars[obj_1],
1149 left_proj.getFirst(),
1150 right_proj.getFirst()));
1152 obj_val_cons.push_back( createObjValCons(orig_vars[obj_2],
1154 last_proj.getSecond(),
1155 left_proj.getSecond()));
1156 for (
auto c : obj_val_cons) {
1160 auto ref_point = std::make_pair(left_proj.getFirst() - 1., last_proj.getSecond() - 1.);
1163 auto beta_2 = (right_proj.getFirst() - ref_point.first) / (left_proj.getSecond() - ref_point.second);
1165 auto var_trans_cons = vector<SCIP_CONS*>{};
1166 var_trans_cons.push_back( createNewVarTransformCons(z,
1172 var_trans_cons.push_back(createNewVarTransformCons(z,
1177 for (
auto c : var_trans_cons) {
1182 std::unique_ptr<TwoDProj> new_proj(
nullptr);
1186 auto res = getOptimalResult();
1187 new_proj = global::make_unique<TwoDProj>(res.second, obj_1, obj_2);
1194 std::cerr <<
"Numerical troubles between " << left_proj <<
" and " << right_proj <<
"\n";
1195 std::cerr <<
"Continuing with next subproblem.\n";
1196 nondom_projs.update();
1199 std::cerr <<
"Unexpected SCIP status in solveWeightedTchebycheff: " +
1206 for (
auto c : var_trans_cons) {
1212 if (nondom_projs.epsilonDominates(left_proj, *new_proj) ||
1213 nondom_projs.epsilonDominates(right_proj, *new_proj)) {
1214 nondom_projs.update();
1219 SCIPdelVar(scip_, z, addressof(var_deleted));
1220 assert (var_deleted);
1222 computeNondomProjResult(obj_val_cons.front(),
1223 obj_val_cons.back(),
1224 new_proj->getFirst(),
1225 new_proj->getSecond(),
1229 auto nd_proj =
TwoDProj(bounded_.back().second, obj_1, obj_2);
1231 nondom_projs.update(std::move(nd_proj), bounded_.back());
1241 for (
auto c : obj_val_cons) {
1250 SCIPdelVar(scip_, z, addressof(var_deleted));
1251 assert (var_deleted);
1254 return nondom_projs.getNondomProjOutcomes();
1263 return polyscip_status_;
1271 return bounded_.size();
1279 return unbounded_.size();
1289 if (!with_presolving)
1292 setWeightedObjective(zero_weight);
1294 if (!with_presolving)
1297 setWeightedObjective(weight);
1299 if (with_presolving)
1300 separateINFORUNBD(weight,
false);
1302 cout <<
"INFORUNBD Status for problem with zero objective and no presolving.\n";
1307 cout <<
"UNBOUNDED Status for problem with zero objective.\n";
1326 else if (is_sub_prob_) {
1341 SCIP_RETCODE Polyscip::handleUnboundedStatus(
bool check_if_new_result) {
1353 auto result = getResult(
false);
1354 if (!check_if_new_result || outcomeIsNew(result.second,
false)) {
1355 unbounded_.push_back(std::move(result));
1359 cout <<
"not added to results.\n";
1370 Result Polyscip::getResult(
bool outcome_is_bounded,
SCIP_SOL *primal_sol) {
1376 for (
auto i=0; i<no_vars; ++i) {
1377 auto var_sol_val = outcome_is_bounded ?
SCIPgetSolVal(scip_, primal_sol, vars[i]) :
1383 for (
size_t index=0; index!=no_objs_; ++index) {
1384 var_obj_vals[index] = objs_probdata->getObjVal(vars[i], index, var_sol_val);
1386 std::transform(begin(outcome), end(outcome),
1387 begin(var_obj_vals),
1389 std::plus<ValueType>());
1393 return {sol, outcome};
1400 Result Polyscip::getOptimalResult() {
1402 assert (best_sol !=
nullptr);
1407 throw std::runtime_error(
"SCIPcreateFiniteSolCopy: return code != SCIP_OKAY.\n");
1408 if (!same_obj_val) {
1411 if (diff > 1.0e-5) {
1412 std::cerr <<
"absolute value difference after calling SCIPcreateFiniteSolCopy: " << diff <<
"\n";
1414 throw std::runtime_error(
"SCIPcreateFiniteSolCopy: unacceptable difference in objective values.");
1417 assert (finite_sol !=
nullptr);
1418 auto new_result = getResult(
true, finite_sol);
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;
1443 bool Polyscip::outcomeIsNew(
const OutcomeType& outcome,
1444 ResultContainer::const_iterator beg,
1445 ResultContainer::const_iterator last)
const {
1447 return std::none_of(beg, last, [
eps, &outcome](
const Result& res)
1449 return outcomesCoincide(outcome, res.second,
eps);
1461 assert (a.size() == b.size());
1462 return std::equal(begin(a), end(a), begin(b),
1465 return std::fabs(v-w) < epsilon;
1475 auto remaining_time = std::max(cmd_line_args_.
getTimeLimit() -
1492 assert (obj_probdata !=
nullptr);
1495 if (weight ==
WeightType(weight.size(), 0.)) {
1496 for (
auto i = 0; i < no_vars; ++i) {
1501 for (
auto i = 0; i < no_vars; ++i) {
1502 auto val = obj_probdata->getWeightedObjVal(vars[i], weight);
1506 fractpart = std::modf(val, &intpart);
1507 fractpart = round(pow(10,no_decimals)*fractpart)/pow(10,no_decimals);
1508 val = intpart + fractpart;
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_,
1529 auto untested_weight = weight_space_poly_->getUntestedWeight();
1530 SCIP_CALL(setWeightedObjective(untested_weight));
1534 scip_status = separateINFORUNBD(untested_weight);
1537 auto res = getOptimalResult();
1538 auto weighted_outcome = std::inner_product(begin(untested_weight),
1539 end(untested_weight),
1542 if (weighted_outcome + cmd_line_args_.
getEpsilon() <
1543 weight_space_poly_->getUntestedVertexWOV(untested_weight)) {
1544 weight_space_poly_->incorporateNewOutcome(scip_,
1547 bounded_.push_back(std::move(res));
1550 weight_space_poly_->incorporateKnownOutcome(untested_weight);
1555 weight_space_poly_->incorporateNewOutcome(scip_,
1557 unbounded_.back().second,
1561 SCIP_CALL(handleNonOptNonUnbdStatus(scip_status));
1575 auto new_line =
false;
1576 for (
const auto& result : bounded_) {
1578 outputOutcome(result.second, os);
1582 printSol(result.first, os);
1589 for (
const auto& result : unbounded_) {
1591 outputOutcome(result.second, os,
"Ray = ");
1595 printSol(result.first, os);
1609 void Polyscip::printSol(
const SolType& sol, ostream& os)
const {
1610 for (
const auto& elem : sol)
1611 os << elem.first <<
"=" << elem.second <<
" ";
1620 void Polyscip::outputOutcome(
const OutcomeType &outcome, std::ostream &os,
const std::string desc)
const {
1634 bool Polyscip::filenameIsOkay(
const string& filename) {
1635 std::ifstream file(filename.c_str());
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());
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);
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;
1675 assert (obj_probdata !=
nullptr);
1676 assert (checked_obj >= 1 && checked_obj < obj_probdata->getNoObjs());
1681 throw std::runtime_error(
"no SCIP_OKAY for SCIPlpiCreate\n.");
1684 auto obj = vector<SCIP_Real>(checked_obj, 1.);
1685 auto lb = vector<SCIP_Real>(checked_obj, 0.);
1687 auto no_nonzero = begin_nonzeros.at(checked_obj);
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]));
1702 for (
auto i=0; i<no_rows; ++i)
1703 lhs[i] = obj_probdata->getObjCoeff(vars[i], checked_obj);
1722 throw std::runtime_error(
"no SCIP_OKAY for SCIPlpiLoadColLP\n");
1726 throw std::runtime_error(
"no SCIP_OKAY for SCIPlpiSolvePrimal\n");
1729 is_redundant =
true;
1737 throw std::runtime_error(
"no SCIP_OKAY for SCIPlpiFree\n");
1739 return is_redundant;
1753 assert (obj_probdata !=
nullptr);
1754 no_objs_ = obj_probdata->getNoObjs();
1757 only_weight_space_phase_ =
true;
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));
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();
1772 cout << obj_ind <<
". objective is zero objective!\n";
1776 auto nonzero_inds = vector<int>(size, 0);
1777 std::transform(begin(nonzero_vars),
1779 begin(nonzero_inds),
1781 std::sort(begin(nonzero_inds), end(nonzero_inds));
1783 auto nonzero_vals = vector<SCIP_Real>(size, 0.);
1784 std::transform(begin(nonzero_inds),
1786 begin(nonzero_vals),
1787 [&](
int var_ind) {
return obj_probdata->getObjCoeff(vars[var_ind], obj_ind); });
1791 printObjective(obj_ind, nonzero_inds, nonzero_vals);
1793 obj_to_nonzero_inds.push_back(std::move(nonzero_inds));
1794 obj_to_nonzero_vals.push_back(std::move(nonzero_vals));
1796 if (obj_ind > 0 && objIsRedundant(begin_nonzeros,
1797 obj_to_nonzero_inds,
1798 obj_to_nonzero_vals,
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";
1811 obj_probdata->negateAllCoeffs();
1814 cout <<
"Objective sense: ";
1816 cout <<
"MAXIMIZE\n";
1818 cout <<
"MINIMIZE\n";
1819 cout <<
"Number of objectives: " << no_objs_ <<
"\n";
1830 size_t prefix = prob_file.find_last_of(
"/"),
1831 suffix = prob_file.find_last_of(
"."),
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";
1837 if (write_path.back() !=
'/')
1838 write_path.push_back(
'/');
1839 std::ofstream solfs(write_path + file_name);
1840 if (solfs.is_open()) {
1843 cout <<
"#Solution file " << file_name
1844 <<
" written to: " << write_path <<
"\n";
1847 std::cerr <<
"ERROR writing solution file\n.";
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) {
1863 else if (std::equal(begin(curr->second),
1866 std::less_equal<ValueType>())) {
1867 outputOutcome(it->second, std::cout,
"outcome ");
1868 outputOutcome(curr->second, std::cout,
" dominated by ");
1881 for (
auto cur=begin(bounded_); cur!=end(bounded_); ++cur) {
1882 if (isDominatedOrEqual(cur, begin(bounded_), end(bounded_)))
bool outputOutcomes() const
SCIP_RETCODE SCIPlpiFree(SCIP_LPI **lpi)
NondomProjections(double epsilon, const ResultContainer &supported, std::size_t first, std::size_t second)
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
bool epsilonDominates(const TwoDProj &lhs, const TwoDProj &rhs) const
std::pair< ValueType, ValueType > Interval
Interval I = [a,b)
SCIP_Real SCIPgetPrimalbound(SCIP *scip)
Class storing multiple objectives of given problem instance.
bool dominatedPointsFound() const
SCIP_RETCODE SCIPlpiSolvePrimal(SCIP_LPI *lpi)
SCIP_RETCODE SCIPdelCons(SCIP *scip, SCIP_CONS *cons)
int SCIPgetNOrigVars(SCIP *scip)
SCIP_Real ValueType
Type for computed values.
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
bool hasParameterFile() const
Class for .mop file reader.
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
SCIP_RETCODE SCIPstopClock(SCIP *scip, SCIP_CLOCK *clck)
bool isFeasible(double epsilon) const
std::vector< OutcomeType > getNondomProjOutcomes() const
SCIP_Real SCIPinfinity(SCIP *scip)
enum SCIP_Retcode SCIP_RETCODE
Status if an error occured.
SCIP_RETCODE SCIPsetPresolving(SCIP *scip, SCIP_PARAMSETTING paramsetting, SCIP_Bool quiet)
int SCIPvarGetProbindex(SCIP_VAR *var)
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
SCIP_VAR ** SCIPgetOrigVars(SCIP *scip)
ValueType getFirst() const
General types used for PolySCIP.
SCIP_RETCODE SCIPfreeClock(SCIP *scip, SCIP_CLOCK **clck)
Target narrow_cast(Source v)
doubledescription::DoubleDescriptionMethod DDMethod
abbreviation
std::vector< ValueType > OutcomeType
Type for points, rays in objective space.
SCIP_RETCODE SCIPcreate(SCIP **scip)
SCIP_Bool SCIPisTransformed(SCIP *scip)
SCIP_RETCODE SCIPsetRealParam(SCIP *scip, const char *name, SCIP_Real value)
SCIP_RETCODE SCIPlpiCreate(SCIP_LPI **lpi, SCIP_MESSAGEHDLR *messagehdlr, const char *name, SCIP_OBJSEN objsen)
std::pair< SolType, OutcomeType > Result
A result comprises of a solution/ray in feasible space and a corresponding outcome in objective space...
std::size_t numberOfBoundedResults() const
std::size_t numberofUnboundedResults() const
Status while results of weight space polyhedron are computed.
SCIP_RETCODE SCIPdelVar(SCIP *scip, SCIP_VAR *var, SCIP_Bool *deleted)
C++ wrapper for default SCIP plugins.
SCIP_RETCODE SCIPsetObjsense(SCIP *scip, SCIP_OBJSENSE objsense)
SCIP_RETCODE SCIPsolve(SCIP *scip)
SCIP_RETCODE SCIPcreateClock(SCIP *scip, SCIP_CLOCK **clck)
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
bool isSupersetOf(const RectangularBox &other) const
Status after lexicographic optimal results were computed.
SCIP_RETCODE SCIPreadProb(SCIP *scip, const char *filename, const char *extension)
Status if problem was solved successfully.
RectangularBox(const std::vector< Interval > &box)
SCIP_STATUS SCIPgetStatus(SCIP *scip)
SCIP_Bool SCIPlpiIsPrimalInfeasible(SCIP_LPI *lpi)
SCIP_RETCODE readProblem()
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)
bool isSubsetOf(const RectangularBox &other) const
Class representing the 1-skeleton of the weight space polyhedron.
friend std::ostream & operator<<(std::ostream &os, const RectangularBox &box)
TimeLimitType getTimeLimit() const
friend std::ostream & operator<<(std::ostream &os, const TwoDProj &proj)
SCIP_Bool SCIPhasPrimalRay(SCIP *scip)
SCIP_RETCODE SCIPchgVarObj(SCIP *scip, SCIP_VAR *var, SCIP_Real newobj)
bool hasTimeLimit() const
SCIP_RETCODE SCIPfreeTransform(SCIP *scip)
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
Command line arguments for PolySCIP.
SCIP_Real SCIPlpiInfinity(SCIP_LPI *lpi)
std::vector< Result > ResultContainer
Container type for results.
enum SCIP_Status SCIP_STATUS
Status after problem instance was read successfully.
SCIP_Real SCIPgetClockTime(SCIP *scip, SCIP_CLOCK *clck)
double getEpsilon() const
Status while computing 2-dimensional non-dominated projection results.
Interval getInterval(std::size_t index) const
SCIP_RETCODE SCIPfreeSol(SCIP *scip, SCIP_SOL **sol)
SCIP_RETCODE SCIPincludeObjReader(SCIP *scip, scip::ObjReader *objreader, SCIP_Bool deleteobject)
int roundWeightedObjCoeff() const
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
void printResults(std::ostream &os=std::cout) const
SCIP_RETCODE SCIPchgRhsLinear(SCIP *scip, SCIP_CONS *cons, SCIP_Real rhs)
Class representing a two-dimensional projection of an outcome.
Initial status after calling public constructor.
int SCIPgetNOrigContVars(SCIP *scip)
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...
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Status if given time limit was reached.
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
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)
SCIP_RETCODE computeNondomPoints()
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
SCIP_RETCODE SCIPreadParams(SCIP *scip, const char *filename)
std::string getProblemFile() const
std::pair< std::size_t, std::size_t > ObjPair
Pair of objectives indices.
Class representing non-dominated projections.
SCIP_Bool SCIPlpiIsPrimalFeasible(SCIP_LPI *lpi)
bool isDominated(const OutcomeType &outcome) const
ValueType getSecond() const
bool isDisjointFrom(const RectangularBox &other) const
TwoDProj getLastProj() const
Polyscip(int argc, const char *const *argv)
A rectangular box R = [a_1,e_1) x ... x [a_k,e_k) is a k-ary Cartesian product of half-open intervals...
SCIP_OBJSENSE SCIPgetObjsense(SCIP *scip)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
std::string getParameterFile() const
PolyscipStatus getStatus() const
SCIP_RETCODE SCIPchgLhsLinear(SCIP *scip, SCIP_CONS *cons, SCIP_Real lhs)
void printStatus(std::ostream &os=std::cout) const
bool onlyExtremal() const
std::vector< std::pair< std::string, ValueType > > SolType
Type for solutions in feasible space.
#define SCIP_CALL_ABORT(x)
std::string getWritePath() const
SCIP_RETCODE SCIPstartClock(SCIP *scip, SCIP_CLOCK *clck)
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
PolySCIP command line arguments.
Algorithm for transforming h-representation to v-representation.
SCIP_Real SCIPgetPrimalRayVal(SCIP *scip, SCIP_VAR *var)
std::vector< RectangularBox > getDisjointPartsFrom(double delta, const RectangularBox &other) const
SCIP_RETCODE SCIPfree(SCIP **scip)
SCIP_RETCODE SCIPcreateFiniteSolCopy(SCIP *scip, SCIP_SOL **sol, SCIP_SOL *sourcesol, SCIP_Bool *success)
TwoDProj(const OutcomeType &outcome, std::size_t first, std::size_t second)