SolCheck: fail to read 'solution_round', 'solution_precision' #200

AMPL does not put them into environment actually

CMakeLists.txt

@@ -327,7 +327,9 @@ add_prefix(MP_HEADERS include/mp/
   valcvt-base.h valcvt-node.h valcvt-link.h valcvt.h
   rstparser.h safeint.h sol.h
   solver.h solver-opt.h solver-base.h solver-io.h solver-app-base.h solver-app.h
-  suffix.h utils-file.h utils-hash.h utils-hash-stream.h utils-string.h)
+  suffix.h
+  utils-file.h utils-hash.h utils-hash-stream.h
+  utils-string.h utils-math.h)
 
 add_prefix(MP_FLAT_HEADERS include/mp/flat/
   backend_flat.h

include/mp/flat/constr_keeper.h

@@ -9,6 +9,7 @@
 #include "mp/common.h"
 #include "mp/format.h"
 #include "mp/env.h"
+#include "mp/utils-math.h"
 
 #include "mp/flat/model_api_base.h"
 #include "mp/flat/constr_hash.h"
@@ -62,12 +63,14 @@ class VarInfo {
   VarInfo(double ft,
           std::vector<double> x,
           ArrayRef<var::Type> type,
-          ArrayRef<double> lb, ArrayRef<double> ub)
+          ArrayRef<double> lb, ArrayRef<double> ub,
+          const char* sol_rnd, const char* sol_prec)
     : feastol_(ft), x_(std::move(x)), x_ref_(x_),
       type_(type), lb_(lb), ub_(ub) {
     assert(x_.size()>=type_.size());  // feasrelax can add more
     assert(type_.size()==lb_.size());
     assert(type_.size()==ub_.size());
+    apply_precision_options(sol_rnd, sol_prec);
   }
   /// Access variable value
   double operator[]( int i ) const {
@@ -103,17 +106,47 @@ class VarInfo {
 
   /// Feasibility tolerance
   double feastol() const { return feastol_; }
+  /// sol_rnd as string
+  std::string solution_round() const
+  { return sol_rnd_ < 100 ? std::to_string(sol_rnd_) : ""; }
+  /// sol_rnd as string
+  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) {
+    try {                 // Apply sol_rnd
+      if (sol_rnd) {
+        sol_rnd_ = std::stoi(sol_rnd);
+        auto scale = std::pow(10, sol_rnd_);
+        auto scale_rec = 1.0/scale;
+        for (auto& x: x_)
+          x = std::round(x * scale) * scale_rec;
+      }
+    } catch (...) { sol_rnd_=100; }     // Could add a warning
+    try {                 // Apply sol_prec
+      if (sol_prec) {
+        sol_prec_ = std::stoi(sol_prec);
+        for (auto& x: x_)
+          x = round_to_digits(x, sol_prec_);
+      }
+    } catch (...) { sol_prec_=100; }     // Could add a warning
+  }
+
 private:
   double feastol_;
 
-  const std::vector<double> x_;   // can be rounded, etc.
+  std::vector<double> x_;   // can be rounded, etc.
   const ArrayRef<double> x_ref_;
   const ArrayRef<var::Type> type_;
   const ArrayRef<double> lb_;
   const ArrayRef<double> ub_;
+  int sol_rnd_=100;    // AMPL option solution_round, if used
+  int sol_prec_=100;   // AMPL option solution_precision, if used
 };
 
 
@@ -126,8 +159,10 @@ struct SolCheck {
            ArrayRef<var::Type> vtype,
            ArrayRef<double> lb,  ArrayRef<double> ub,
            double feastol, double inttol,
+           const char* sol_rnd, const char* sol_prec,
            bool reportBridged)
-    : x_(feastol, x, vtype, lb, ub), y_(duals), obj_(obj),
+    : x_(feastol, x, vtype, lb, ub, sol_rnd, sol_prec),
+      y_(duals), obj_(obj),
       feastol_(feastol), inttol_(inttol),
       reportBridgedCons_(reportBridged) { }
   /// Any violations?

include/mp/flat/converter.h

@@ -619,6 +619,10 @@ class FlatConverter :
                  GetModel().var_ub_vec(),
                  options_.solfeastol_,
                  options_.solinttol_,
+                 options_.dont_use_sol_round_
+                   ? "" : std::getenv("solution_round"),
+                 options_.dont_use_sol_prec_
+                   ? "" : std::getenv("solution_precision"),
                  options_.reprefcons_);
     CheckVars(chk);
     CheckCons(chk);
@@ -674,9 +678,12 @@ class FlatConverter :
     fmt::MemoryWriter wrt;
     if (chk.HasAnyConViols()) {
       wrt.write(
-            "Constraint violations "
-            "(sol:chk:feastol={}, sol:chk:inttol={}):\n",
-            options_.solfeastol_, options_.solinttol_);
+            "Constraint violations\n"
+            "    (sol:chk:feastol={}, sol:chk:inttol={},\n"
+            "     solution_round='{}', solution_precision='{}'):\n",
+            options_.solfeastol_, options_.solinttol_,
+            chk.x_ext().solution_round(),
+            chk.x_ext().solution_precision());
       Gen1Viol(chk.VarViolBnds().at(0), wrt,
                "  - {} original variable(s) violate bounds,\n"
                "      max by {}");
@@ -694,8 +701,11 @@ class FlatConverter :
     GenConViol(chk.ConViolLog(), wrt, "Logical");
     if (chk.HasAnyObjViols()) {
       wrt.write("Objective value violations"
-                "(sol:chk:feastol={})\n",
-                options_.solfeastol_);
+                "    (sol:chk:feastol={},\n"
+                "     solution_round='{}', solution_precision='{}'):\n",
+                options_.solfeastol_,
+                chk.x_ext().solution_round(),
+                chk.x_ext().solution_precision());
       Gen1Viol(chk.ObjViols(), wrt,
                "  - {} objective value(s) violated,\n      max by {}");
     }
@@ -1044,6 +1054,8 @@ class FlatConverter :
     double solfeastol_ = 1e-6;
     double solinttol_ = 1e-5;
     bool reprefcons_ = false;
+    bool dont_use_sol_round_ = false;
+    bool dont_use_sol_prec_ = false;
   };
   Options options_;
 
@@ -1139,6 +1151,12 @@ class FlatConverter :
     GetEnv().AddOption("sol:chk:refcons chk:refcons",
         "Report violations of reformulated constraints.",
         options_.reprefcons_, false, true);
+    GetEnv().AddOption("sol:chk:noround chk:noround chk:no_solution_round",
+        "Don't use AMPL solution_round option when checking.",
+        options_.dont_use_sol_round_, false, true);
+    GetEnv().AddOption("sol:chk:noprec chk:noprec chk:no_solution_precision",
+        "Don't use AMPL solution_precision option when checking.",
+        options_.dont_use_sol_prec_, false, true);
   }
 
 

include/mp/utils-math.h

@@ -0,0 +1,22 @@
+#ifndef MP_UTILS_MATH_H
+#define MP_UTILS_MATH_H
+
+#include <cmath>
+
+namespace mp {
+
+/// https://stackoverflow.com/questions/13094224/a-c-routine-to-round-a-float-to-n-significant-digits
+template <class F>
+F round_to_digits(F value, int digits) {
+    if (value == 0.0) // otherwise it will return 'nan'
+        return 0.0;   // due to the log10() of zero
+
+    F factor = std::pow(10.0,
+                        digits
+                        - std::ceil(std::log10(std::fabs(value))));
+    return std::round(value * factor) / factor;
+}
+
+}  // namespace mp
+
+#endif // MP_UTILS_MATH_H