From 0a0e96f63f6371b256a1fbd3424cf05749a537c3 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Wed, 30 Aug 2023 02:33:38 +1000 Subject: [PATCH] SolCheck with recomputed aux vars #200 This comes close toAMPL's slacks computation --- include/mp/flat/constr_algebraic.h | 5 +- include/mp/flat/constr_base.h | 14 ++-- include/mp/flat/constr_eval.h | 29 ++++--- include/mp/flat/constr_general.h | 3 +- include/mp/flat/constr_keeper.h | 120 ++++++++++++++++++++++++----- include/mp/flat/converter.h | 46 ++++++++++- include/mp/flat/expr_affine.h | 3 +- include/mp/flat/expr_algebraic.h | 3 +- include/mp/flat/expr_quadratic.h | 6 +- 9 files changed, 181 insertions(+), 48 deletions(-) diff --git a/include/mp/flat/constr_algebraic.h b/include/mp/flat/constr_algebraic.h index 337ba6d6b..f8ac74039 100644 --- a/include/mp/flat/constr_algebraic.h +++ b/include/mp/flat/constr_algebraic.h @@ -71,9 +71,10 @@ class AlgebraicConstraint : } /// Compute violation. - /// Negative if if holds with slack. + /// Negative if holds with slack. + template double - ComputeViolation(const ArrayRef& x) const { + ComputeViolation(const VarInfo& x) const { auto bd = Body::ComputeValue(x); return std::max( // same for strict cmp? RhsOrRange::lb() - bd, bd - RhsOrRange::ub()); diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index ed4bfb006..7374383be 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -39,7 +39,8 @@ class BasicConstraint { /// For functional constraints, result variable index int GetResultVar() const { return -1; } /// Compute violation - double ComputeViolation(const ArrayRef& ) const + template + double ComputeViolation(const VarInfo& ) const { return 0.0; } private: @@ -165,15 +166,10 @@ class CustomFunctionalConstraint : /// Dummy template: compute result of functional constraint. /// Should be implemented for proper functional specializations, /// but has no sense for conic constraints, for example. -template -double ComputeValue( - const CustomFunctionalConstraint& , - const VarVec& ) { +template +double ComputeValue(const Con& , const VarVec& ) { MP_RAISE(fmt::format("ComputeValue({}) not implemented.", - typeid( - CustomFunctionalConstraint).name())); + typeid(Con).name())); return 0.0; } diff --git a/include/mp/flat/constr_eval.h b/include/mp/flat/constr_eval.h index 6e83a2996..b5526fed4 100644 --- a/include/mp/flat/constr_eval.h +++ b/include/mp/flat/constr_eval.h @@ -254,24 +254,31 @@ double ComputeValue(const AtanhConstraint& con, const VarVec& x) { return std::atanh(x[con.GetArguments()[0]]); } +/// Compute result of the LinearFuncCon constraint. +template +double ComputeValue( + const LinearFunctionalConstraint& con, const VarVec& x) { + return con.GetAffineExpr().ComputeValue(x); +} + +/// Compute result of the QuadrFuncCon constraint. +template +double ComputeValue( + const QuadraticFunctionalConstraint& con, const VarVec& x) { + return con.GetQuadExpr().ComputeValue(x); +} + + + /// Compute result of a conditional constraint. +/// Just return bool(viol(subcon) <= 0). /// This is not used to compute violation. template double ComputeValue( const ConditionalConstraint& con, const VarVec& x) { auto viol = con.GetConstraint().ComputeViolation(x); bool ccon_valid = viol<=0.0; - bool has_arg = x[con.GetResultVar()] >= 0.5; - switch (con.GetContext().GetValue()) { - case Context::CTX_MIX: - return has_arg == ccon_valid; - case Context::CTX_POS: - return has_arg <= ccon_valid; - case Context::CTX_NEG: - return has_arg >= ccon_valid; - default: - return INFINITY; - } + return ccon_valid; } /// Compute violation of the QuadraticCone constraint. diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index ce14183a9..99ca9c586 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -43,7 +43,8 @@ class IndicatorConstraint: public BasicConstraint { const Con& get_constraint() const { return con_; } /// Compute violation - double ComputeViolation(const ArrayRef& x) const { + template + double ComputeViolation(const VarInfo& x) const { assert(b_<(int)x.size()); auto xb = x[b_]; if (std::round(xb) == bv_) // Implication needed diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h index 136c7803c..2b4de2cbc 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -55,31 +55,84 @@ struct ViolSummary { template using ViolSummArray = std::array; +class VarInfoRecomp; + +/// Function prototype to recompute +/// variable at index \a i. +using VarsRecomputeFn + = std::function< + double(int i, const VarInfoRecomp& x) >; + +/// Var vector managing recomputation +class VarVecRecomp { +public: + /// Construct + VarVecRecomp(std::vector x, + VarsRecomputeFn rec_fn) + : x_(std::move(x)), is_recomp_(x_.size()), + recomp_fn_(rec_fn) { assert(recomp_fn_); } + /// Set p_var_info_recomp_ + void set_p_var_info(const VarInfoRecomp* p) const + { p_var_info_recomp_ = p; } + /// Number of vars + int size() const { return (int)x_.size(); } + /// Access variable value. + /// Recompute if not yet. + double operator[]( int i ) const { + assert(i>=0 && i<(int)x_.size()); + assert(p_var_info_recomp_ != nullptr); + if (!is_recomp_[i]) { + x_[i] = recomp_fn_(i, *p_var_info_recomp_); + is_recomp_[i] = true; + } + return x_[i]; + } + /// Expose begin() + std::vector::iterator begin() { return x_.begin(); } + /// Expose end() + std::vector::iterator end() { return x_.end(); } + /// Move out x + std::vector& get_x() const { return x_; } + +private: + mutable std::vector x_; + mutable std::vector is_recomp_; + VarsRecomputeFn recomp_fn_; + mutable const VarInfoRecomp* p_var_info_recomp_{nullptr}; +}; + + +/// Static var vector +using VarVecStatic = std::vector; + /// Variable information used by solution check -class VarInfo { +template +class VarInfoImpl { public: /// Constructor - VarInfo(double ft, - std::vector x, + VarInfoImpl(double ft, + VarVec x, ArrayRef type, ArrayRef lb, ArrayRef ub, const char* sol_rnd, const char* sol_prec) - : feastol_(ft), x_(std::move(x)), x_ref_(x_), + : feastol_(ft), x_(std::move(x)), type_(type), lb_(lb), ub_(ub) { - assert(x_.size()>=type_.size()); // feasrelax can add more + assert((int)x_.size()>=(int)type_.size()); // feasrelax can add more assert(type_.size()==lb_.size()); assert(type_.size()==ub_.size()); - apply_precision_options(sol_rnd, sol_prec); + apply_precision_options(sol_rnd, sol_prec); // after recomp? } + /// Number of vars + int size() const { return (int)x_.size(); } /// Access variable value double operator[]( int i ) const { assert(i>=0 && i<(int)x_.size()); return x_[i]; } - /// Access whole vectorref - operator const ArrayRef& () const - { return x_ref(); } + /// Access VarVec + const VarVec& get_x() const { return x_; } + /// Access integrality condition bool is_var_int(int i) const { assert(i>=0 && i<(int)type_.size()); @@ -113,9 +166,6 @@ class VarInfo { std::string solution_precision() const { return sol_prec_ < 100 ? std::to_string(sol_prec_) : ""; } - /// x() as ArrayRef - const ArrayRef& x_ref() const { return x_ref_; } - protected: void apply_precision_options( const char* sol_rnd, const char* sol_prec) { @@ -140,8 +190,7 @@ class VarInfo { private: double feastol_; - std::vector x_; // can be rounded, etc. - const ArrayRef x_ref_; + VarVec x_; // can be rounded, recomputed, etc. const ArrayRef type_; const ArrayRef lb_; const ArrayRef ub_; @@ -150,6 +199,20 @@ class VarInfo { }; +/// VarInfoRecompTypedef +using VarInfoRecompTypedef = VarInfoImpl; + +/// Define VarInfoRecomp +class VarInfoRecomp : public VarInfoRecompTypedef { +public: + /// Inherit ctor's + using VarInfoRecompTypedef::VarInfoRecompTypedef; +}; + +/// VarInfoStatic +using VarInfoStatic = VarInfoImpl; + + /// Solution check data struct SolCheck { /// Construct @@ -183,9 +246,7 @@ struct SolCheck { const std::string& GetReport() const { return report_; } /// VarInfo, can be used like x() for templates - const VarInfo& x_ext() const { return x_; } - /// x as array - const ArrayRef& x() { return x_.x_ref(); } + const VarInfoStatic& x_ext() const { return x_; } /// x[i] double x(int i) const { return x_[i]; } /// Feasibility tolerance @@ -216,7 +277,7 @@ struct SolCheck { { report_ = std::move(rep); } private: - VarInfo x_; + VarInfoStatic x_; const pre::ValueMapDbl& y_; ArrayRef obj_; double feastol_; @@ -391,9 +452,16 @@ class BasicConstraintKeeper { /// Mark as unused. Use index only. virtual void MarkAsUnused(int i) = 0; + /// Is constraint \a i unused? + virtual bool IsUnused(int i) const = 0; + /// Copy names from ValueNodes virtual void CopyNamesFromValueNodes() = 0; + /// Compute result for constraint \a i + /// (for functional constraints) + virtual double ComputeValue(int i, const VarInfoRecomp& ) = 0; + /// Compute violations virtual void ComputeViolations(SolCheck& ) = 0; @@ -524,7 +592,8 @@ class BasicFlatConverter { /// Specialize ConstraintKeeper for a given constraint type /// to store an array of such constraints template -class ConstraintKeeper : public BasicConstraintKeeper { +class ConstraintKeeper final + : public BasicConstraintKeeper { public: /// Constructor, adds this CK to the provided ConstraintManager /// Requires the CM to be already constructed @@ -748,6 +817,11 @@ class ConstraintKeeper : public BasicConstraintKeeper { MarkAsUnused(cons_.at(i), i); } + /// Is constraint \a i unused? + bool IsUnused(int i) const override { + return cons_.at(i).IsUnused(); + } + /// Copy names from ValueNodes void CopyNamesFromValueNodes() override { const auto& vn = GetValueNode().GetStrVec(); @@ -766,6 +840,14 @@ class ConstraintKeeper : public BasicConstraintKeeper { MarkAsBridged(cons_[i], i); } + /// Compute result for constraint \a i + /// (for functional constraints). + double + ComputeValue(int i, const VarInfoRecomp& vir) override { + assert(cons_[i].con_.GetResultVar() >= 0); + return mp::ComputeValue(cons_[i].con_, vir); + } + /// Compute violations for this constraint type. /// We do it for redefined (intermediate) ones too. void ComputeViolations(SolCheck& chk) override { diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h index d20c2479c..aebad1cb6 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -607,12 +607,54 @@ class FlatConverter : } - /// Check unpostsolved solution + /// in various ways bool CheckSolution( ArrayRef x, const pre::ValueMapDbl& duals, ArrayRef obj) { + DoCheckSol(x, duals, obj, "SolutionCheck_Aux"); + auto x1 = RecomputeAuxVars(x); + DoCheckSol(x1, duals, obj, "SolutionCheck"); + return true; + } + + /// Functor to recompute auxiliary var \a i + VarsRecomputeFn recomp_fn + = [this](int i, const VarInfoRecomp& x) { + if (MPCD( HasInitExpression(i) )) { + const auto& iexpr = MPCD( GetInitExpression(i) ); + assert(iexpr.GetCK() + ->GetResultVar(iexpr.GetIndex()) == i); + if (!iexpr.GetCK()->IsUnused(iexpr.GetIndex())) + return iexpr.GetCK() + ->ComputeValue(iexpr.GetIndex(), x); + } + return x.get_x().get_x()[i]; // no recomputation + }; + + /// Recompute auxiliary variables + ArrayRef RecomputeAuxVars(ArrayRef x) { + VarInfoRecomp vir { + options_.solfeastol_, + {x, recomp_fn}, + GetModel().var_type_vec(), + GetModel().var_lb_vec(), + GetModel().var_ub_vec(), + nullptr, nullptr + }; + vir.get_x().set_p_var_info(&vir); + for (auto i=vir.size(); i--; ) + vir[i]; // touch the variable to be recomputed + return std::move(vir.get_x().get_x()); + } + + /// Check single unpostsolved solution + bool DoCheckSol( + ArrayRef x, + const pre::ValueMapDbl& duals, + ArrayRef obj, + const char* warning_header) { SolCheck chk(x, duals, obj, GetModel().var_type_vec(), GetModel().var_lb_vec(), @@ -641,7 +683,7 @@ class FlatConverter : chk.GetReport()); else AddWarning( - "SolutionCheck", chk.GetReport(), + warning_header, chk.GetReport(), true); // replace for multiple solutions } return !chk.HasAnyViols(); diff --git a/include/mp/flat/expr_affine.h b/include/mp/flat/expr_affine.h index f5ddb133a..dbbad883b 100644 --- a/include/mp/flat/expr_affine.h +++ b/include/mp/flat/expr_affine.h @@ -74,7 +74,8 @@ class LinTerms { const LinTerms& GetLinTerms() const { return *this; } /// Compute value given a dense vector of variable values - long double ComputeValue(const ArrayRef& x) const { + template + long double ComputeValue(const VarInfo& x) const { long double s=0.0; for (size_t i=coefs().size(); i--; ) s += (long double)(coefs()[i]) * x[vars()[i]]; diff --git a/include/mp/flat/expr_algebraic.h b/include/mp/flat/expr_algebraic.h index af7cb681a..466f0d9c9 100644 --- a/include/mp/flat/expr_algebraic.h +++ b/include/mp/flat/expr_algebraic.h @@ -76,7 +76,8 @@ class AlgebraicExpression : public Body { void add_to_constant(double a) { constant_term_ += a; } /// Compute value given a dense vector of variable values - double ComputeValue(const ArrayRef& x) const + template + double ComputeValue(const VarInfo& x) const { return Body::ComputeValue(x) + constant_term(); } /// Negate the ae diff --git a/include/mp/flat/expr_quadratic.h b/include/mp/flat/expr_quadratic.h index bba4b4e7d..d2f426b97 100644 --- a/include/mp/flat/expr_quadratic.h +++ b/include/mp/flat/expr_quadratic.h @@ -54,7 +54,8 @@ class QuadTerms { int var2(int i) const { return vars2_[i]; } /// Compute value given a dense vector of variable values - long double ComputeValue(const ArrayRef& x) const { + template + long double ComputeValue(const VarInfo& x) const { long double s=0.0; for (size_t i=coefs().size(); i--; ) s += (long double)(coefs()[i]) * x[vars1()[i]] * x[vars2()[i]]; @@ -200,7 +201,8 @@ class QuadAndLinTerms : } /// Value at given variable vector - long double ComputeValue(const ArrayRef& x) const { + template + long double ComputeValue(const VarInfo& x) const { return LinTerms::ComputeValue(x) + QuadTerms::ComputeValue(x); }