From efba5c48b7080a78c426738d046c01f4873fe434 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Wed, 23 Aug 2023 12:59:38 +1000 Subject: [PATCH] INFINITY (vs numeric_limits::infinity()) Should be more portable Start ComputeValue()'s for functional constraints #200 --- include/mp/backend-std.h | 5 +- include/mp/common.h | 5 +- include/mp/flat/constr_base.h | 45 ++++++++++++++++++ include/mp/flat/constr_keeper.h | 2 + include/mp/flat/constr_viol.h | 72 +++++++++++++++++++++++++++++ include/mp/flat/preprocess.h | 1 - include/mp/flat/redef/conic/cones.h | 1 - include/mp/nl-reader.h | 4 +- solvers/mosek/mosekmodelapi.cc | 18 ++++---- src/expr-writer.h | 4 +- test/common-test.cc | 2 +- test/function-test.cc | 3 +- test/nl-reader-test.cc | 5 +- test/solvers/nl-solver-test.h | 6 +-- 14 files changed, 147 insertions(+), 26 deletions(-) create mode 100644 include/mp/flat/constr_viol.h diff --git a/include/mp/backend-std.h b/include/mp/backend-std.h index 18c211b2c..fe37586b5 100644 --- a/include/mp/backend-std.h +++ b/include/mp/backend-std.h @@ -22,7 +22,6 @@ #define STD_BACKEND_H_ #include -#include #include #include "mp/clock.h" @@ -584,10 +583,10 @@ class StdBackend : /// AMPL's inf static constexpr double AMPLInf() - { return std::numeric_limits::infinity(); } + { return INFINITY; } /// AMPL's -inf static constexpr double AMPLMinusInf() - { return -std::numeric_limits::infinity(); } + { return -INFINITY; } public: diff --git a/include/mp/common.h b/include/mp/common.h index 8f9ed0fbb..a251ca77b 100644 --- a/include/mp/common.h +++ b/include/mp/common.h @@ -23,6 +23,7 @@ #ifndef MP_COMMON_H_ #define MP_COMMON_H_ +#include #include // for std::size_t #include "mp/error.h" // for MP_ASSERT @@ -86,13 +87,13 @@ class ComplInfo { /** Constraint lower bound. */ double con_lb() const { return (flags_ & INF_LB) != 0 ? - -std::numeric_limits::infinity() : 0; + -INFINITY : 0; } /** Constraint upper bound. */ double con_ub() const { return (flags_ & INF_UB) != 0 ? - std::numeric_limits::infinity() : 0; + INFINITY : 0; } }; diff --git a/include/mp/flat/constr_base.h b/include/mp/flat/constr_base.h index 42caf8897..0d7216ef7 100644 --- a/include/mp/flat/constr_base.h +++ b/include/mp/flat/constr_base.h @@ -8,6 +8,9 @@ #include #include +#include "mp/format.h" +#include "mp/error.h" + #include "mp/flat/context.h" #include "mp/arrayref.h" @@ -104,6 +107,7 @@ using DblParamArray3 = ParamArrayN; /// Variable-length parameter array using DblParamArray = std::vector; + /// A functional constraint with given arguments /// and further info as parameters /// @param Args: arguments type @@ -151,9 +155,50 @@ class CustomFunctionalConstraint : const Parameters& GetParameters() const { return params_; } /// Get Parameters& Parameters& GetParameters() { return params_; } + + /// Compute violation + template + double ComputeViolation(const VarVec& x) const; }; +/// Dummy template: compute result of functional constraint. +/// Should be implemented for proper functional specializations, +/// but has no sense for conic constraints, for example. +template +double ComputeValue( + const CustomFunctionalConstraint& , + const VarVec& ) { + MP_RAISE(fmt::format("ComputeValue({}) not implemented.", + typeid( + CustomFunctionalConstraint).name())); + return 0.0; +} + +/// Compute violation of functional constraint. +/// Can only be used for proper functional constraints, +/// should be redefined for cones. +template +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)); +} + +template +template +double +CustomFunctionalConstraint +::ComputeViolation(const VarVec& x) const +{ return mp::ComputeViolation(*this, x); } + + /// A base class for numerical functional constraint. /// It provides default properties of such a constraint class NumericFunctionalConstraintTraits { diff --git a/include/mp/flat/constr_keeper.h b/include/mp/flat/constr_keeper.h index eda4ceed6..5cc1fd854 100644 --- a/include/mp/flat/constr_keeper.h +++ b/include/mp/flat/constr_keeper.h @@ -9,10 +9,12 @@ #include "mp/common.h" #include "mp/format.h" #include "mp/env.h" + #include "mp/flat/model_api_base.h" #include "mp/flat/constr_hash.h" #include "mp/flat/redef/redef_base.h" #include "mp/valcvt-node.h" +#include "mp/flat/constr_viol.h" namespace mp { diff --git a/include/mp/flat/constr_viol.h b/include/mp/flat/constr_viol.h new file mode 100644 index 000000000..34279116d --- /dev/null +++ b/include/mp/flat/constr_viol.h @@ -0,0 +1,72 @@ +#ifndef CONSTR_VIOL_H +#define CONSTR_VIOL_H + +/** + * Violations of (mainly functional) constraints + */ + +#include + +#include "constr_functional.h" +#include "constr_general.h" + +namespace mp { + +/// Compute result of the max constraint. +template +double ComputeValue(const MaxConstraint& con, const VarVec& x) { + double result = -INFINITY; + for (auto i: con.GetArguments()) { + if (result < x[i]) + result = x[i]; + } + return result; +} + +/// Compute result of the min constraint. +template +double ComputeValue(const MinConstraint& con, const VarVec& x) { + double result = INFINITY; + for (auto i: con.GetArguments()) { + if (result > x[i]) + result = x[i]; + } + return result; +} + +/// Compute result of the abs constraint. +template +double ComputeValue(const AbsConstraint& con, const VarVec& x) { + return std::fabs(x[con.GetArguments()[0]]); +} + +/// Compute result of the and constraint. +template +double ComputeValue(const AndConstraint& con, const VarVec& x) { + for (auto i: con.GetArguments()) { + if (x[i] < 0.5) + return 0.0; + } + return 1.0; +} + +/// Compute result of the or constraint. +template +double ComputeValue(const OrConstraint& con, const VarVec& x) { + for (auto i: con.GetArguments()) { + if (x[i] >= 0.5) + return 1.0; + } + return 0.0; +} + +/// Compute result of the not constraint. +template +double ComputeValue(const NotConstraint& con, const VarVec& x) { + return x[con.GetArguments()[0]] < 0.5; +} + + +} // namespace mp + +#endif // CONSTR_VIOL_H diff --git a/include/mp/flat/preprocess.h b/include/mp/flat/preprocess.h index e6e27d754..5b267ef45 100644 --- a/include/mp/flat/preprocess.h +++ b/include/mp/flat/preprocess.h @@ -1,7 +1,6 @@ #ifndef PREPROCESS_H #define PREPROCESS_H -#include #include #include diff --git a/include/mp/flat/redef/conic/cones.h b/include/mp/flat/redef/conic/cones.h index 7f5140e8a..acd20702c 100644 --- a/include/mp/flat/redef/conic/cones.h +++ b/include/mp/flat/redef/conic/cones.h @@ -8,7 +8,6 @@ #define MP_FLAT_CVT_CONES_H #include -#include #include #include diff --git a/include/mp/nl-reader.h b/include/mp/nl-reader.h index ea2b90834..f75d5d3d3 100644 --- a/include/mp/nl-reader.h +++ b/include/mp/nl-reader.h @@ -48,7 +48,7 @@ #include #include #include -#include +#include #include namespace mp { @@ -1823,7 +1823,7 @@ void NLReader::ReadBounds() { double lb = 0, ub = 0; BoundHandler bh(*this); int num_bounds = bh.num_items(); - double infinity = std::numeric_limits::infinity(); + double infinity = INFINITY; for (int i = 0; i < num_bounds; ++i) { switch (reader_.ReadChar() - '0') { case RANGE: diff --git a/solvers/mosek/mosekmodelapi.cc b/solvers/mosek/mosekmodelapi.cc index 99fa5803b..1bed24698 100644 --- a/solvers/mosek/mosekmodelapi.cc +++ b/solvers/mosek/mosekmodelapi.cc @@ -1,3 +1,5 @@ +#include + #include "mosekmodelapi.h" @@ -29,8 +31,8 @@ void MosekModelAPI::AddVariables(const VarArrayDef& v) { MOSEK_CCALL(MSK_appendvars(lp(), v.size())); for (size_t i = v.size(); i--; ) { MSKboundkeye k; - bool freelb = v.plb()[i] == -std::numeric_limits::infinity(); - bool freeub = v.pub()[i] == std::numeric_limits::infinity(); + bool freelb = v.plb()[i] == -INFINITY; + bool freeub = v.pub()[i] == INFINITY; if (freelb && freeub) k = MSK_BK_FR; else if (freelb == freeub) @@ -119,16 +121,16 @@ void MosekModelAPI::AddConstraint(const LinConRange& lc) double ub = lc.ub(); MSKboundkey_enum key; - if (lb == -std::numeric_limits::infinity()) + if (lb == -INFINITY) { - if (ub == std::numeric_limits::infinity()) + if (ub == INFINITY) key = MSK_BK_FR; else key = MSK_BK_UP; } else { - if (ub == std::numeric_limits::infinity()) + if (ub == INFINITY) key = MSK_BK_LO; else if (lb == ub) key = MSK_BK_FX; @@ -178,16 +180,16 @@ void MosekModelAPI::AddConstraint(const QuadConRange& qc) { double ub = qc.ub(); MSKboundkey_enum key; - if (lb == -std::numeric_limits::infinity()) + if (lb == -INFINITY) { - if (ub == std::numeric_limits::infinity()) + if (ub == INFINITY) key = MSK_BK_FR; else key = MSK_BK_UP; } else { - if (ub == std::numeric_limits::infinity()) + if (ub == INFINITY) key = MSK_BK_LO; else if (lb == ub) key = MSK_BK_FX; diff --git a/src/expr-writer.h b/src/expr-writer.h index 0db9a34e9..695a138a8 100644 --- a/src/expr-writer.h +++ b/src/expr-writer.h @@ -23,6 +23,8 @@ #ifndef MP_EXPR_WRITER_H_ #define MP_EXPR_WRITER_H_ +#include + #include "mp/basic-expr-visitor.h" namespace mp { @@ -396,7 +398,7 @@ void WriteExpr(fmt::Writer &w, const LinearExpr &linear, template void Write(fmt::Writer &w, const Problem &p) { // Write variables. - double inf = std::numeric_limits::infinity(); + double inf = INFINITY; int num_vars = p.num_vars(); for (int i = 0; i < num_vars; ++i) { w << "var x" << (i + 1); diff --git a/test/common-test.cc b/test/common-test.cc index 7eeecb293..45b687cc7 100644 --- a/test/common-test.cc +++ b/test/common-test.cc @@ -40,7 +40,7 @@ TEST(CommmonTest, ComplInfo) { EXPECT_NE(0, ComplInfo::INF_LB); EXPECT_NE(0, ComplInfo::INF_UB); EXPECT_NE(ComplInfo::INF_LB, ComplInfo::INF_UB); - double inf = std::numeric_limits::infinity(); + double inf = INFINITY; for (std::size_t i = 0, n = sizeof(flags) / sizeof(*flags); i < n; ++i) { int f = flags[i]; ComplInfo info(f); diff --git a/test/function-test.cc b/test/function-test.cc index 6b3173a3f..d126b8dc9 100644 --- a/test/function-test.cc +++ b/test/function-test.cc @@ -21,7 +21,6 @@ */ #include -#include #include #include #include @@ -565,7 +564,7 @@ TEST(FunctionTest, DifferentiatorDetectsNaN) { Differentiator diff; EXPECT_EQ(0, std::bind2nd(ptr_fun(Hypot), 0)(0)); EXPECT_NE(0, isnan(diff(std::bind2nd(ptr_fun(Hypot), 0), 0))); - EXPECT_EQ(-std::numeric_limits::infinity(), std::log(0.0)); + EXPECT_EQ(-INFINITY, std::log(0.0)); EXPECT_NE(0, isnan(diff(GetDoubleFun(std::log), 0))); } diff --git a/test/nl-reader-test.cc b/test/nl-reader-test.cc index 897593fca..cc38b862c 100644 --- a/test/nl-reader-test.cc +++ b/test/nl-reader-test.cc @@ -25,6 +25,7 @@ // Include the source file to test the implementation. #include "../src/nl-reader.cc" +#include #include #include @@ -488,7 +489,7 @@ class TestNLHandler { void WriteBounds(char type, int index, double lb, double ub) { WriteSep(); - double infinity = std::numeric_limits::infinity(); + double infinity = INFINITY; if (lb != -infinity && lb != ub) log << lb << " <= "; log << type << index; @@ -1922,7 +1923,7 @@ MATCHER_P2(MatchComplInfo, lb, ub, "") { } TEST_F(NLProblemBuilderTest, OnComplementarity) { - double inf = std::numeric_limits::infinity(); + double inf = INFINITY; EXPECT_CALL(builder, SetComplementarity(66, 77, MatchComplInfo(0, inf))); adapter.OnComplementarity(66, 77, mp::ComplInfo(1)); } diff --git a/test/solvers/nl-solver-test.h b/test/solvers/nl-solver-test.h index 623bd723e..dcbc7ca7f 100644 --- a/test/solvers/nl-solver-test.h +++ b/test/solvers/nl-solver-test.h @@ -394,7 +394,7 @@ EvalResult NLSolverTest::Solve( info.num_vars = info.num_nl_integer_vars_in_cons = 4; info.num_logical_cons = 1; pb.SetInfo(info); - double inf = std::numeric_limits::infinity(); + double inf = INFINITY; pb.AddVar(need_result ? -inf : 0, need_result ? inf : 0, var::INTEGER); pb.AddVar(var1, var1, var::INTEGER); pb.AddVar(var2, var2, var::INTEGER); @@ -413,7 +413,7 @@ EvalResult NLSolverTest::Eval( info.num_con_nonzeros = 1; info.num_common_exprs_in_objs = 1; SetInfo(pb, info); - auto inf = std::numeric_limits::infinity(); + auto inf = INFINITY; pb.AddVar(-inf, inf, mp::var::INTEGER); pb.AddVar(var1, var1, mp::var::INTEGER); pb.AddVar(var2, var2, mp::var::INTEGER); @@ -1495,7 +1495,7 @@ TEST_F(NLSolverTest, ZeroUB) { auto info = mp::ProblemInfo(); info.num_vars = info.num_objs = 1; pb.SetInfo(info); - pb.AddVar(-std::numeric_limits::infinity(), 0, var::CONTINUOUS); + pb.AddVar(-INFINITY, 0, var::CONTINUOUS); pb.AddObj(obj::MAX, NumericExpr()); pb.obj(0).set_linear_expr(1).AddTerm(0, 1); EXPECT_EQ(0, Solve(pb).obj_value());