diff --git a/include/mp/flat/constr_algebraic.h b/include/mp/flat/constr_algebraic.h index 0aba91ef2..34aa3564c 100644 --- a/include/mp/flat/constr_algebraic.h +++ b/include/mp/flat/constr_algebraic.h @@ -75,16 +75,23 @@ class AlgebraicConstraint : /// report violation amount /// (negative if holds with slack.) /// In logical mode, report 1/0 - /// (what's the violation amount of 0 for >0?) + /// (what's the violation amount + /// for "x>0" when x==0?) template - double + Violation ComputeViolation(const VarInfo& x, bool logical=false) const { double bd = Body::ComputeValue(x); if (!logical) { - return std::max( // same for strict cmp? - RhsOrRange::lb() - bd, bd - RhsOrRange::ub()); + if (RhsOrRange::lb() > bd) + return {RhsOrRange::lb() - bd, RhsOrRange::lb()}; + if (bd > RhsOrRange::ub()) + return {bd - RhsOrRange::ub(), RhsOrRange::ub()}; + return + {std::max( // negative. Same for strict cmp? + RhsOrRange::lb() - bd, bd - RhsOrRange::ub()), + 0.0}; } - return !RhsOrRange::is_valid(bd); + return {double(!RhsOrRange::is_valid(bd)), 1.0}; } /// Sorting and merging terms, some solvers require diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index f9e666470..af94db1ce 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -16,6 +16,25 @@ namespace mp { +/// Constraint/obj violation +struct Violation { + double viol_; // abs violation: >0 if really violated + double valX_; // value compared to + /// Compute whether violated + relative violation. + /// Both absolute and relative should be violated + /// (relative only if refvar!=0.) + std::pair Check( + double epsabs, double epsrel) const { + double violRel {0.0}; + if (viol_ > epsabs + && (0.0==std::fabs(valX_) + || (violRel=std::fabs(viol_/valX_))>epsrel)) + return {true, violRel}; + return {false, 0.0}; + } +}; + + /// Custom constraints to derive from, so that overloaded default settings work class BasicConstraint { public: @@ -40,8 +59,8 @@ class BasicConstraint { int GetResultVar() const { return -1; } /// Compute violation template - double ComputeViolation(const VarInfo& ) const - { return 0.0; } + Violation ComputeViolation(const VarInfo& ) const + { return {0.0, 0.0}; } private: std::string name_; @@ -159,7 +178,7 @@ class CustomFunctionalConstraint : /// Compute violation template - double ComputeViolation(const VarVec& x) const; + Violation ComputeViolation(const VarVec& x) const; }; @@ -178,7 +197,7 @@ double ComputeValue(const Con& , const VarVec& ) { /// should be redefined for cones. template -double ComputeViolation( +Violation ComputeViolation( const CustomFunctionalConstraint& c, const VarVec& x) { auto resvar = c.GetResultVar(); @@ -186,18 +205,18 @@ double ComputeViolation( auto viol = x[resvar] - ComputeValue(c, x); switch (c.GetContext().GetValue()) { case Context::CTX_MIX: - return std::fabs(viol); + return {std::fabs(viol), x[resvar]}; case Context::CTX_POS: - return viol; + return {viol, x[resvar]}; case Context::CTX_NEG: - return -viol; + return {-viol, x[resvar]}; default: - return INFINITY; + return {INFINITY, 0.0}; } } return // recomputed var minus solver's - std::fabs(x[resvar] - x.raw(resvar)) - + std::max(0.0, x.bounds_viol(resvar)); + { std::fabs(x[resvar] - x.raw(resvar)) + + std::max(0.0, x.bounds_viol(resvar)), x[resvar]}; } @@ -294,19 +313,25 @@ class ConditionalConstraint : /// If the subconstr holds but should not, /// return the opposite gap. Similar to Indicator. template - double ComputeViolation(const VarVec& x) { + Violation ComputeViolation(const VarVec& x) { auto viol = GetConstraint().ComputeViolation(x); - bool ccon_valid = viol<=0.0; + bool ccon_valid = viol.viol_<=0.0; bool has_arg = x[GetResultVar()] >= 0.5; switch (this->GetContext().GetValue()) { case Context::CTX_MIX: // Viol is non-positive if holds - return has_arg == ccon_valid ? 0.0 : std::fabs(viol); + if (has_arg == ccon_valid) + return {0.0, 0.0}; + return {std::fabs(viol.viol_), viol.valX_}; case Context::CTX_POS: - return has_arg <= ccon_valid ? 0.0 : viol; + if (has_arg <= ccon_valid) + return {0.0, 0.0}; + return {viol.viol_, viol.valX_}; case Context::CTX_NEG: - return has_arg >= ccon_valid ? 0.0 : -viol; + if (has_arg >= ccon_valid) + return {0.0, 0.0}; + return {-viol.viol_, viol.valX_}; default: - return INFINITY; + return {INFINITY, 0.0}; } } }; diff --git a/include/mp/flat/constr_eval.h b/include/mp/flat/constr_eval.h index 95232b047..4524b0d3a 100644 --- a/include/mp/flat/constr_eval.h +++ b/include/mp/flat/constr_eval.h @@ -282,13 +282,13 @@ template double ComputeValue( const ConditionalConstraint& con, const VarVec& x) { auto viol = con.GetConstraint().ComputeViolation(x, true); - bool ccon_valid = viol<=0.0; + bool ccon_valid = viol.viol_<=0.0; return ccon_valid; } /// Compute violation of the QuadraticCone constraint. template -double ComputeViolation( +Violation ComputeViolation( const QuadraticConeConstraint& con, const VarVec& x) { const auto& args = con.GetArguments(); const auto& params = con.GetParameters(); @@ -296,12 +296,13 @@ double ComputeViolation( double sum = 0.0; for (auto i=args.size(); --i; ) sum += std::pow( params[i]*x[args[i]], 2.0 ); - return std::sqrt(sum) - params[0]*x[args[0]]; + return {std::sqrt(sum) - params[0]*x[args[0]], + sum}; } /// Compute violation of the RotatedQuadraticCone constraint. template -double ComputeViolation( +Violation ComputeViolation( const RotatedQuadraticConeConstraint& con, const VarVec& x) { const auto& args = con.GetArguments(); const auto& params = con.GetParameters(); @@ -309,22 +310,24 @@ double ComputeViolation( double sum = 0.0; for (auto i=args.size(); --i>1; ) sum += std::pow( params[i]*x[args[i]], 2.0 ); - return sum - - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]]; + return {sum + - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]], + sum}; } /// Compute violation of the ExponentialCone constraint. template // ax >= by exp(cz / (by)) -double ComputeViolation( // where ax, by >= 0 +Violation ComputeViolation( // where ax, by >= 0 const ExponentialConeConstraint& con, const VarVec& x) { const auto& args = con.GetArguments(); const auto& params = con.GetParameters(); auto ax = params[0]*x[args[0]]; auto by = params[1]*x[args[1]]; if (0.0==std::fabs(by)) - return -ax; + return {-ax, 0.0}; auto cz = params[2]*x[args[2]]; - return by * std::exp(cz / by) - ax; + return {by * std::exp(cz / by) - ax, + by * std::exp(cz / by)}; } /// Compute result of the PL constraint. @@ -354,7 +357,7 @@ double ComputeValue(const PLConstraint& con, const VarVec& x) { template template -double +Violation CustomFunctionalConstraint ::ComputeViolation(const VarVec& x) const { return mp::ComputeViolation(*this, x); } diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 99ca9c586..5cfde3865 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -44,12 +44,12 @@ class IndicatorConstraint: public BasicConstraint { /// Compute violation template - double ComputeViolation(const VarInfo& x) const { + Violation ComputeViolation(const VarInfo& 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; + return {0.0, 0.0}; } @@ -151,7 +151,7 @@ class SOS_1or2_Constraint: public BasicConstraint { /// Compute violation template - double ComputeViolation(const VarInfo& x) const { + Violation ComputeViolation(const VarInfo& x) const { if (1==type) return ComputeViolationSOS1(x); return ComputeViolationSOS2(x); @@ -159,18 +159,18 @@ class SOS_1or2_Constraint: public BasicConstraint { /// Compute violation template - double ComputeViolationSOS1(const VarInfo& x) const { + Violation ComputeViolationSOS1(const VarInfo& x) const { int nnz=0; for (int i=(int)v_.size(); i--; ) { if (x.is_nonzero(v_[i])) ++nnz; } - return std::max(0, nnz-1); + return {(double)std::max(0, nnz-1), 0.0}; } /// Compute violation template - double ComputeViolationSOS2(const VarInfo& x) const { + Violation ComputeViolationSOS2(const VarInfo& x) const { int npos=0, pos1=-1, pos2=-1; int posDist=1; // should be 1 for (int i=(int)v_.size(); i--; ) { @@ -184,8 +184,8 @@ class SOS_1or2_Constraint: public BasicConstraint { } } } - return std::max(0, npos-2) - + std::abs(1 - posDist); + return {double(std::max(0, npos-2) + + std::abs(1 - posDist)), 0.0}; } @@ -247,13 +247,13 @@ class ComplementarityConstraint : public BasicConstraint { /// Compute violation template - double ComputeViolation(const VarInfo& x) const { + Violation ComputeViolation(const VarInfo& x) const { auto ve = compl_expr_.ComputeValue(x); if (x.is_at_lb(compl_var_)) - return -ve; + return {-ve, 0.0}; else if (x.is_at_ub(compl_var_)) - return ve; - return std::fabs(ve); + return {ve, 0.0}; + return {std::fabs(ve), 0.0}; } diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h index 95d5a5cde..b0c3f36f4 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -32,24 +32,36 @@ 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 + /// Check if this violation should be counted. void CheckViol( - double val, double eps, const char* nm) { - if (val > eps) - CountViol(val, nm); + Violation viol, + double epsabs, double epsrel, + const char* nm) { + auto chk = viol.Check(epsabs, epsrel); + if (chk.first) + CountViol(viol, chk.second, nm); } /// Count violation - void CountViol(double val, const char* nm) { + void CountViol( + Violation viol, double violRel, const char* nm) { ++N_; - if (epsMax_ < val) - epsMax_ = val; - name_ = nm; + if (epsAbsMax_ < viol.viol_) { + epsAbsMax_ = viol.viol_; + nameAbs_ = nm; + } + if (epsRelMax_ < violRel) { + epsRelMax_ = violRel; + nameRel_ = nm; + } } int N_ {0}; - double epsMax_ {0.0}; - const char* name_ {nullptr}; + double epsAbsMax_ {0.0}; + const char* nameAbs_ {nullptr}; + double epsRelMax_ {0.0}; + const char* nameRel_ {nullptr}; }; + /// Array of violation summaries. /// For different kinds, e.g., original / aux vars. template @@ -245,13 +257,13 @@ struct SolCheck { ArrayRef x_raw, ArrayRef vtype, ArrayRef lb, ArrayRef ub, - double feastol, double , //inttol, + double feastol, double feastolrel, const char* sol_rnd, const char* sol_prec, bool recomp_vals, int chk_mode) : x_(feastol, recomp_vals, x, x_raw, vtype, lb, ub, sol_rnd, sol_prec), obj_(obj), - feastol_(feastol), + feastol_(feastol), feastolrel_(feastolrel), fRecomputedVals_(recomp_vals), check_mode_(chk_mode) { } /// Any violations? @@ -279,8 +291,10 @@ struct SolCheck { const ArrayRef& obj_vals() const { return obj_; } - /// Feasibility tolerance + /// Absolute feasibility tolerance double GetFeasTol() const { return feastol_; } + /// Relative feasibility tolerance + double GetFeasTolRel() const { return feastolrel_; } /// Using recomputed aux vars? bool if_recomputed() const { return fRecomputedVals_; } @@ -313,6 +327,7 @@ struct SolCheck { VarInfoStatic x_; ArrayRef obj_; double feastol_; + double feastolrel_; bool fRecomputedVals_; int check_mode_; @@ -906,7 +921,9 @@ class ConstraintKeeper final c_class = 4; // intermediate if (c_class & chk.check_mode()) { auto viol = cons_[i].con_.ComputeViolation(x); - if (viol > chk.GetFeasTol()) { + auto cr = viol.Check( + chk.GetFeasTol(), chk.GetFeasTolRel()); + if (cr.first) { if (!conviolarray) conviolarray = // lazy map access &conviolmap[GetShortTypeName()]; @@ -917,7 +934,7 @@ class ConstraintKeeper final ? 2 : 1; assert(index < (int)conviolarray->size()); (*conviolarray)[index].CountViol( - viol, cons_[i].con_.name()); + viol, cr.second, cons_[i].con_.name()); } } } diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h index 1e459855d..59e2bad54 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -617,8 +617,8 @@ class FlatConverter : ArrayRef obj) { bool result = true; std::string err_msg; - try { // protect - std::vector x_back; // to extract x used + try { // protect + std::vector x_back = x; if (options_.solcheckmode_ & (1+2+4+8+16)) { if (!DoCheckSol(x, duals, obj, {}, x_back, false)) result = false; @@ -680,7 +680,11 @@ class FlatConverter : return std::move(vir.get_x().get_x()); } - /// Check single unpostsolved solution + /// Check single unpostsolved solution. + /// @param x_back: solution vector from realistic mode. + /// It can be changed by AMPL solution_... options. + /// Its auxiliary vars are compared + /// with recomputed expression values. bool DoCheckSol( ArrayRef x, const pre::ValueMapDbl& duals, @@ -693,7 +697,7 @@ class FlatConverter : GetModel().var_lb_vec(), GetModel().var_ub_vec(), options_.solfeastol_, - options_.solinttol_, + options_.solfeastolrel_, options_.dont_use_sol_round_ ? "" : std::getenv("solution_round"), options_.dont_use_sol_prec_ @@ -709,6 +713,9 @@ class FlatConverter : if (chk.check_mode() & 16) CheckObjs(chk); GenerateViolationsReport(chk); + // Should messages for realistic and idealistic + // modes be coordinated? + // I.e., when no expressions. // What if this is an intermediate solution? // Should be fine - warning by default, // fail if requested explicitly. @@ -726,7 +733,7 @@ class FlatConverter : chk.GetReport(), true); // replace for multiple solutions } - x_back = chk.x_ext().get_x(); + x_back = chk.x_ext().get_x(); // to reuse 'realistic' vector return !chk.HasAnyViols(); } @@ -737,17 +744,18 @@ class FlatConverter : // no aux vars for idealistic mode if (!aux || !chk.if_recomputed()) { chk.VarViolBnds().at(aux).CheckViol( - MPCD( lb(i) ) - x, - options_.solfeastol_, + {MPCD( lb(i) ) - x, MPCD( lb(i) )}, + options_.solfeastol_, options_.solfeastolrel_, GetModel().var_name(i)); chk.VarViolBnds().at(aux).CheckViol( - x - MPCD( ub(i) ), - options_.solfeastol_, + {x - MPCD( ub(i) ), MPCD( ub(i) )}, + options_.solfeastol_, options_.solfeastolrel_, 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)), + std::round(x) }, + options_.solinttol_, INFINITY, GetModel().var_name(i)); } } @@ -765,10 +773,10 @@ class FlatConverter : for (auto i =std::min(objs.size(), chk.obj_vals().size()); i--; ) { + auto val1 = ComputeValue(objs[i], chk.x_ext()); chk.ObjViols().CheckViol( - std::fabs(chk.obj_vals()[i] - - ComputeValue(objs[i], chk.x_ext())), - options_.solfeastol_, + {std::fabs(chk.obj_vals()[i] - val1), val1}, + options_.solfeastol_, options_.solfeastolrel_, objs[i].name()); } } @@ -777,9 +785,10 @@ class FlatConverter : fmt::MemoryWriter wrt; if (chk.HasAnyViols()) wrt.write( - " [ sol:chk:feastol={}, sol:chk:inttol={},\n" + " [ sol:chk:feastol={}, :feastolrel={}, :inttol={},\n" " solution_round='{}', solution_precision='{}' ]\n", - options_.solfeastol_, options_.solinttol_, + options_.solfeastol_, options_.solfeastolrel_, + options_.solinttol_, chk.x_ext().solution_round(), chk.x_ext().solution_precision()); if (chk.HasAnyConViols()) { @@ -802,22 +811,42 @@ class FlatConverter : chk.SetReport( wrt.str() ); } + /// Generate message about 1 violation. + /// @param f_max: whether we need to print + /// the maximal violations. void Gen1Viol( const ViolSummary& vs, fmt::MemoryWriter& wrt, bool f_max, const std::string& format) { if (vs.N_) { wrt.write(format, vs.N_); - if (f_max) - wrt.write(",\n up to {:.0E}", vs.epsMax_); - if (vs.name_ && *vs.name_ != '\0') { - if (!f_max) - wrt.write("\n "); - wrt.write(" (item '{}')", vs.name_); - } + auto vmaxabs = Gen1ViolMax( + f_max, vs.epsAbsMax_, vs.nameAbs_, false); + auto vmaxrel = Gen1ViolMax( + f_max, vs.epsRelMax_, vs.nameRel_, true); + if (vmaxabs.size() || vmaxrel.size()) + wrt.write(",\n {}", vmaxabs); + if (vmaxabs.size() && vmaxrel.size()) + wrt.write(", "); + wrt.write("{}", vmaxrel); wrt.write("\n"); } } + /// Stringify 1 maximal violation + std::string Gen1ViolMax( + bool f_max, double viol, const char* nm, bool relVsAbs) { + fmt::MemoryWriter wrt; + if (f_max) + wrt.write("up to {:.0E} ({}", + viol, relVsAbs ? "rel" : "abs"); + if (nm && *nm != '\0') { + wrt.write(f_max ? ", " : "("); + wrt.write("item '{}'", nm); + } else if (f_max) + wrt.write(")"); + return wrt.str(); + } + void GenConViol( const std::map< std::string, ViolSummArray<3> >& cvmap, fmt::MemoryWriter& wrt, int alg_log) { @@ -1146,6 +1175,7 @@ class FlatConverter : int solcheckmode_ = 1+2+16+512; bool solcheckfail_ = false; double solfeastol_ = 1e-6; + double solfeastolrel_ = 1e-6; double solinttol_ = 1e-5; bool dont_use_sol_round_ = false; bool dont_use_sol_prec_ = false; @@ -1246,16 +1276,20 @@ class FlatConverter : "\n" "Default: 1+2+16+512.", options_.solcheckmode_, 0, 1024); - GetEnv().AddOption("sol:chk:feastol sol:chk:eps sol:eps chk:eps", - "Solution checking tolerance for objective values, variable " + GetEnv().AddOption("sol:chk:feastol sol:chk:eps chk:eps chk:tol", + "Absolute tolerance to check objective values, variable " "and constraint bounds. Default 1e-6.", options_.solfeastol_, 0.0, 1e100); + GetEnv().AddOption("sol:chk:feastolrel sol:chk:epsrel chk:epsrel chk:tolrel", + "Relative tolerance to check objective values, variable " + "and constraint bounds. Default 1e-6.", + options_.solfeastolrel_, 0.0, 1e100); GetEnv().AddOption("sol:chk:inttol sol:chk:inteps sol:inteps chk:inteps", "Solution checking tolerance for variables' integrality. " "Default 1e-5.", options_.solinttol_, 0.0, 1e100); GetEnv().AddOption("sol:chk:fail chk:fail checkfail", - "Fail on constraint violations.", + "Fail on solution checking violations.", options_.solcheckfail_, false, true); GetEnv().AddOption("sol:chk:noround chk:noround chk:no_solution_round", "Don't use AMPL solution_round option when checking.",