From 576cf27cc7a673ea59577307a914202a6c6d25f3 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Tue, 22 Aug 2023 14:17:10 +1000 Subject: [PATCH] Check algebraic and indicator constraints #200 --- include/mp/backend-app.h | 3 +- include/mp/flat/constr_algebraic.h | 8 ++ include/mp/flat/constr_base.h | 3 +- include/mp/flat/constr_general.h | 16 +++- include/mp/flat/constr_keeper.h | 54 +++++++++---- include/mp/flat/converter.h | 80 +++++++++++++++++-- include/mp/flat/converter_model.h | 6 ++ include/mp/flat/expr_affine.h | 6 +- include/mp/flat/expr_quadratic.h | 8 +- include/mp/valcvt.h | 9 ++- .../categorized/fast/tech/infeas_int_01.mod | 21 +++++ .../categorized/fast/tech/modellist.json | 12 +++ 12 files changed, 186 insertions(+), 40 deletions(-) create mode 100644 test/end2end/cases/categorized/fast/tech/infeas_int_01.mod diff --git a/include/mp/backend-app.h b/include/mp/backend-app.h index 0409b8ec3..3e17abcca 100644 --- a/include/mp/backend-app.h +++ b/include/mp/backend-app.h @@ -88,7 +88,8 @@ int BackendApp::Run(char **argv) { // if the solution handler is available. GetBackend().ReportError( er.exit_code()>=0 ? er.exit_code() : 500, - er.what()); + std::string(GetBackend().long_name()) + ": " + + er.what()); } catch (const std::exception& ex) { // For std::exception, which can be thrown by anything, // we try to print the result into .sol file, diff --git a/include/mp/flat/constr_algebraic.h b/include/mp/flat/constr_algebraic.h index 4893447a6..c45882f75 100644 --- a/include/mp/flat/constr_algebraic.h +++ b/include/mp/flat/constr_algebraic.h @@ -70,6 +70,14 @@ class AlgebraicConstraint : return Body::ComputeValue(x) - RhsOrRange::lb(); } + /// Compute violation + double + ComputeViolation(const ArrayRef& x) const { + auto bd = Body::ComputeValue(x); + return std::max( + RhsOrRange::lb() - bd, bd - RhsOrRange::ub()); + } + /// Sorting and merging terms, some solvers require void sort_terms() { Body::sort_terms(); } diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index fc22850d6..42caf8897 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -36,7 +36,8 @@ class BasicConstraint { /// For functional constraints, result variable index int GetResultVar() const { return -1; } /// Compute violation - double ComputeViolation(const ArrayRef& ) { return 0.0; } + double ComputeViolation(const ArrayRef& ) const + { return 0.0; } private: std::string name_; diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 9e16a9361..49bab9786 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -38,13 +38,23 @@ class IndicatorConstraint: public BasicConstraint { /// Getters int get_binary_var() const { return b_; } int get_binary_value() const { return bv_; } - bool is_binary_value_1() const { return 1==get_binary_value(); } + bool is_binary_value_1() const + { return 1==get_binary_value(); } const Con& get_constraint() const { return con_; } + /// Compute violation + double ComputeViolation(const ArrayRef& x) const { + assert(b_<(int)x.size()); + auto xb = x[b_]; + if (std::round(xb) == bv_) // Implication needed + return con_.ComputeViolation(x); + return 0.0; + } + private: - const int b_=-1; // the indicator variable - const int bv_=1; // the value, 0/1 + const int b_=-1; // the indicator variable + const int bv_=1; // the value, 0/1 const Con con_; }; diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h index 99095e773..eda4ceed6 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -30,21 +30,21 @@ static const mp::OptionValueInfo values_item_acceptance[] = { /// Violation summary for a class of vars/cons/objs struct ViolSummary { /// Check if this violation should be counted - void CheckViol(double val, double eps) { - if (val > eps) { - ++N_; - if (epsMax_ < val) - epsMax_ = val; - } + void CheckViol( + double val, double eps, const char* nm) { + if (val > eps) + CountViol(val, nm); } /// Count violation - void CountViol(double val) { + void CountViol(double val, const char* nm) { ++N_; if (epsMax_ < val) epsMax_ = val; + name_ = nm; } int N_ {0}; double epsMax_ {0.0}; + const char* name_ {nullptr}; }; /// Array of violation summaries. @@ -62,7 +62,19 @@ struct SolCheck { : x_(x), y_(duals), obj_(obj), feastol_(feastol), inttol_(inttol) { } /// Any violations? - bool HasAnyViols() const { return hasAnyViol_; } + bool HasAnyViols() const + { return HasAnyConViols() || HasAnyObjViols(); } + /// Any constraint violations? + bool HasAnyConViols() const { + return viol_var_bnds_[0].N_ || viol_var_bnds_[1].N_ + || viol_var_int_[0].N_ || viol_var_int_[1].N_ + || viol_cons_alg_.size() + || viol_cons_log_.size(); + } + /// Any objective value violations? + bool HasAnyObjViols() const + { return viol_obj_.N_; } + /// Summary const std::string& GetReport() const { return report_; } @@ -88,6 +100,13 @@ struct SolCheck { std::map< std::string, ViolSummArray<3> >& ConViolLog() { return viol_cons_log_; } + /// Obj viols + ViolSummary& ObjViols() { return viol_obj_; } + + /// Set report + void SetReport(std::string rep) + { report_ = std::move(rep); } + private: ArrayRef x_; const pre::ValueMapDbl& y_; @@ -95,7 +114,6 @@ struct SolCheck { double feastol_; double inttol_; - bool hasAnyViol_ = false; std::string report_; /// Variable bounds: orig, aux @@ -595,23 +613,27 @@ class ConstraintKeeper : public BasicConstraintKeeper { } /// Compute violations for this constraint type. - /// We do it for redefined ones too. + /// We do it for redefined (intermediate) ones too. void ComputeViolations(SolCheck& chk) { if (cons_.size()) { auto& conviolmap = cons_.front().con_.IsLogical() ? - chk.ConViolAlg() : - chk.ConViolLog(); - auto& conviolarray = - conviolmap[cons_.front().con_.GetTypeName()]; + chk.ConViolLog() : + chk.ConViolAlg(); const auto& x = chk.x(); + ViolSummArray<3>* conviolarray {nullptr}; for (int i=(int)cons_.size(); i--; ) { auto viol = cons_[i].con_.ComputeViolation(x); if (viol > chk.GetFeasTol()) { + if (!conviolarray) + conviolarray = + &conviolmap[cons_.front().con_.GetTypeName()]; /// Solver-side? /// TODO also original NL constraints (index 0) - int index = cons_[i].IsDeleted() ? 2 : 1; - conviolarray[index].CountViol(viol); + int index = cons_[i].IsDeleted() ? 1 : 2; + assert(index < (int)conviolarray->size()); + (*conviolarray)[index].CountViol( + viol, cons_[i].con_.name()); } } } diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h index 5d8a0f91b..57e81776b 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -622,7 +622,8 @@ class FlatConverter : // For now, do this via warnings? if (chk.HasAnyViols()) { if (options_.solcheckfail_) - MP_RAISE(chk.GetReport()); + MP_RAISE_WITH_CODE(520, // numeric error + chk.GetReport()); else AddWarning("SolutionCheck", chk.GetReport()); } @@ -634,13 +635,18 @@ class FlatConverter : auto x = chk.x(i); bool aux = !MPCD( is_var_original(i) ); chk.VarViolBnds().at(aux).CheckViol( - MPCD( lb(i) ) - x, options_.solfeastol_); + MPCD( lb(i) ) - x, + options_.solfeastol_, + GetModel().var_name(i)); chk.VarViolBnds().at(aux).CheckViol( - x - MPCD( ub(i) ), options_.solfeastol_); + x - MPCD( ub(i) ), + options_.solfeastol_, + GetModel().var_name(i)); if (is_var_integer(i)) chk.VarViolIntty().at(aux).CheckViol( - std::fabs(x - std::round(x)), - options_.solinttol_); + std::fabs(x - std::round(x)), + options_.solinttol_, + GetModel().var_name(i)); } } @@ -649,10 +655,68 @@ class FlatConverter : GetModel().ComputeViolations(chk); } - void CheckObjs(SolCheck& ) { + void CheckObjs(SolCheck& ) { } + + void GenerateViolationsReport(SolCheck& chk) { + fmt::MemoryWriter wrt; + if (chk.HasAnyConViols()) { + wrt.write( + "Constraint violations " + "(sol:chk:feastol={}, sol:chk:inttol={}):\n", + options_.solfeastol_, options_.solinttol_); + Gen1Viol(chk.VarViolBnds().at(0), wrt, + " - {} original variable(s) violate bounds, max by {}"); + Gen1Viol(chk.VarViolBnds().at(1), wrt, + " - {} auxiliary variable(s) violate bounds, max by {}"); + Gen1Viol(chk.VarViolIntty().at(0), wrt, + " - {} original variable(s) violate integrality, max by {}"); + Gen1Viol(chk.VarViolIntty().at(1), wrt, + " - {} auxiliary variable(s) violate integrality, max by {}"); + } + GenConViol(chk.ConViolAlg(), wrt, "Algebraic"); + GenConViol(chk.ConViolLog(), wrt, "Logical"); + if (chk.HasAnyObjViols()) { + wrt.write("Objective value violations" + "(sol:chk:feastol={})\n", + options_.solfeastol_); + Gen1Viol(chk.ObjViols(), wrt, + " - {} objective value(s) violated, max by {}"); + } + chk.SetReport( wrt.str() ); + } + + void Gen1Viol( + const ViolSummary& vs, fmt::MemoryWriter& wrt, + const std::string& format) { + if (vs.N_) { + wrt.write(format, vs.N_, vs.epsMax_); + if (vs.name_ && *vs.name_ != '\0') + wrt.write(" (item '{}')", vs.name_); + wrt.write("\n"); + } } - void GenerateViolationsReport(SolCheck& chk) { } + void GenConViol( + const std::map< std::string, ViolSummArray<3> >& cvmap, + fmt::MemoryWriter& wrt, const std::string& classnm) { + if (cvmap.size()) { + wrt.write(classnm + " constraint violations:\n"); + for (const auto& cva: cvmap) { + Gen1Viol(cva.second.at(0), wrt, + " - {} original constraint(s) of type '" + + std::string(cva.first) + + "' violate bounds, max by {}"); + Gen1Viol(cva.second.at(1), wrt, + " - {} reformulated constraint(s) of type '" + + std::string(cva.first) + + "' violate bounds, max by {}"); + Gen1Viol(cva.second.at(2), wrt, + " - {} auxiliary constraint(s) of type '" + + std::string(cva.first) + + "' violate bounds, max by {}"); + } + } + } //////////////////////////// UTILITIES ///////////////////////////////// /// @@ -961,7 +1025,7 @@ class FlatConverter : bool solcheckfail_ = false; double solfeastol_ = 1e-6; - double solinttol_ = 1e-6; + double solinttol_ = 1e-5; }; Options options_; diff --git a/include/mp/flat/converter_model.h b/include/mp/flat/converter_model.h index 39e25972e..e5265e83e 100644 --- a/include/mp/flat/converter_model.h +++ b/include/mp/flat/converter_model.h @@ -61,6 +61,12 @@ class FlatModel var_names_storage_ = names; } + /// To be used only after names presolve + const char* var_name(int i) const { + return i<(int)var_names_storage_.size() + ? var_names_storage_[i].c_str() : nullptr; + } + int num_vars() const { assert(check_vars()); return (int)var_lb_.size(); } diff --git a/include/mp/flat/expr_affine.h b/include/mp/flat/expr_affine.h index db2a1d29b..980c4e2af 100644 --- a/include/mp/flat/expr_affine.h +++ b/include/mp/flat/expr_affine.h @@ -74,10 +74,10 @@ class LinTerms { const LinTerms& GetLinTerms() const { return *this; } /// Compute value given a dense vector of variable values - double ComputeValue(ArrayRef x) const { - double s=0.0; + long double ComputeValue(const ArrayRef& x) const { + long double s=0.0; for (size_t i=coefs().size(); i--; ) - s += coefs()[i] * x[vars()[i]]; + s += (long double)(coefs()[i]) * x[vars()[i]]; return s; } diff --git a/include/mp/flat/expr_quadratic.h b/include/mp/flat/expr_quadratic.h index 6ca919e75..bba4b4e7d 100644 --- a/include/mp/flat/expr_quadratic.h +++ b/include/mp/flat/expr_quadratic.h @@ -54,10 +54,10 @@ class QuadTerms { int var2(int i) const { return vars2_[i]; } /// Compute value given a dense vector of variable values - double ComputeValue(ArrayRef x) const { - double s=0.0; + long double ComputeValue(const ArrayRef& x) const { + long double s=0.0; for (size_t i=coefs().size(); i--; ) - s += coefs()[i] * x[vars1()[i]] * x[vars2()[i]]; + s += (long double)(coefs()[i]) * x[vars1()[i]] * x[vars2()[i]]; return s; } @@ -200,7 +200,7 @@ class QuadAndLinTerms : } /// Value at given variable vector - double ComputeValue(ArrayRef x) const { + long double ComputeValue(const ArrayRef& x) const { return LinTerms::ComputeValue(x) + QuadTerms::ComputeValue(x); } diff --git a/include/mp/valcvt.h b/include/mp/valcvt.h index 9bb65757c..2a69b8330 100644 --- a/include/mp/valcvt.h +++ b/include/mp/valcvt.h @@ -296,10 +296,11 @@ class ValuePresolver : public ValuePresolverImpl { const auto& mx = mv.GetVarValues(); const auto& mo = mv.GetObjValues(); if (mx.IsSingleKey() - && mx().size() // solution available - && mo.IsSingleKey()) { - solchk_(mx(), - mv.GetConValues(), mo()); + && mx().size()) { // solution available + ArrayRef objs; + if (mo.IsSingleKey()) + objs = mo(); + solchk_(mx(), mv.GetConValues(), objs); } } return ValuePresolverImpl::PostsolveSolution(mv); diff --git a/test/end2end/cases/categorized/fast/tech/infeas_int_01.mod b/test/end2end/cases/categorized/fast/tech/infeas_int_01.mod new file mode 100644 index 000000000..b70feeb13 --- /dev/null +++ b/test/end2end/cases/categorized/fast/tech/infeas_int_01.mod @@ -0,0 +1,21 @@ + +# ------------------------------------------------------------- +# IIS Int 01 +# infeas_int_01.mod +# ------------------------------------------------------------- + +var x integer; +var y integer; + +minimize TotalSum: + x - 2*y; + +subj to C1: + -x + 21*y >= 2; + +subj to C2: + -3*x + 2*y <= 1; + +subj to C3: + 20*x + y <= 20; + diff --git a/test/end2end/cases/categorized/fast/tech/modellist.json b/test/end2end/cases/categorized/fast/tech/modellist.json index 9698b5bab..8d4bdad88 100644 --- a/test/end2end/cases/categorized/fast/tech/modellist.json +++ b/test/end2end/cases/categorized/fast/tech/modellist.json @@ -25,5 +25,17 @@ "files" : ["diet.mod", "diet.dat"], "options": { "ANYSOLVER_options": "outlev 1 barrier outlev 0" }, "objective": 88.2 + }, + { + "name" : "Check sol:chk:fail", + "tags" : ["linear", "feasrelax"], + "files" : ["infeas_int_01.mod"], + "options": { + "ANYSOLVER_options": + "sol:chk:fail sol:chk:feastol=1e-7 sol:chk:inttol=1e-4 feasrelax=1" + }, + "values": { + "solve_result_num": 520 + } } ]