From bab5d21246b29ffd7082c72f7c895e1a7b3e0510 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Mon, 28 Aug 2023 17:52:44 +1000 Subject: [PATCH] Check SOS #200 --- include/mp/flat/constr_base.h | 1 - include/mp/flat/constr_general.h | 40 ++++++++++++++++++++ include/mp/flat/constr_keeper.h | 62 ++++++++++++++++++++++++++++--- include/mp/flat/converter.h | 2 +- include/mp/flat/converter_model.h | 3 ++ 5 files changed, 101 insertions(+), 7 deletions(-) diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index 3aa8229ce..36dcc03bc 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -186,7 +186,6 @@ double ComputeViolation( const CustomFunctionalConstraint& 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 index 41b0d6505..916cfd170 100644 --- 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 + double ComputeViolation(const VarInfo& x) const { + if (1==type) + return ComputeViolationSOS1(x); + return ComputeViolationSOS2(x); + } + + /// Compute violation + template + 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 + 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 index 069a7e3c8..d1e3453d2 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -54,14 +54,64 @@ struct ViolSummary { template using ViolSummArray = std::array; + +/// Variable information used by solution check +class VarInfo { +public: + /// Constructor + VarInfo(double ft, + std::vector x, + ArrayRef 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& () 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& x_ref() const { return x_ref_; } + +private: + double feastol_; + + const std::vector x_; // can be rounded, etc. + const ArrayRef x_ref_; + const ArrayRef type_; +}; + + /// Solution check data struct SolCheck { /// Construct SolCheck(ArrayRef x, const pre::ValueMapDbl& duals, ArrayRef obj, + ArrayRef 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& x() { return x_; } + /// VarInfo, can be used like x() for templates + const VarInfo& x_ext() const { return x_; } + /// x as array + const ArrayRef& 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 x_; + VarInfo x_; const pre::ValueMapDbl& y_; ArrayRef 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 index c68393e97..652c9943e 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -613,7 +613,7 @@ class FlatConverter : ArrayRef x, const pre::ValueMapDbl& duals, ArrayRef 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 index e5265e83e..689797067 100644 --- 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_vec() const + { return var_type_; } + template double lb_array(const VarArray& va) const { double result = Inf();