SolCheck: relative tol #200

include/mp/flat/constr_algebraic.h

@@ -75,16 +75,23 @@ class AlgebraicConstraint :
   /// report violation amount
   /// (negative if holds with slack.)
   /// In logical mode, report 1/0
-  /// (what's the violation amount of 0 for >0?)
+  /// (what's the violation amount
+  ///  for "x>0" when x==0?)
   template <class VarInfo>
-  double
+  Violation
   ComputeViolation(const VarInfo& x, bool logical=false) const {
     double bd = Body::ComputeValue(x);
     if (!logical) {
-      return std::max(      // same for strict cmp?
-          RhsOrRange::lb() - bd, bd - RhsOrRange::ub());
+      if (RhsOrRange::lb() > bd)
+        return {RhsOrRange::lb() - bd, RhsOrRange::lb()};
+      if (bd > RhsOrRange::ub())
+        return {bd - RhsOrRange::ub(), RhsOrRange::ub()};
+      return
+      {std::max( // negative. Same for strict cmp?
+                 RhsOrRange::lb() - bd, bd - RhsOrRange::ub()),
+            0.0};
     }
-    return !RhsOrRange::is_valid(bd);
+    return {double(!RhsOrRange::is_valid(bd)), 1.0};
   }
 
   /// Sorting and merging terms, some solvers require

include/mp/flat/constr_base.h

