Check algebraic and indicator constraints #200

include/mp/backend-app.h

@@ -88,7 +88,8 @@ int BackendApp::Run(char **argv) {
     // if the solution handler is available.
     GetBackend().ReportError(
           er.exit_code()>=0 ? er.exit_code() : 500,
-          er.what());
+          std::string(GetBackend().long_name()) + ":  "
+          + er.what());
   } catch (const std::exception& ex) {
     // For std::exception, which can be thrown by anything,
     // we try to print the result into .sol file,

include/mp/flat/constr_algebraic.h

@@ -70,6 +70,14 @@ class AlgebraicConstraint :
     return Body::ComputeValue(x) - RhsOrRange::lb();
   }
 
+  /// Compute violation
+  double
+  ComputeViolation(const ArrayRef<double>& x) const {
+    auto bd = Body::ComputeValue(x);
+    return std::max(
+          RhsOrRange::lb() - bd, bd - RhsOrRange::ub());
+  }
+
   /// Sorting and merging terms, some solvers require
   void sort_terms() { Body::sort_terms(); }
 

include/mp/flat/constr_base.h

@@ -36,7 +36,8 @@ class BasicConstraint {
   /// For functional constraints, result variable index
   int GetResultVar() const { return -1; }
   /// Compute violation
-  double ComputeViolation(const ArrayRef<double>& ) { return 0.0; }
+  double ComputeViolation(const ArrayRef<double>& ) const
+  { return 0.0; }
 
 private:
   std::string name_;

include/mp/flat/constr_general.h

@@ -38,13 +38,23 @@ class IndicatorConstraint: public BasicConstraint {
   /// Getters
   int get_binary_var() const { return b_; }
   int get_binary_value() const { return bv_; }
-  bool is_binary_value_1() const { return 1==get_binary_value(); }
+  bool is_binary_value_1() const
+  { return 1==get_binary_value(); }
   const Con& get_constraint() const { return con_; }
 
+  /// Compute violation
+  double ComputeViolation(const ArrayRef<double>& x) const {
+    assert(b_<(int)x.size());
+    auto xb = x[b_];
+    if (std::round(xb) == bv_)      // Implication needed
+      return con_.ComputeViolation(x);
+    return 0.0;
+  }
+
 
 private:
-  const int b_=-1;                            // the indicator variable
-  const int bv_=1;                            // the value, 0/1
+  const int b_=-1;                  // the indicator variable
+  const int bv_=1;                  // the value, 0/1
   const Con con_;
 };
 

include/mp/flat/constr_keeper.h

@@ -30,21 +30,21 @@ static const mp::OptionValueInfo values_item_acceptance[] = {
 /// 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;
-    }
+  void CheckViol(
+      double val, double eps, const char* nm) {
+    if (val > eps)
+      CountViol(val, nm);
   }
   /// Count violation
-  void CountViol(double val) {
+  void CountViol(double val, const char* nm) {
     ++N_;
     if (epsMax_ < val)
       epsMax_ = val;
+    name_ = nm;
   }
   int N_ {0};
   double epsMax_ {0.0};
+  const char* name_ {nullptr};
 };
 
 /// Array of violation summaries.
@@ -62,7 +62,19 @@ struct SolCheck {
     : x_(x), y_(duals), obj_(obj),
       feastol_(feastol), inttol_(inttol) { }
   /// Any violations?
-  bool HasAnyViols() const { return hasAnyViol_; }
+  bool HasAnyViols() const
+  { return HasAnyConViols() || HasAnyObjViols(); }
+  /// Any constraint violations?
+  bool HasAnyConViols() const {
+    return viol_var_bnds_[0].N_ || viol_var_bnds_[1].N_
+        || viol_var_int_[0].N_ || viol_var_int_[1].N_
+        || viol_cons_alg_.size()
+        || viol_cons_log_.size();
+  }
+  /// Any objective value violations?
+  bool HasAnyObjViols() const
+  { return viol_obj_.N_; }
+
   /// Summary
   const std::string& GetReport() const { return report_; }
 
@@ -88,14 +100,20 @@ struct SolCheck {
   std::map< std::string, ViolSummArray<3> >&
   ConViolLog() { return viol_cons_log_; }
 
+  /// Obj viols
+  ViolSummary& ObjViols() { return viol_obj_; }
+
+  /// Set report
+  void SetReport(std::string rep)
+  { report_ = std::move(rep); }
+
 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
@@ -595,23 +613,27 @@ class ConstraintKeeper : public BasicConstraintKeeper {
 	}
 
   /// Compute violations for this constraint type.
-  /// We do it for redefined ones too.
+  /// We do it for redefined (intermediate) 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()];
+            chk.ConViolLog() :
+            chk.ConViolAlg();
       const auto& x = chk.x();
+      ViolSummArray<3>* conviolarray {nullptr};
       for (int i=(int)cons_.size(); i--; ) {
         auto viol = cons_[i].con_.ComputeViolation(x);
         if (viol > chk.GetFeasTol()) {
+          if (!conviolarray)
+            conviolarray =
+                &conviolmap[cons_.front().con_.GetTypeName()];
           /// Solver-side?
           /// TODO also original NL constraints (index 0)
-          int index = cons_[i].IsDeleted() ? 2 : 1;
-          conviolarray[index].CountViol(viol);
+          int index = cons_[i].IsDeleted() ? 1 : 2;
+          assert(index < (int)conviolarray->size());
+          (*conviolarray)[index].CountViol(
+                viol, cons_[i].con_.name());
         }
       }
     }

include/mp/flat/converter.h

@@ -622,7 +622,8 @@ class FlatConverter :
     // For now, do this via warnings?
     if (chk.HasAnyViols()) {
       if (options_.solcheckfail_)
-        MP_RAISE(chk.GetReport());
+        MP_RAISE_WITH_CODE(520,    // numeric error
+                           chk.GetReport());
       else
         AddWarning("SolutionCheck", chk.GetReport());
     }
@@ -634,13 +635,18 @@ class FlatConverter :
       auto x = chk.x(i);
       bool aux = !MPCD( is_var_original(i) );
       chk.VarViolBnds().at(aux).CheckViol(
-                MPCD( lb(i) ) - x, options_.solfeastol_);
+            MPCD( lb(i) ) - x,
+            options_.solfeastol_,
+            GetModel().var_name(i));
       chk.VarViolBnds().at(aux).CheckViol(
-                x - MPCD( ub(i) ), options_.solfeastol_);
+            x - MPCD( ub(i) ),
+            options_.solfeastol_,
+            GetModel().var_name(i));
       if (is_var_integer(i))
         chk.VarViolIntty().at(aux).CheckViol(
-                  std::fabs(x - std::round(x)),
-                  options_.solinttol_);
+              std::fabs(x - std::round(x)),
+              options_.solinttol_,
+              GetModel().var_name(i));
     }
   }
 
