From c60543a9370a657451486704acdd673db567c474 Mon Sep 17 00:00:00 2001
From: Gleb Belov <gleb@ampl.com>
Date: Mon, 28 Aug 2023 21:22:07 +1000
Subject: [PATCH] Check Complementarity #200

---
 include/mp/flat/constr_general.h  | 11 +++++++++++
 include/mp/flat/constr_keeper.h   | 21 ++++++++++++++++++---
 include/mp/flat/converter.h       |  5 ++++-
 include/mp/flat/converter_model.h |  5 +++++
 include/mp/flat/expr_affine.h     |  2 +-
 include/mp/flat/expr_algebraic.h  |  6 ++++++
 6 files changed, 45 insertions(+), 5 deletions(-)

diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h
index 916cfd170..ce14183a9 100644
--- a/include/mp/flat/constr_general.h
+++ b/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_;
diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h
index d1e3453d2..8e36a40a8 100644
--- a/include/mp/flat/constr_keeper.h
+++ b/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
diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h
index 652c9943e..bd9001fa2 100644
--- a/include/mp/flat/converter.h
+++ b/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);
diff --git a/include/mp/flat/converter_model.h b/include/mp/flat/converter_model.h
index 689797067..7c0539fd1 100644
--- a/include/mp/flat/converter_model.h
+++ b/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_; }
 
diff --git a/include/mp/flat/expr_affine.h b/include/mp/flat/expr_affine.h
index 980c4e2af..f5ddb133a 100644
--- a/include/mp/flat/expr_affine.h
+++ b/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
diff --git a/include/mp/flat/expr_algebraic.h b/include/mp/flat/expr_algebraic.h
index cd70a50d0..af7cb681a 100644
--- a/include/mp/flat/expr_algebraic.h
+++ b/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();