@@ -16,6 +16,25 @@
 
 namespace mp {
 
+/// Constraint/obj violation
+struct Violation {
+  double viol_;    // abs violation: >0 if really violated
+  double valX_;    // value compared to
+  /// Compute whether violated + relative violation.
+  /// Both absolute and relative should be violated
+  /// (relative only if refvar!=0.)
+  std::pair<bool, double> Check(
+      double epsabs, double epsrel) const {
+    double violRel {0.0};
+    if (viol_ > epsabs
+        && (0.0==std::fabs(valX_)
+            || (violRel=std::fabs(viol_/valX_))>epsrel))
+      return {true, violRel};
+    return {false, 0.0};
+  }
+};
+
+
 /// Custom constraints to derive from, so that overloaded default settings work
 class BasicConstraint {
 public:
@@ -40,8 +59,8 @@ class BasicConstraint {
   int GetResultVar() const { return -1; }
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolation(const VarInfo& ) const
-  { return 0.0; }
+  Violation ComputeViolation(const VarInfo& ) const
+  { return {0.0, 0.0}; }
 
 private:
   std::string name_;
@@ -159,7 +178,7 @@ class CustomFunctionalConstraint :
 
   /// Compute violation
   template <class VarVec>
-  double ComputeViolation(const VarVec& x) const;
+  Violation ComputeViolation(const VarVec& x) const;
 };
 
 
@@ -178,26 +197,26 @@ double ComputeValue(const Con& , const VarVec& ) {
 /// should be redefined for cones.
 template <class Args, class Params,
           class NumOrLogic, class Id, class VarVec>
-double ComputeViolation(
+Violation ComputeViolation(
     const CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>& c,
     const VarVec& x) {
   auto resvar = c.GetResultVar();
   if (!x.recomp_vals()) {    // solver's var values: normal check
     auto viol = x[resvar] - ComputeValue(c, x);
     switch (c.GetContext().GetValue()) {
     case Context::CTX_MIX:
-      return std::fabs(viol);
+      return {std::fabs(viol), x[resvar]};
     case Context::CTX_POS:
-      return viol;
+      return {viol, x[resvar]};
     case Context::CTX_NEG:
-      return -viol;
+      return {-viol, x[resvar]};
     default:
-      return INFINITY;
+      return {INFINITY, 0.0};
     }
   }
   return                              // recomputed var minus solver's
-      std::fabs(x[resvar] - x.raw(resvar))
-      + std::max(0.0, x.bounds_viol(resvar));
+  { std::fabs(x[resvar] - x.raw(resvar))
+        + std::max(0.0, x.bounds_viol(resvar)), x[resvar]};
 }
 
 
@@ -294,19 +313,25 @@ class ConditionalConstraint :
   /// If the subconstr holds but should not,
   /// return the opposite gap. Similar to Indicator.
   template <class VarVec>
-  double ComputeViolation(const VarVec& x) {
+  Violation ComputeViolation(const VarVec& x) {
     auto viol = GetConstraint().ComputeViolation(x);
-    bool ccon_valid = viol<=0.0;
+    bool ccon_valid = viol.viol_<=0.0;
     bool has_arg = x[GetResultVar()] >= 0.5;
     switch (this->GetContext().GetValue()) {
     case Context::CTX_MIX:    // Viol is non-positive if holds
-      return has_arg == ccon_valid ? 0.0 : std::fabs(viol);
+      if (has_arg == ccon_valid)
+        return {0.0, 0.0};
+      return {std::fabs(viol.viol_), viol.valX_};
     case Context::CTX_POS:
-      return has_arg <= ccon_valid ? 0.0 : viol;
+      if (has_arg <= ccon_valid)
+        return {0.0, 0.0};
+      return {viol.viol_, viol.valX_};
     case Context::CTX_NEG:
-      return has_arg >= ccon_valid ? 0.0 : -viol;
+      if (has_arg >= ccon_valid)
+        return {0.0, 0.0};
+      return {-viol.viol_, viol.valX_};
     default:
-      return INFINITY;
+      return {INFINITY, 0.0};
     }
   }
 };

include/mp/flat/constr_eval.h

@@ -282,49 +282,52 @@ template <class Con, class VarVec>
 double ComputeValue(
     const ConditionalConstraint<Con>& con, const VarVec& x) {
   auto viol = con.GetConstraint().ComputeViolation(x, true);
-  bool ccon_valid = viol<=0.0;
+  bool ccon_valid = viol.viol_<=0.0;
   return ccon_valid;
 }
 
 /// Compute violation of the QuadraticCone constraint.
 template <class VarVec>
-double ComputeViolation(
+Violation ComputeViolation(
     const QuadraticConeConstraint& con, const VarVec& x) {
   const auto& args = con.GetArguments();
   const auto& params = con.GetParameters();
   assert(args.size()==params.size());
   double sum = 0.0;
   for (auto i=args.size(); --i; )
     sum += std::pow( params[i]*x[args[i]], 2.0 );
-  return std::sqrt(sum) - params[0]*x[args[0]];
+  return {std::sqrt(sum) - params[0]*x[args[0]],
+        sum};
 }
 
 /// Compute violation of the RotatedQuadraticCone constraint.
 template <class VarVec>
-double ComputeViolation(
+Violation ComputeViolation(
     const RotatedQuadraticConeConstraint& con, const VarVec& x) {
   const auto& args = con.GetArguments();
   const auto& params = con.GetParameters();
   assert(args.size()==params.size());
   double sum = 0.0;
   for (auto i=args.size(); --i>1; )
     sum += std::pow( params[i]*x[args[i]], 2.0 );
-  return sum
-      - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]];
+  return {sum
+        - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]],
+        sum};
 }
 
 /// Compute violation of the ExponentialCone constraint.
 template <class VarVec>        // ax >= by exp(cz / (by))
-double ComputeViolation(       // where ax, by >= 0
+Violation ComputeViolation(       // where ax, by >= 0
     const ExponentialConeConstraint& con, const VarVec& x) {
   const auto& args = con.GetArguments();
   const auto& params = con.GetParameters();
   auto ax = params[0]*x[args[0]];
   auto by = params[1]*x[args[1]];
   if (0.0==std::fabs(by))
-    return -ax;
+    return {-ax, 0.0};
   auto cz = params[2]*x[args[2]];
-  return by * std::exp(cz / by) - ax;
+  return {by * std::exp(cz / by) - ax,
+        by * std::exp(cz / by)};
 }
 
 /// Compute result of the PL constraint.
@@ -354,7 +357,7 @@ double ComputeValue(const PLConstraint& con, const VarVec& x) {
 template <class Args, class Params,
           class NumOrLogic, class Id>
 template <class VarVec>
-double
+Violation
 CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>
 ::ComputeViolation(const VarVec& x) const
 { return mp::ComputeViolation(*this, x); }

include/mp/flat/constr_general.h

@@ -44,12 +44,12 @@ class IndicatorConstraint: public BasicConstraint {
 
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolation(const VarInfo& x) const {
+  Violation ComputeViolation(const VarInfo& 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;
+    return {0.0, 0.0};
   }
 
 
@@ -151,26 +151,26 @@ class SOS_1or2_Constraint: public BasicConstraint {
 
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolation(const VarInfo& x) const {
+  Violation ComputeViolation(const VarInfo& x) const {
     if (1==type)
       return ComputeViolationSOS1(x);
     return ComputeViolationSOS2(x);
   }
 
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolationSOS1(const VarInfo& x) const {
+  Violation ComputeViolationSOS1(const VarInfo& x) const {
     int nnz=0;
     for (int i=(int)v_.size(); i--; ) {
       if (x.is_nonzero(v_[i]))
         ++nnz;
     }
-    return std::max(0, nnz-1);
+    return {(double)std::max(0, nnz-1), 0.0};
   }
 
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolationSOS2(const VarInfo& x) const {
+  Violation ComputeViolationSOS2(const VarInfo& x) const {
     int npos=0, pos1=-1, pos2=-1;
     int posDist=1;     // should be 1
     for (int i=(int)v_.size(); i--; ) {
@@ -184,8 +184,8 @@ class SOS_1or2_Constraint: public BasicConstraint {
         }
       }
     }
-    return std::max(0, npos-2)
-        + std::abs(1 - posDist);
+    return {double(std::max(0, npos-2)
+        + std::abs(1 - posDist)), 0.0};
   }
 
 
@@ -247,13 +247,13 @@ class ComplementarityConstraint : public BasicConstraint {
 
   /// Compute violation
   template <class VarInfo>
-  double ComputeViolation(const VarInfo& x) const {
+  Violation ComputeViolation(const VarInfo& x) const {
     auto ve = compl_expr_.ComputeValue(x);
     if (x.is_at_lb(compl_var_))
-      return -ve;
+      return {-ve, 0.0};
     else if (x.is_at_ub(compl_var_))
-      return ve;
-    return std::fabs(ve);
+      return {ve, 0.0};
+    return {std::fabs(ve), 0.0};
   }
 
 

include/mp/flat/constr_keeper.h

@@ -32,24 +32,36 @@ 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
+  /// Check if this violation should be counted.
   void CheckViol(
-      double val, double eps, const char* nm) {
-    if (val > eps)
-      CountViol(val, nm);
+      Violation viol,
+      double epsabs, double epsrel,
+      const char* nm) {
+    auto chk = viol.Check(epsabs, epsrel);
+    if (chk.first)
+      CountViol(viol, chk.second, nm);
   }
   /// Count violation
-  void CountViol(double val, const char* nm) {
+  void CountViol(
+      Violation viol, double violRel, const char* nm) {
     ++N_;
-    if (epsMax_ < val)
-      epsMax_ = val;
-    name_ = nm;
+    if (epsAbsMax_ < viol.viol_) {
+      epsAbsMax_ = viol.viol_;
+      nameAbs_ = nm;
+    }
+    if (epsRelMax_ < violRel) {
+      epsRelMax_ = violRel;
+      nameRel_ = nm;
+    }
   }
   int N_ {0};
-  double epsMax_ {0.0};
-  const char* name_ {nullptr};
+  double epsAbsMax_ {0.0};
+  const char* nameAbs_ {nullptr};
+  double epsRelMax_ {0.0};
+  const char* nameRel_ {nullptr};
 };
 
+
 /// Array of violation summaries.
 /// For different kinds, e.g., original / aux vars.
 template <int Nkinds>
@@ -245,13 +257,13 @@ struct SolCheck {
            ArrayRef<double> x_raw,
            ArrayRef<var::Type> vtype,
            ArrayRef<double> lb,  ArrayRef<double> ub,
-           double feastol, double , //inttol,
+           double feastol, double feastolrel,
            const char* sol_rnd, const char* sol_prec,
            bool recomp_vals, int chk_mode)
     : x_(feastol, recomp_vals,
          x, x_raw, vtype, lb, ub, sol_rnd, sol_prec),
       obj_(obj),
-      feastol_(feastol),
+      feastol_(feastol), feastolrel_(feastolrel),
       fRecomputedVals_(recomp_vals),
       check_mode_(chk_mode) { }
   /// Any violations?
@@ -279,8 +291,10 @@ struct SolCheck {
   const ArrayRef<double>& obj_vals() const
   { return obj_; }
 
-  /// Feasibility tolerance
+  /// Absolute feasibility tolerance
   double GetFeasTol() const { return feastol_; }
+  /// Relative feasibility tolerance
+  double GetFeasTolRel() const { return feastolrel_; }
 
   /// Using recomputed aux vars?
   bool if_recomputed() const { return fRecomputedVals_; }
@@ -313,6 +327,7 @@ struct SolCheck {
   VarInfoStatic x_;
   ArrayRef<double> obj_;
   double feastol_;
+  double feastolrel_;
   bool fRecomputedVals_;
   int check_mode_;
 
@@ -906,7 +921,9 @@ class ConstraintKeeper final
             c_class = 4;      // intermediate
           if (c_class & chk.check_mode()) {
             auto viol = cons_[i].con_.ComputeViolation(x);
-            if (viol > chk.GetFeasTol()) {
+            auto cr = viol.Check(
+                  chk.GetFeasTol(), chk.GetFeasTolRel());
+            if (cr.first) {
               if (!conviolarray)
                 conviolarray =         // lazy map access
                     &conviolmap[GetShortTypeName()];
@@ -917,7 +934,7 @@ class ConstraintKeeper final
                                           ? 2 : 1;
               assert(index < (int)conviolarray->size());
               (*conviolarray)[index].CountViol(
-                    viol, cons_[i].con_.name());
+                    viol, cr.second, cons_[i].con_.name());
             }
           }
         }