From 81d8fb1a662f4e30251fbcdfc4c08b20417c9824 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Tue, 22 Aug 2023 10:14:28 +1000 Subject: [PATCH] Toward #200: options sol:chk:... --- include/mp/flat/constr_algebraic.h | 5 +- include/mp/flat/constr_base.h | 3 + include/mp/flat/constr_general.h | 9 +++ include/mp/flat/constr_keeper.h | 124 ++++++++++++++++++++++++++++- include/mp/flat/converter.h | 89 ++++++++++++++++++++- include/mp/flat/converter_model.h | 11 ++- include/mp/valcvt.h | 32 +++++++- 7 files changed, 262 insertions(+), 11 deletions(-) diff --git a/include/mp/flat/constr_algebraic.h b/include/mp/flat/constr_algebraic.h index 42ef7d810..4893447a6 100644 --- a/include/mp/flat/constr_algebraic.h +++ b/include/mp/flat/constr_algebraic.h @@ -31,6 +31,9 @@ class AlgebraicConstraint : return name; } + /// Is logical? + static bool IsLogical() { return false; } + /// BodyType using BodyType = Body; @@ -39,7 +42,7 @@ class AlgebraicConstraint : /// Constructor. /// By default (\a fSort = true), it sorts terms. - /// Pass \a fSort = false to skip if you complement the terms list + /// Pass \a fSort = false to skip if you populate the terms /// but do sorting later. /// @param le: linear / linear + higher-order terms /// @param rr: rhs or range diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index 9bd70589f..fc22850d6 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -9,6 +9,7 @@ #include #include "mp/flat/context.h" +#include "mp/arrayref.h" namespace mp { @@ -34,6 +35,8 @@ class BasicConstraint { void AddContext(Context ) const { } /// For functional constraints, result variable index int GetResultVar() const { return -1; } + /// Compute violation + double ComputeViolation(const ArrayRef& ) { return 0.0; } private: std::string name_; diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 4a855b159..9e16a9361 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -27,6 +27,9 @@ class IndicatorConstraint: public BasicConstraint { return name; } + /// Is logical? + static bool IsLogical() { return true; } + /// Constructor IndicatorConstraint(int b, int bv, Con con) noexcept : b_(b), bv_(bv), con_(std::move(con)) { assert(check()); } @@ -99,6 +102,9 @@ class SOS_1or2_Constraint: public BasicConstraint { static const char* GetTypeName() { return 1==type ? name1_ : name2_; } + /// Is logical? + static bool IsLogical() { return true; } + int get_sos_type() const { return type; } int size() const { return (int)v_.size(); } /// Returns vector of variables, sorted by weights @@ -175,6 +181,9 @@ class ComplementarityConstraint : public BasicConstraint { return name; } + /// Is logical? + static bool IsLogical() { return true; } + /// Constructor ComplementarityConstraint(ExprType expr, int var) : compl_expr_(std::move(expr)), compl_var_(var) { } diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h index ea026ba88..99095e773 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -26,6 +26,91 @@ static const mp::OptionValueInfo values_item_acceptance[] = { { "2", "Accepted natively and preferred", 2} }; + +/// 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; + } + } + /// Count violation + void CountViol(double val) { + ++N_; + if (epsMax_ < val) + epsMax_ = val; + } + int N_ {0}; + double epsMax_ {0.0}; +}; + +/// Array of violation summaries. +/// For different kinds, e.g., original / aux vars. +template +using ViolSummArray = std::array; + +/// Solution check data +struct SolCheck { + /// Construct + SolCheck(ArrayRef x, + const pre::ValueMapDbl& duals, + ArrayRef obj, + double feastol, double inttol) + : x_(x), y_(duals), obj_(obj), + feastol_(feastol), inttol_(inttol) { } + /// Any violations? + bool HasAnyViols() const { return hasAnyViol_; } + /// Summary + const std::string& GetReport() const { return report_; } + + /// x + ArrayRef& x() { return x_; } + /// x[i] + double x(int i) const { return x_[i]; } + /// Feasibility tolerance + double GetFeasTol() const { return feastol_; } + + /// Var bnd violations + ViolSummArray<2>& VarViolBnds() { return viol_var_bnds_; } + /// Var int-ty violations + ViolSummArray<2>& VarViolIntty() { return viol_var_int_; } + + /// Constraints: algebraic. + /// Map by constraint type. + /// Values: for original, intermediate, + /// and solver-side constraints. + std::map< std::string, ViolSummArray<3> >& + ConViolAlg() { return viol_cons_alg_; } + /// Constraints: logical. + std::map< std::string, ViolSummArray<3> >& + ConViolLog() { return viol_cons_log_; } + +private: + ArrayRef x_; + const pre::ValueMapDbl& y_; + ArrayRef obj_; + double feastol_; + double inttol_; + + bool hasAnyViol_ = false; + std::string report_; + + /// Variable bounds: orig, aux + ViolSummArray<2> viol_var_bnds_; + /// Variable integrality: orig, aux + ViolSummArray<2> viol_var_int_; + /// Constraints: algebraic. + std::map< std::string, ViolSummArray<3> > viol_cons_alg_; + /// Constraints: logical. + std::map< std::string, ViolSummArray<3> > viol_cons_log_; + /// Objectives + ViolSummary viol_obj_; +}; + + /// Interface for an array of constraints of certain type class BasicConstraintKeeper { public: @@ -174,6 +259,9 @@ class BasicConstraintKeeper { /// Copy names from ValueNodes virtual void CopyNamesFromValueNodes() = 0; + /// Compute violations + virtual void ComputeViolations(SolCheck& ) = 0; + protected: int& GetAccLevRef() { return acceptance_level_; } @@ -419,11 +507,16 @@ class ConstraintKeeper : public BasicConstraintKeeper { struct Container { Container(Constraint&& c) noexcept : con_(std::move(c)) { } + bool IsDeleted() const { return IsBridged(); } bool IsBridged() const { return is_bridged_; } void MarkAsBridged() { is_bridged_=true; } + /// Depth in the redef tree + int GetDepth() const { return depth_; } + Constraint con_; bool is_bridged_ = false; + int depth_ {0}; }; /// Convert all new constraints of this type @@ -492,7 +585,7 @@ class ConstraintKeeper : public BasicConstraintKeeper { } /// ForEachActive(). - /// Deletes every constraint where fn() returned true. + /// Deletes every constraint where fn() returns true. template void ForEachActive(Fn fn) { for (int i=0; i<(int)cons_.size(); ++i) @@ -501,6 +594,29 @@ class ConstraintKeeper : public BasicConstraintKeeper { MarkAsDeleted(cons_[i], i); } + /// Compute violations for this constraint type. + /// We do it for redefined 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()]; + const auto& x = chk.x(); + for (int i=(int)cons_.size(); i--; ) { + auto viol = cons_[i].con_.ComputeViolation(x); + if (viol > chk.GetFeasTol()) { + /// Solver-side? + /// TODO also original NL constraints (index 0) + int index = cons_[i].IsDeleted() ? 2 : 1; + conviolarray[index].CountViol(viol); + } + } + } + } + protected: /// Add all non-converted items to ModelAPI @@ -725,6 +841,12 @@ class ConstraintManager { for (const auto& ck: con_keepers_) ck.second.AddUnbridgedToBackend(be); } + + /// Compute violations + void ComputeViolations(SolCheck& chk) { + for (const auto& ck: con_keepers_) + ck.second.ComputeViolations(chk); + } }; } // namespace mp diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h index fbcc696d6..5d8a0f91b 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -601,6 +601,59 @@ class FlatConverter : } + + /// Check unpostsolved solution + bool CheckSolution( + ArrayRef x, + const pre::ValueMapDbl& duals, + ArrayRef obj) { + SolCheck chk(x, duals, obj, + options_.solfeastol_, options_.solinttol_); + CheckVars(chk); + CheckCons(chk); + CheckObjs(chk); + GenerateViolationsReport(chk); + // What if this is an intermediate solution? + // Should be fine - warning by default, + // fail if requested explicitly. + // If warning, we should add the report + // to that solutions' solve message, and + // a summary in the final solve message. + // For now, do this via warnings? + if (chk.HasAnyViols()) { + if (options_.solcheckfail_) + MP_RAISE(chk.GetReport()); + else + AddWarning("SolutionCheck", chk.GetReport()); + } + return !chk.HasAnyViols(); + } + + void CheckVars(SolCheck& chk) { + for (auto i=num_vars(); i--; ) { + auto x = chk.x(i); + bool aux = !MPCD( is_var_original(i) ); + chk.VarViolBnds().at(aux).CheckViol( + MPCD( lb(i) ) - x, options_.solfeastol_); + chk.VarViolBnds().at(aux).CheckViol( + x - MPCD( ub(i) ), options_.solfeastol_); + if (is_var_integer(i)) + chk.VarViolIntty().at(aux).CheckViol( + std::fabs(x - std::round(x)), + options_.solinttol_); + } + } + + /// Includes logical constraints. + void CheckCons(SolCheck& chk) { + GetModel().ComputeViolations(chk); + } + + void CheckObjs(SolCheck& ) { + } + + void GenerateViolationsReport(SolCheck& chk) { } + //////////////////////////// UTILITIES ///////////////////////////////// /// public: @@ -669,8 +722,11 @@ class FlatConverter : public: - /// Shortcut num_vars() - int num_vars() const { return MPCD(GetModel()).num_vars(); } + /// Shortcut num_vars() + int num_vars() const { return MPCD(GetModel()).num_vars(); } + /// Shortcut is_var_original() + int is_var_original(int i) const + { return MPCD(GetModel()).is_var_original(i); } /// Shortcut lb(var) double lb(int var) const { return this->GetModel().lb(var); } /// Shortcut ub(var) @@ -874,7 +930,7 @@ class FlatConverter : return ModelAPI::AcceptsNonconvexQC(); } - /// Whether the ModelAPI accepts quadragtic cones + /// Whether the ModelAPI accepts quadratic cones int ModelAPIAcceptsQuadraticCones() { return std::max( @@ -902,6 +958,10 @@ class FlatConverter : int passExpCones_ = 0; int relax_ = 0; + + bool solcheckfail_ = false; + double solfeastol_ = 1e-6; + double solinttol_ = 1e-6; }; Options options_; @@ -975,6 +1035,19 @@ class FlatConverter : GetEnv().AddOption("alg:relax relax", "0*/1: Whether to relax integrality of variables.", options_.relax_, 0, 1); + GetEnv().AddOption("sol:chk:fail solcheckfail checkfail chk:fail", + "Fail if solution violates tolerances " + "(normally only warning).", + options_.solcheckfail_, false, true); + GetEnv().AddOption("sol:chk:feastol sol:chk:eps sol:eps chk:eps", + "Solution checking tolerance for objective values, variable " + "and constraint bounds. Default 1e-6. " + "Violated logical constraints are always reported.", + options_.solfeastol_, 0.0, 1e100); + GetEnv().AddOption("sol:chk:inttol sol:chk:inteps sol:inteps chk:inteps", + "Solution checking tolerance for variables' integrality. " + "Default 1e-6. ", + options_.solinttol_, 0.0, 1e100); } @@ -1040,7 +1113,15 @@ class FlatConverter : /// ValuePresolver: should be init before constraint keepers /// and links pre::ValuePresolver value_presolver_ - {GetModel(), GetEnv(), graph_exporter_fn_}; + { + GetModel(), GetEnv(), graph_exporter_fn_, + [this]( + ArrayRef x, + const pre::ValueMapDbl& y, + ArrayRef obj) -> bool { + return this->CheckSolution(x, y, obj); + } + }; pre::CopyLink copy_link_ { GetValuePresolver() }; // the copy links pre::One2ManyLink one2many_link_ { GetValuePresolver() }; // the 1-to-many links pre::NodeRange auto_link_src_item_; // the source item for autolinking diff --git a/include/mp/flat/converter_model.h b/include/mp/flat/converter_model.h index c88118023..39e25972e 100644 --- a/include/mp/flat/converter_model.h +++ b/include/mp/flat/converter_model.h @@ -45,7 +45,8 @@ class FlatModel return var_type_.size()-1; } - /// Add several variables + /// Add several variables. + /// This is only done once when we add original variables. void AddVars__basic(const VarBndVec& lbs, const VarBndVec& ubs, const VarTypeVec& types) { assert(check_vars()); @@ -53,6 +54,7 @@ class FlatModel var_ub_.insert(var_ub_.end(), ubs.begin(), ubs.end()); var_type_.insert(var_type_.end(), types.begin(), types.end()); assert(check_vars()); + num_vars_orig_ = var_lb_.size(); } void AddVarNames(const std::vector& names) { @@ -62,6 +64,9 @@ class FlatModel int num_vars() const { assert(check_vars()); return (int)var_lb_.size(); } + bool is_var_original(int i) const + { return i /// Variables' names - /// mutable VarNameVec var_names_; std::vector var_names_storage_; + /// Number of original NL variables + int num_vars_orig_ {0}; /// Objectives ObjList objs_; diff --git a/include/mp/valcvt.h b/include/mp/valcvt.h index 3922b3dac..9bb65757c 100644 --- a/include/mp/valcvt.h +++ b/include/mp/valcvt.h @@ -251,12 +251,22 @@ class ValuePresolverImpl : public BasicValuePresolver { }; +/// How to call a solution checker +using SolCheckerCall = bool( + ArrayRef x, + const ValueMapDbl& y, + ArrayRef obj); +/// Function wrapper +using SolCheckerType = std::function; + + /// Final ValuePresolver. /// It specializes PresolveSolution() to update fixed variables class ValuePresolver : public ValuePresolverImpl { public: - ValuePresolver(BasicFlatModel& m, Env& env, ExporterFn bts={}) - : ValuePresolverImpl(env, bts), model_(m) { } + ValuePresolver(BasicFlatModel& m, Env& env, + ExporterFn bts={}, SolCheckerType sc={}) + : ValuePresolverImpl(env, bts), model_(m), solchk_(sc) { } /// Override PresolveSolution(). /// Move warm start values into the bounds. @@ -278,9 +288,27 @@ class ValuePresolver : public ValuePresolverImpl { return result; } + /// Override PostsolveSolution(). + /// Check solution if checker provided. + MVOverEl PostsolveSolution ( + const MVOverEl& mv) override { + if (solchk_) { + const auto& mx = mv.GetVarValues(); + const auto& mo = mv.GetObjValues(); + if (mx.IsSingleKey() + && mx().size() // solution available + && mo.IsSingleKey()) { + solchk_(mx(), + mv.GetConValues(), mo()); + } + } + return ValuePresolverImpl::PostsolveSolution(mv); + } + private: BasicFlatModel& model_; + SolCheckerType solchk_; };