Toward #200: options sol:chk:...

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

include/mp/flat/constr_base.h

@@ -9,6 +9,7 @@
 #include <typeinfo>
 
 #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<double>& ) { return 0.0; }
 
 private:
   std::string name_;

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) { }

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 <int Nkinds>
+using ViolSummArray = std::array<ViolSummary, Nkinds>;
+
+/// Solution check data
+struct SolCheck {
+  /// Construct
+  SolCheck(ArrayRef<double> x,
+           const pre::ValueMapDbl& duals,
+           ArrayRef<double> 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<double>& 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<double> x_;
+  const pre::ValueMapDbl& y_;
+  ArrayRef<double> 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 <class Fn>
 	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

include/mp/flat/converter.h

@@ -601,6 +601,59 @@ class FlatConverter :
   }
 
 
+
+  /// Check unpostsolved solution
+  bool CheckSolution(
+      ArrayRef<double> x,
+      const pre::ValueMapDbl& duals,
+      ArrayRef<double> 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<double> x,
+        const pre::ValueMapDbl& y,
+        ArrayRef<double> 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