Check objective values #200

Disable idealistic check by default, 'experimental'

include/mp/flat/constr_algebraic.h

@@ -71,13 +71,19 @@ class AlgebraicConstraint :
   }
 
   /// Compute violation.
-  /// Negative if holds with slack.
+  /// For checking solver values,
+  /// report violation amount
+  /// (negative if holds with slack.)
+  /// In idealistic mode, report 1/0
+  /// (what's the violation amount of 0 for >0?)
   template <class VarInfo>
   double
   ComputeViolation(const VarInfo& x) const {
     auto bd = Body::ComputeValue(x);
-    return std::max(      // same for strict cmp?
+    if (!x.idealistic())
+      return std::max(      // same for strict cmp?
           RhsOrRange::lb() - bd, bd - RhsOrRange::ub());
+    return !RhsOrRange::is_valid(bd);
   }
 
   /// Sorting and merging terms, some solvers require
@@ -115,6 +121,9 @@ class AlgConRange {
   /// operator==
   bool equals(const AlgConRange& r) const
   { return lb()==r.lb() && ub()==r.ub(); }
+  /// validity of the body value
+  bool is_valid(double bv) const
+  { return bv >= lb() && bv <= ub(); }
 private:
   double lb_, ub_;
 };
@@ -153,6 +162,17 @@ class AlgConRhs {
   /// operator==
   bool equals(const AlgConRhs& r) const
   { return rhs()==r.rhs(); }
+  /// validity of the body value
+  bool is_valid(double bv) const {
+    switch (kind_) {
+    case -2: return bv <  rhs();
+    case -1: return bv <= rhs();
+    case  0: return bv == rhs();
+    case  1: return bv >= rhs();
+    case  2: return bv >  rhs();
+    default: MP_RAISE("wrong comparison kind");
+    }
+  }
 private:
   double rhs_;
 };

include/mp/flat/constr_base.h

@@ -195,7 +195,9 @@ double ComputeViolation(
       return INFINITY;
     }
   }
-  return x.bounds_viol(resvar);  // recomputed vars: bounds viol
+  return                              // recomputed var minus solver's
+      std::fabs(x[resvar] - x.raw(resvar))
+      + x.bounds_viol(resvar);
 }
 
 

include/mp/flat/constr_keeper.h