@@ -649,10 +655,68 @@ class FlatConverter :
     GetModel().ComputeViolations(chk);
   }
 
-  void CheckObjs(SolCheck& ) {
+  void CheckObjs(SolCheck& ) { }
+
+  void GenerateViolationsReport(SolCheck& chk) {
+    fmt::MemoryWriter wrt;
+    if (chk.HasAnyConViols()) {
+      wrt.write(
+            "Constraint violations "
+            "(sol:chk:feastol={}, sol:chk:inttol={}):\n",
+            options_.solfeastol_, options_.solinttol_);
+      Gen1Viol(chk.VarViolBnds().at(0), wrt,
+               "  - {} original variable(s) violate bounds, max by {}");
+      Gen1Viol(chk.VarViolBnds().at(1), wrt,
+               "  - {} auxiliary variable(s) violate bounds, max by {}");
+      Gen1Viol(chk.VarViolIntty().at(0), wrt,
+               "  - {} original variable(s) violate integrality, max by {}");
+      Gen1Viol(chk.VarViolIntty().at(1), wrt,
+               "  - {} auxiliary variable(s) violate integrality, max by {}");
+    }
+    GenConViol(chk.ConViolAlg(), wrt, "Algebraic");
+    GenConViol(chk.ConViolLog(), wrt, "Logical");
+    if (chk.HasAnyObjViols()) {
+      wrt.write("Objective value violations"
+                "(sol:chk:feastol={})\n",
+                options_.solfeastol_);
+      Gen1Viol(chk.ObjViols(), wrt,
+               "  - {} objective value(s) violated, max by {}");
+    }
+    chk.SetReport( wrt.str() );
+  }
+
+  void Gen1Viol(
+      const ViolSummary& vs, fmt::MemoryWriter& wrt,
+      const std::string& format) {
+    if (vs.N_) {
+      wrt.write(format, vs.N_, vs.epsMax_);
+      if (vs.name_ && *vs.name_ != '\0')
+        wrt.write(" (item '{}')", vs.name_);
+      wrt.write("\n");
+    }
   }
 
-  void GenerateViolationsReport(SolCheck& chk) { }
+  void GenConViol(
+      const std::map< std::string, ViolSummArray<3> >& cvmap,
+      fmt::MemoryWriter& wrt, const std::string& classnm) {
+    if (cvmap.size()) {
+      wrt.write(classnm + " constraint violations:\n");
+      for (const auto& cva: cvmap) {
+        Gen1Viol(cva.second.at(0), wrt,
+                 "  - {} original constraint(s) of type '"
+                 + std::string(cva.first)
+                 + "' violate bounds, max by {}");
+        Gen1Viol(cva.second.at(1), wrt,
+                 "  - {} reformulated constraint(s) of type '"
+                 + std::string(cva.first)
+                 + "' violate bounds, max by {}");
+        Gen1Viol(cva.second.at(2), wrt,
+                 "  - {} auxiliary constraint(s) of type '"
+                 + std::string(cva.first)
+                 + "' violate bounds, max by {}");
+      }
+    }
+  }
 
   //////////////////////////// UTILITIES /////////////////////////////////
   ///
@@ -961,7 +1025,7 @@ class FlatConverter :
 
     bool solcheckfail_ = false;
     double solfeastol_ = 1e-6;
-    double solinttol_ = 1e-6;
+    double solinttol_ = 1e-5;
   };
   Options options_;
 

