SolCheck with recomputed aux vars #200

This comes close toAMPL's slacks computation

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 <class VarInfo>
   double
-  ComputeViolation(const ArrayRef<double>& 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());

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<double>& ) const
+  template <class VarInfo>
+  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 <class Args, class Params,
-          class NumOrLogic, class Id, class VarVec>
-double ComputeValue(
-    const CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>& ,
-    const VarVec& ) {
+template <class Con, class VarVec>
+double ComputeValue(const Con& , const VarVec& ) {
   MP_RAISE(fmt::format("ComputeValue({}) not implemented.",
-                       typeid(
-                         CustomFunctionalConstraint<Args,
-                         Params, NumOrLogic, Id>).name()));
+                       typeid(Con).name()));
   return 0.0;
 }
 

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 <class VarVec>
+double ComputeValue(
+    const LinearFunctionalConstraint& con, const VarVec& x) {
+  return con.GetAffineExpr().ComputeValue(x);
+}
+
+/// Compute result of the QuadrFuncCon constraint.
+template <class VarVec>
+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 <class Con, class VarVec>
 double ComputeValue(
     const ConditionalConstraint<Con>& 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.

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<double>& x) const {
+  template <class VarInfo>
+  double ComputeViolation(const VarInfo& x) const {
     assert(b_<(int)x.size());
     auto xb = x[b_];
     if (std::round(xb) == bv_)      // Implication needed

include/mp/flat/constr_keeper.h

@@ -55,31 +55,84 @@ struct ViolSummary {
 template <int Nkinds>
 using ViolSummArray = std::array<ViolSummary, Nkinds>;
 
+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<double> 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<double>::iterator begin() { return x_.begin(); }
+  /// Expose end()
+  std::vector<double>::iterator end() { return x_.end(); }
+  /// Move out x
+  std::vector<double>& get_x() const { return x_; }
+
+private:
+  mutable std::vector<double> x_;
+  mutable std::vector<bool> is_recomp_;
+  VarsRecomputeFn recomp_fn_;
+  mutable const VarInfoRecomp* p_var_info_recomp_{nullptr};
+};
+
+
+/// Static var vector
+using VarVecStatic = std::vector<double>;
+
 
 /// Variable information used by solution check
-class VarInfo {
+template <class VarVec>
+class VarInfoImpl {
 public:
   /// Constructor
-  VarInfo(double ft,
-          std::vector<double> x,
+  VarInfoImpl(double ft,
+          VarVec x,
           ArrayRef<var::Type> type,
           ArrayRef<double> lb, ArrayRef<double> 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<double>& () 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<double>& 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<double> x_;   // can be rounded, etc.
-  const ArrayRef<double> x_ref_;
+  VarVec x_;   // can be rounded, recomputed, etc.
   const ArrayRef<var::Type> type_;
   const ArrayRef<double> lb_;
   const ArrayRef<double> ub_;
@@ -150,6 +199,20 @@ class VarInfo {
 };
 
 
+/// VarInfoRecompTypedef
+using VarInfoRecompTypedef = VarInfoImpl<VarVecRecomp>;
+
+/// Define VarInfoRecomp
+class VarInfoRecomp : public VarInfoRecompTypedef {
+public:
+  /// Inherit ctor's
+  using VarInfoRecompTypedef::VarInfoRecompTypedef;
+};
+
+/// VarInfoStatic
+using VarInfoStatic = VarInfoImpl<VarVecStatic>;
+
+
 /// 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<double>& 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<double> 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 Converter, class Backend, class Constraint>
-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 {

include/mp/flat/converter.h

@@ -607,12 +607,54 @@ class FlatConverter :
   }
 
 
-
   /// Check unpostsolved solution
+  /// in various ways
   bool CheckSolution(
       ArrayRef<double> x,
       const pre::ValueMapDbl& duals,
       ArrayRef<double> 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<double> RecomputeAuxVars(ArrayRef<double> 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<double> x,
+      const pre::ValueMapDbl& duals,
+      ArrayRef<double> 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();

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<double>& x) const {
+  template <class VarInfo>
+  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]];

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<double>& x) const
+  template <class VarInfo>
+  double ComputeValue(const VarInfo& x) const
   { return Body::ComputeValue(x) + constant_term(); }
 
   /// Negate the ae