@@ -113,11 +113,13 @@ class VarInfoImpl {
   /// Constructor
   VarInfoImpl(double ft, bool recomp_vals,
           VarVec x,
+          ArrayRef<double> x_raw,
           ArrayRef<var::Type> type,
           ArrayRef<double> lb, ArrayRef<double> ub,
           const char* sol_rnd, const char* sol_prec)
     : feastol_(ft), recomp_vals_(recomp_vals),
-      x_(std::move(x)), type_(type), lb_(lb), ub_(ub) {
+      x_(std::move(x)), x_raw_(x_raw),
+      type_(type), lb_(lb), ub_(ub) {
     assert((int)x_.size()>=(int)type_.size());  // feasrelax can add more
     assert(type_.size()==lb_.size());
     assert(type_.size()==ub_.size());
@@ -132,6 +134,13 @@ class VarInfoImpl {
   }
   /// Access VarVec
   const VarVec& get_x() const { return x_; }
+  /// Access raw variables
+  double raw(int i) const {
+    assert(i < x_raw_.size()
+           && "Can only access raw solver values "
+              "in idealistic mode and they should be available");
+    return x_raw_[i];
+  }
 
   /// Access integrality condition
   bool is_var_int(int i) const {
@@ -167,6 +176,9 @@ class VarInfoImpl {
   double feastol() const { return feastol_; }
   /// Using recomputed auxiliary vars?
   bool recomp_vals() const { return recomp_vals_; }
+  /// Using idealistic checking of solution
+  /// (without tolerances)?
+  bool idealistic() const { return recomp_vals(); }
   /// sol_rnd as string
   std::string solution_round() const
   { return sol_rnd_ < 100 ? std::to_string(sol_rnd_) : ""; }
@@ -200,6 +212,7 @@ class VarInfoImpl {
   bool recomp_vals_;   // variables are recomputed
 
   VarVec x_;   // can be rounded, recomputed, etc.
+  ArrayRef<double> x_raw_;     // solver values
   const ArrayRef<var::Type> type_;
   const ArrayRef<double> lb_;
   const ArrayRef<double> ub_;
@@ -229,13 +242,14 @@ struct SolCheck {
   SolCheck(ArrayRef<double> x,
            const pre::ValueMapDbl& duals,
            ArrayRef<double> obj,
+           ArrayRef<double> x_raw,
            ArrayRef<var::Type> vtype,
            ArrayRef<double> lb,  ArrayRef<double> ub,
            double feastol, double inttol,
            const char* sol_rnd, const char* sol_prec,
            bool recomp_vals, int chk_mode)
     : x_(feastol, recomp_vals,
-         x, vtype, lb, ub, sol_rnd, sol_prec),
+         x, x_raw, vtype, lb, ub, sol_rnd, sol_prec),
       y_(duals), obj_(obj),
       feastol_(feastol), inttol_(inttol),
       fRecomputedVals_(recomp_vals),
@@ -261,6 +275,10 @@ struct SolCheck {
   const VarInfoStatic& x_ext() const { return x_; }
   /// x[i]
   double x(int i) const { return x_[i]; }
+  /// objective values
+  const ArrayRef<double>& obj_vals() const
+  { return obj_; }
+
   /// Feasibility tolerance
   double GetFeasTol() const { return feastol_; }
 

include/mp/flat/converter.h

@@ -618,13 +618,14 @@ class FlatConverter :
     bool result = true;
     std::string err_msg;
     try {            // protect
+      std::vector<double> x_back; // to extract x used
       if (options_.solcheckmode_ & (1+2+4+8+16)) {
-        if (!DoCheckSol(x, duals, obj, false))
+        if (!DoCheckSol(x, duals, obj, {}, x_back, false))
           result = false;
       }
       if (options_.solcheckmode_ & (32+64+128+256+512)) {
         auto x1 = RecomputeAuxVars(x);
-        if (!DoCheckSol(x1, duals, obj, true))
+        if (!DoCheckSol(x1, duals, obj, x_back, x_back, true))
           result = false;
       }
     } catch (const mp::Error& err) {
@@ -664,6 +665,7 @@ class FlatConverter :
     VarInfoRecomp vir {
       options_.solfeastol_, false,
       {x, recomp_fn},
+      {},              // now raw values
       GetModel().var_type_vec(),
           GetModel().var_lb_vec(),
           GetModel().var_ub_vec(),
@@ -680,8 +682,10 @@ class FlatConverter :
       ArrayRef<double> x,
       const pre::ValueMapDbl& duals,
       ArrayRef<double> obj,
+      ArrayRef<double> x_raw,
+      std::vector<double>& x_back,
       bool if_recomp_vals) {
-    SolCheck chk(x, duals, obj,
+    SolCheck chk(x, duals, obj, x_raw,
                  GetModel().var_type_vec(),
                  GetModel().var_lb_vec(),
                  GetModel().var_ub_vec(),
@@ -719,6 +723,7 @@ class FlatConverter :
               chk.GetReport(),
               true);  // replace for multiple solutions
     }
+    x_back = chk.x_ext().get_x();
     return !chk.HasAnyViols();
   }
 
@@ -750,7 +755,20 @@ class FlatConverter :
     GetModel().ComputeViolations(chk);
   }
 
-  void CheckObjs(SolCheck& ) { }
+  void CheckObjs(SolCheck& chk) {
+    const auto& objs = GetModel().get_objectives();
+    // Solvers might have dummy obj.
+    // Unbounded problems might have no obj value.
+    for (auto i
+         =std::min(objs.size(), chk.obj_vals().size());
+         i--; ) {
+      chk.ObjViols().CheckViol(
+            std::fabs(chk.obj_vals()[i]
+                      - ComputeValue(objs[i], chk.x_ext())),
+            options_.solfeastol_,
+            objs[i].name());
+    }
+  }
 
   void GenerateViolationsReport(SolCheck& chk) {
     fmt::MemoryWriter wrt;
@@ -1128,7 +1146,7 @@ class FlatConverter :
 
     int relax_ = 0;
 
-    int solcheckmode_ = 1+2+16+512;
+    int solcheckmode_ = 1+2+16;
     bool solcheckfail_ = false;
     double solfeastol_ = 1e-6;
     double solinttol_ = 1e-5;
@@ -1223,12 +1241,12 @@ class FlatConverter :
         "| 32, 64, 128, 256, 512 - similar, but "
         "      non-linear expressions are recomputed "
         "      (vs using their values reported by the solver.) "
-        "      This is an idealistic check, because "
+        "      *Experimental.* This is an idealistic check, because "
         "      it does not consider possible tolerances "
         "      applied by the solver when computing "
         "      expression values.\n"
                              "\n"
-                             "Default: 1+2+16+512.",
+                             "Default: 1+2+16.",
         options_.solcheckmode_, 0, 1024);
     GetEnv().AddOption("sol:chk:feastol sol:chk:eps sol:eps chk:eps",
         "Solution checking tolerance for objective values, variable "

include/mp/flat/converter_model.h

@@ -174,7 +174,7 @@ class FlatModel
   }
 
   ///////////////////////////// OBJECTIVES ////////////////////////////
-protected:
+public:
   using ObjList = std::vector<QuadraticObjective>;
   const ObjList& get_objectives() const { return objs_; }
   ObjList& get_objectives() { return objs_; }

include/mp/flat/obj_std.h

@@ -6,14 +6,14 @@
 
 #include "mp/common.h"
 #include "mp/flat/expr_quadratic.h"
+#include "mp/flat/constr_eval.h"
 
 namespace mp {
 
 /// Linear objective incl. sense and name
 class LinearObjective {
   obj::Type sense_;
-  std::vector<double> coefs_;
-  std::vector<int> vars_;
+  LinTerms lt_;
   std::string name_;
 public:
   /// Construct
@@ -22,26 +22,28 @@ class LinearObjective {
   LinearObjective(obj::Type s, CoefVec&& c, VarVec&& v,
                   std::string nm = {}) noexcept :
     sense_(s),
-    coefs_(std::forward<CoefVec>(c)), vars_(std::forward<VarVec>(v)),
+    lt_(std::forward<CoefVec>(c), std::forward<VarVec>(v)),
     name_(std::move(nm)){ }
   /// Get sense
   obj::Type obj_sense() const { return sense_; }
+  /// Get lin terms
+  const LinTerms& GetLinTerms() const { return lt_; }
   /// Get N terms
-  int num_terms() const { assert(check()); return (int)vars_.size(); }
+  int num_terms() const { assert(check()); return lt_.size(); }
   /// Validate
-  bool check() const { return coefs_.size()==vars_.size(); }
+  bool check() const { return lt_.check(); }
   /// Coefs vector
-  const std::vector<double>& coefs() const { return coefs_; }
+  const std::vector<double>& coefs() const { return lt_.coefs(); }
   /// Var vector
-  const std::vector<int>& vars() const { return vars_; }
+  const std::vector<int>& vars() const { return lt_.vars(); }
   /// Name
   const char* name() const { return name_.c_str(); }
   /// Set name
   void set_name(std::string nm) { name_ = std::move(nm); }
 
   /// Testing API
   bool operator==(const LinearObjective& lc) const {
-    return sense_==lc.sense_ && coefs_==lc.coefs_ && vars_==lc.vars_;
+    return sense_==lc.sense_ && lt_==lc.lt_;
   }
 };
 
@@ -67,6 +69,15 @@ class QuadraticObjective : public LinearObjective {
   }
 };
 
+/// Compute value of an objective.
+template <class VarVec>
+double ComputeValue(
+    const QuadraticObjective& obj, const VarVec& x) {
+  return
+      obj.GetLinTerms().ComputeValue(x)
+      + obj.GetQPTerms().ComputeValue(x);
+}
+
 } // namespace mp
 
 #endif // STD_OBJ_H

test/end2end/cases/categorized/fast/tech/modellist.json

@@ -32,7 +32,7 @@
     "files" : ["infeas_int_01.mod"],
     "options": {
       "ANYSOLVER_options":
-        "sol:chk:fail sol:chk:feastol=1e-7 sol:chk:inttol=1e-4 feasrelax=1"
+        "sol:chk:fail sol:chk:mode=1023 sol:chk:feastol=1e-7 sol:chk:inttol=1e-4 feasrelax=1"
     },
     "values": {
       "solve_result_num": 520