include/mp/flat/converter_model.h

@@ -61,6 +61,12 @@ class FlatModel
     var_names_storage_ = names;
   }
 
+  /// To be used only after names presolve
+  const char* var_name(int i) const {
+    return i<(int)var_names_storage_.size()
+        ? var_names_storage_[i].c_str() : nullptr;
+  }
+
   int num_vars() const
   { assert(check_vars()); return (int)var_lb_.size(); }
 

include/mp/flat/expr_affine.h

@@ -74,10 +74,10 @@ class LinTerms {
   const LinTerms& GetLinTerms() const { return *this; }
 
   /// Compute value given a dense vector of variable values
-  double ComputeValue(ArrayRef<double> x) const {
-    double s=0.0;
+  long double ComputeValue(const ArrayRef<double>& x) const {
+    long double s=0.0;
     for (size_t i=coefs().size(); i--; )
-      s += coefs()[i] * x[vars()[i]];
+      s += (long double)(coefs()[i]) * x[vars()[i]];
     return s;
   }
 

include/mp/flat/expr_quadratic.h

@@ -54,10 +54,10 @@ class QuadTerms {
   int var2(int i) const { return vars2_[i]; }
 
   /// Compute value given a dense vector of variable values
-  double ComputeValue(ArrayRef<double> x) const {
-    double s=0.0;
+  long double ComputeValue(const ArrayRef<double>& x) const {
+    long double s=0.0;
     for (size_t i=coefs().size(); i--; )
-      s += coefs()[i] * x[vars1()[i]] * x[vars2()[i]];
+      s += (long double)(coefs()[i]) * x[vars1()[i]] * x[vars2()[i]];
     return s;
   }
 
@@ -200,7 +200,7 @@ class QuadAndLinTerms :
   }
 
   /// Value at given variable vector
-  double ComputeValue(ArrayRef<double> x) const {
+  long double ComputeValue(const ArrayRef<double>& x) const {
     return LinTerms::ComputeValue(x) + QuadTerms::ComputeValue(x);
   }
 

include/mp/valcvt.h

@@ -296,10 +296,11 @@ class ValuePresolver : public ValuePresolverImpl {
       const auto& mx = mv.GetVarValues();
       const auto& mo = mv.GetObjValues();
       if (mx.IsSingleKey()
-          && mx().size()          // solution available
-          && mo.IsSingleKey()) {
-        solchk_(mx(),
-                mv.GetConValues(), mo());
+          && mx().size()) {          // solution available
+        ArrayRef<double> objs;
+        if (mo.IsSingleKey())
+          objs = mo();
+        solchk_(mx(), mv.GetConValues(), objs);
       }
     }
     return ValuePresolverImpl::PostsolveSolution(mv);

test/end2end/cases/categorized/fast/tech/infeas_int_01.mod

@@ -0,0 +1,21 @@
+
+# -------------------------------------------------------------
+# IIS Int 01
+# infeas_int_01.mod
+# -------------------------------------------------------------
+
+var x integer;
+var y integer;
+
+minimize TotalSum:
+    x - 2*y;
+
+subj to C1:
+       -x + 21*y >=  2;
+
+subj to C2:
+     -3*x +  2*y <=  1;
+
+subj to C3:
+     20*x +    y <= 20;
+

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

@@ -25,5 +25,17 @@
     "files" : ["diet.mod", "diet.dat"],
     "options": { "ANYSOLVER_options": "outlev 1 barrier outlev 0" },
     "objective": 88.2
+  },
+  {
+    "name" : "Check sol:chk:fail",
+    "tags" : ["linear", "feasrelax"],
+    "files" : ["infeas_int_01.mod"],
+    "options": {
+      "ANYSOLVER_options":
+        "sol:chk:fail sol:chk:feastol=1e-7 sol:chk:inttol=1e-4 feasrelax=1"
+    },
+    "values": {
+      "solve_result_num": 520
+    }
   }
 ]