Check SOS #200

ampl · Aug 28, 2023 · bab5d21 · bab5d21
1 parent 5252a55
commit bab5d21
Show file tree

Hide file tree

Showing 5 changed files with 101 additions and 7 deletions.
diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h
@@ -186,7 +186,6 @@ double ComputeViolation(
     const CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>& c,
     const VarVec& x) {
   auto resvar = c.GetResultVar();
-  assert(resvar < (int)x.size());
   return std::fabs(x[resvar] - ComputeValue(c, x));
 }
 

diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h
@@ -144,9 +144,49 @@ class SOS_1or2_Constraint: public BasicConstraint {
     sort();
     assert(check());
   }
+  /// Check data
   bool check() const { return type>=1 && type<=2 &&
                      v_.size()==w_.size(); }
 
+  /// Compute violation
+  template <class VarInfo>
+  double 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 {
+    int nnz=0;
+    for (int i=(int)v_.size(); i--; ) {
+      if (x.is_nonzero(v_[i]))
+        ++nnz;
+    }
+    return std::max(0, nnz-1);
+  }
+
+  /// Compute violation
+  template <class VarInfo>
+  double 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--; ) {
+      if (x.is_positive(v_[i])) {
+        ++npos;
+        if (pos1<0)
+          pos1=i;
+        else if (pos2<0) {
+          pos2=i;
+          posDist = pos1-pos2;
+        }
+      }
+    }
+    return std::max(0, npos-2)
+        + std::abs(1 - posDist);
+  }
+
 
 protected:
   /// Sort by weights

diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h
@@ -54,14 +54,64 @@ struct ViolSummary {
 template <int Nkinds>
 using ViolSummArray = std::array<ViolSummary, Nkinds>;
 
+
+/// Variable information used by solution check
+class VarInfo {
+public:
+  /// Constructor
+  VarInfo(double ft,
+          std::vector<double> x,
+          ArrayRef<var::Type> type)
+    : feastol_(ft), x_(std::move(x)), x_ref_(x_),
+      type_(type) { }
+  /// Access variable value
+  double operator[]( int i ) const {
+    assert(i>=0 && i<(int)x_.size());
+    return x_[i];
+  }
+  /// Access whole vectorref
+  operator const ArrayRef<double>& () const
+  { return x_ref(); }
+  /// Access integrality condition
+  bool is_var_int(int i) const {
+    assert(i>=0 && i<(int)type_.size());
+    return var::INTEGER==type_[i];
+  }
+  /// Variable value nonzero?
+  bool is_nonzero(int i) const {
+    return
+        std::fabs( (*this)[i] )
+        >= (is_var_int(i) ? 0.5 : feastol());
+  }
+  /// Variable value positive?
+  bool is_positive(int i) const {
+    return
+        (*this)[i]
+        >= (is_var_int(i) ? 0.5 : feastol());
+  }
+  /// Feasibility tolerance
+  double feastol() const { return feastol_; }
+  /// x() as ArrayRef
+  const ArrayRef<double>& x_ref() const { return x_ref_; }
+
+private:
+  double feastol_;
+
+  const std::vector<double> x_;   // can be rounded, etc.
+  const ArrayRef<double> x_ref_;
+  const ArrayRef<var::Type> type_;
+};
+
+
 /// Solution check data
 struct SolCheck {
   /// Construct
   SolCheck(ArrayRef<double> x,
            const pre::ValueMapDbl& duals,
            ArrayRef<double> obj,
+           ArrayRef<var::Type> vtype,
            double feastol, double inttol)
-    : x_(x), y_(duals), obj_(obj),
+    : x_(feastol, x, vtype), y_(duals), obj_(obj),
       feastol_(feastol), inttol_(inttol) { }
   /// Any violations?
   bool HasAnyViols() const
@@ -80,8 +130,10 @@ struct SolCheck {
   /// Summary
   const std::string& GetReport() const { return report_; }
 
-  /// x
-  ArrayRef<double>& x() { return x_; }
+  /// VarInfo, can be used like x() for templates
+  const VarInfo& x_ext() const { return x_; }
+  /// x as array
+  const ArrayRef<double>& x() { return x_.x_ref(); }
   /// x[i]
   double x(int i) const { return x_[i]; }
   /// Feasibility tolerance
@@ -110,7 +162,7 @@ struct SolCheck {
   { report_ = std::move(rep); }
 
 private:
-  ArrayRef<double> x_;
+  VarInfo x_;
   const pre::ValueMapDbl& y_;
   ArrayRef<double> obj_;
   double feastol_;
@@ -667,7 +719,7 @@ class ConstraintKeeper : public BasicConstraintKeeper {
           cons_.front().con_.IsLogical() ?
             chk.ConViolLog() :
             chk.ConViolAlg();
-      const auto& x = chk.x();
+      const auto& x = chk.x_ext();
       ViolSummArray<3>* conviolarray {nullptr};
       for (int i=(int)cons_.size(); i--; ) {
         if (!cons_[i].IsUnused()) {

diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h
@@ -613,7 +613,7 @@ class FlatConverter :
       ArrayRef<double> x,
       const pre::ValueMapDbl& duals,
       ArrayRef<double> obj) {
-    SolCheck chk(x, duals, obj,
+    SolCheck chk(x, duals, obj, GetModel().var_type_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
@@ -88,6 +88,9 @@ class FlatModel
     return var_type_[v];
   }
 
+  const std::vector<var::Type>& var_type_vec() const
+  { return var_type_; }
+
   template <class VarArray>
   double lb_array(const VarArray& va) const {
     double result = Inf();