MO Emulator/Sol Checker: recompute expressions #200 #237 #239

include/mp/flat/constr_2_expr.h

@@ -745,6 +745,18 @@ class Constraints2Expr {
     return false;
   }
 
+  /// Recompute implicit aux vars
+  /// (those corresponding to expressions.)
+  /// needed for MO emulator and sol checker.
+  void RecomputeNLAuxVars(pre::ModelValuesDbl& sol) {
+    if (MPCD( IfWantNLOutput() )) {
+      auto& xx = sol.GetVarValues()();
+      if (xx.size()) {                    // solution available
+        const auto& var_is_proper = MPCD( GetVarProperFlags() );
+        xx = MPD( RecomputeAuxVars(xx, var_is_proper) );
+      }
+    }
+  }
 
 private:
   /// (Argument) variable visitor: mark var as proper

include/mp/flat/converter.h

@@ -1504,8 +1504,11 @@ class FlatConverter :
         return !this->options_.solcheckmode_   // not desired
                || MPD( CheckSolution(x, y, obj, p_extra) );
       },
-      [this](pre::ModelValuesDbl& sol)
-      { MPD( ProcessMOIterationUnpostsolvedSolution(sol) ); }
+      [this](pre::ModelValuesDbl& sol)  // Solution pre-postsolver
+      {
+        MPD( RecomputeNLAuxVars(sol) );
+        MPD( ProcessMOIterationUnpostsolvedSolution(sol) );
+      }
   };
   pre::CopyLink copy_link_ { GetValuePresolver() }; // the copy links
   pre::One2ManyLink one2many_link_ { GetValuePresolver() }; // the 1-to-many links

include/mp/flat/obj_std.h

@@ -114,10 +114,10 @@ void WriteModelItem(fmt::MemoryWriter& wrt, const QuadraticObjective& obj,
 template <class VarVec>
 double ComputeValue(
     const QuadraticObjective& obj, const VarVec& x) {
-  return
-      obj.GetLinTerms().ComputeValue(x)
-         + (obj.HasExpr() ? x[obj.ExprIndex()] : 0.0)
-      + obj.GetQPTerms().ComputeValue(x);
+  double linp = obj.GetLinTerms().ComputeValue(x);
+  double qp = obj.GetQPTerms().ComputeValue(x);
+  double nlp = (obj.HasExpr() ? x[obj.ExprIndex()] : 0.0);
+  return linp + qp + nlp;
 }
 
 } // namespace mp

include/mp/flat/sol_check.h

