Check Complementarity #200

include/mp/flat/constr_general.h

@@ -244,6 +244,17 @@ class ComplementarityConstraint : public BasicConstraint {
   /// Get variable
   int GetVariable() const { return compl_var_; }
 
+  /// Compute violation
+  template <class VarInfo>
+  double ComputeViolation(const VarInfo& x) const {
+    auto ve = compl_expr_.ComputeValue(x);
+    if (x.is_at_lb(compl_var_))
+      return -ve;
+    else if (x.is_at_ub(compl_var_))
+      return ve;
+    return std::fabs(ve);
+  }
+
 
 private:
   ExprType compl_expr_;

include/mp/flat/constr_keeper.h

@@ -61,9 +61,14 @@ class VarInfo {
   /// Constructor
   VarInfo(double ft,
           std::vector<double> x,
-          ArrayRef<var::Type> type)
+          ArrayRef<var::Type> type,
+          ArrayRef<double> lb, ArrayRef<double> ub)
     : feastol_(ft), x_(std::move(x)), x_ref_(x_),
-      type_(type) { }
+      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());
+  }
   /// Access variable value
   double operator[]( int i ) const {
     assert(i>=0 && i<(int)x_.size());
@@ -89,6 +94,13 @@ class VarInfo {
         (*this)[i]
         >= (is_var_int(i) ? 0.5 : feastol());
   }
+  /// Is at lb?
+  bool is_at_lb(int i) const
+  { return (*this)[i] - lb_[i] <= feastol(); }
+  /// Is at ub?
+  bool is_at_ub(int i) const
+  { return -((*this)[i] - ub_[i]) <= feastol(); }
+
   /// Feasibility tolerance
   double feastol() const { return feastol_; }
   /// x() as ArrayRef
@@ -100,6 +112,8 @@ class VarInfo {
   const 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_;
 };
 
 
@@ -110,8 +124,9 @@ struct SolCheck {
            const pre::ValueMapDbl& duals,
            ArrayRef<double> obj,
            ArrayRef<var::Type> vtype,
+           ArrayRef<double> lb,  ArrayRef<double> ub,
            double feastol, double inttol)
-    : x_(feastol, x, vtype), y_(duals), obj_(obj),
+    : x_(feastol, x, vtype, lb, ub), y_(duals), obj_(obj),
       feastol_(feastol), inttol_(inttol) { }
   /// Any violations?
   bool HasAnyViols() const

include/mp/flat/converter.h

@@ -613,7 +613,10 @@ class FlatConverter :
       ArrayRef<double> x,
       const pre::ValueMapDbl& duals,
       ArrayRef<double> obj) {
-    SolCheck chk(x, duals, obj, GetModel().var_type_vec(),
+    SolCheck chk(x, duals, obj,
+                 GetModel().var_type_vec(),
+                 GetModel().var_lb_vec(),
+                 GetModel().var_ub_vec(),
                  options_.solfeastol_, options_.solinttol_);
     CheckVars(chk);
     CheckCons(chk);

include/mp/flat/converter_model.h

@@ -70,6 +70,7 @@ class FlatModel
   int num_vars() const
   { assert(check_vars()); return (int)var_lb_.size(); }
 
+  /// Not auxiliary?
   bool is_var_original(int i) const
   { return i<num_vars_orig_; }
 
@@ -88,6 +89,10 @@ class FlatModel
     return var_type_[v];
   }
 
+  const std::vector<double>& var_lb_vec() const
+  { return var_lb_; }
+  const std::vector<double>& var_ub_vec() const
+  { return var_ub_; }
   const std::vector<var::Type>& var_type_vec() const
   { return var_type_; }
 

include/mp/flat/expr_affine.h

@@ -158,6 +158,6 @@ class LinTerms {
 /// Typedef AffineExpr
 using AffineExpr = AlgebraicExpression<LinTerms>;
 
-} // namepsace mp
+} // namespace mp
 
 #endif // AFFINE_EXPR_H

include/mp/flat/expr_algebraic.h

@@ -9,6 +9,8 @@
 #include <string>
 #include <cmath>
 
+#include "mp/arrayref.h"
+
 namespace mp {
 
 /// A template algebraic expression:
@@ -73,6 +75,10 @@ class AlgebraicExpression : public Body {
   /// Add to constant
   void add_to_constant(double a) { constant_term_ += a; }
 
+  /// Compute value given a dense vector of variable values
+  double ComputeValue(const ArrayRef<double>& x) const
+  { return Body::ComputeValue(x) + constant_term(); }
+
   /// Negate the ae
   void negate() {
     Body::negate();