@@ -16,12 +16,12 @@ namespace mp {
 template <class Impl>
 class SolutionChecker {
 public:
-  /// Check postsolved & unpostsolved solutions
+  /// Check postsolved & unpostsolved (with aux vars) solutions
   /// in various ways.
   /// @param p_extra: !=0 means known infeas solution
   bool CheckSolution(
         ArrayRef<double> x,
-        const pre::ValueMapDbl& duals,
+        const pre::ValueMapDbl& duals,   // unused
         ArrayRef<double> obj,
         void* p_extra) {
     bool fKnownInfeas = (bool)p_extra;
@@ -32,8 +32,6 @@ class SolutionChecker {
     try {                            // protect
       std::vector<double> x_back = x;
       if (MPCD( sol_check_mode() ) & (1+2+4+8+16)) {
-        if (MPCD( IfWantNLOutput() ))
-          x = RecomputeAuxVars(x, true);   // Compute expressions
         msgreal = DoCheckSol(x, duals, obj, {}, x_back, false);
       }
       if (MPCD( sol_check_mode() ) & (32+64+128+256+512)) {
@@ -94,6 +92,7 @@ class SolutionChecker {
     return msgreal.empty() && msgidea.empty();
   }
 
+private:
   /// Functor to recompute auxiliary var \a i
   VarsRecomputeFn recomp_fn
   = [this](int i, const VarInfoRecomp& x) {
@@ -110,33 +109,31 @@ class SolutionChecker {
     return x.get_x().get_x()[i];  // no recomputation
   };
 
+public:
   /// Recompute auxiliary variables
-  /// @param fExprOnly: only NL expressions
+  /// @param is_final: if provided, which vars are final
   ArrayRef<double> RecomputeAuxVars(
-      ArrayRef<double> x, bool fExprOnly=false) {
+      ArrayRef<double> x, std::vector<bool> is_final={}) {
     VarInfoRecomp vir {
       MPCD( sol_feas_tol() ),
           true,        // currently not relevant for recomputation
-      {x, recomp_fn},
-      {},              // now raw values
+      {x, recomp_fn, is_final},
+      {},              // no raw values
       MPCD( GetModel() ).var_type_vec(),
           MPCD( GetModel() ).var_lb_vec(),
           MPCD( GetModel() ).var_ub_vec(),
           MPCD( sol_round() ), MPCD( sol_prec() )
     };
     vir.get_x().set_p_var_info(&vir);
-    if (fExprOnly) {
-      for (auto i=vir.size(); i--; )
-        if ( !MPCD( IsProperVar(i) ) )
-          vir[i];         // touch the variable to be recomputed
-    } else {
-      for (auto i=vir.size(); i--; )
-        vir[i];
-    }
+    for (auto i=vir.size(); i--; )
+      vir[i];         // touch the variable to be recomputed
     return std::move(vir.get_x().get_x());
   }
 
+protected:
   /// Check single unpostsolved solution.
+  /// @param x_raw: for idealistic check,
+  ///   should contain auxiliary values from the solver.
   /// @param x_back: solution vector from realistic mode.
   /// It can be changed by the :round and :prec options.
   /// Its auxiliary vars are compared
@@ -167,7 +164,7 @@ class SolutionChecker {
     if (chk.check_mode() & 16)
       CheckObjs(chk);
     MPD( GenerateViolationsReport(chk, if_recomp_vals) );
-    x_back = chk.x_ext().get_x();   // to reuse 'realistic' vector
+    x_back = std::move(chk.x_ext().get_x()); // to reuse 'realistic' vector
     return chk.GetReport();
   }
 

include/mp/flat/sol_check_data.h

@@ -79,10 +79,15 @@ using VarsRecomputeFn
 class VarVecRecomp {
 public:
   /// Construct
+  /// @param is_final: if provided, which values are final
   VarVecRecomp(std::vector<double> x,
-               VarsRecomputeFn rec_fn)
-      : x_(std::move(x)), is_recomp_(x_.size()),
-      recomp_fn_(rec_fn) { assert(recomp_fn_); }
+               VarsRecomputeFn rec_fn,
+               std::vector<bool> is_final = {})
+      : x_(std::move(x)), is_recomp_(std::move(is_final)),
+      recomp_fn_(rec_fn) {
+    assert(recomp_fn_);
+    is_recomp_.resize(x_.size());   // x_ now
+  }
   /// Set p_var_info_recomp_
   void set_p_var_info(const VarInfoRecomp* p) const
   { p_var_info_recomp_ = p; }
@@ -103,7 +108,7 @@ class VarVecRecomp {
   std::vector<double>::iterator begin() { return x_.begin(); }
   /// Expose end()
   std::vector<double>::iterator end() { return x_.end(); }
-  /// Move out x
+  /// x can be moved out
   std::vector<double>& get_x() const { return x_; }
 
 private:
@@ -119,6 +124,7 @@ using VarVecStatic = std::vector<double>;
 
 
 /// Variable information used by solution check
+/// @param VarVec should provide recomputation info
 template <class VarVec>
 class VarInfoImpl {
 public:
@@ -194,12 +200,13 @@ class VarInfoImpl {
   /// sol_rnd as string
   std::string solution_round() const
   { return sol_rnd_ < 100 ? std::to_string(sol_rnd_) : ""; }
-  /// sol_rnd as string
+  /// sol_prec as string
   std::string solution_precision() const
   { return sol_prec_ < 100 ? std::to_string(sol_prec_) : ""; }
 
 
 protected:
+  /// AMPL options solution_round, solution_prec
   void apply_precision_options(
       int sol_rnd, int sol_prec) {
     try {                 // Apply sol_rnd

include/mp/valcvt.h

@@ -279,9 +279,10 @@ class ValuePresolver : public ValuePresolverImpl {
     return result;
   }
 
-  /// Override PostsolveSolution().
-  /// Check solution if checker provided.
-  /// mv's ExtraData() is passed to the checker.
+  /// Override PostsolveSolution(). As follows:
+  /// 1. Call preposts_() if provided.
+  /// 2. Check solution if checker provided.
+  /// 3. mv's ExtraData() is passed to the checker.
   MVOverEl<double> PostsolveSolution (
       const MVOverEl<double>& mv) override {
     auto mv_copy = mv;