INFINITY (vs numeric_limits<double>::infinity())

Should be more portable

Start ComputeValue()'s for functional constraints #200

include/mp/backend-std.h

@@ -22,7 +22,6 @@
 #define STD_BACKEND_H_
 
 #include <cmath>
-#include <limits>
 #include <functional>
 
 #include "mp/clock.h"
@@ -584,10 +583,10 @@ class StdBackend :
 
   /// AMPL's inf
   static constexpr double AMPLInf()
-  { return std::numeric_limits<double>::infinity(); }
+  { return INFINITY; }
   /// AMPL's -inf
   static constexpr double AMPLMinusInf()
-  { return -std::numeric_limits<double>::infinity(); }
+  { return -INFINITY; }
 
 
 public:

include/mp/common.h

@@ -23,6 +23,7 @@
 #ifndef MP_COMMON_H_
 #define MP_COMMON_H_
 
+#include <cmath>
 #include <cstddef>     // 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<double>::infinity() : 0;
+          -INFINITY : 0;
   }
 
   /** Constraint upper bound. */
   double con_ub() const {
     return (flags_ & INF_UB) != 0 ?
-          std::numeric_limits<double>::infinity() : 0;
+          INFINITY : 0;
   }
 };
 

include/mp/flat/constr_base.h

@@ -8,6 +8,9 @@
 #include <utility>
 #include <typeinfo>
 
+#include "mp/format.h"
+#include "mp/error.h"
+
 #include "mp/flat/context.h"
 #include "mp/arrayref.h"
 
@@ -104,6 +107,7 @@ using DblParamArray3 = ParamArrayN<double, 3>;
 /// Variable-length parameter array
 using DblParamArray = std::vector<double>;
 
+
 /// 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 <class VarVec>
+  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 <class Args, class Params,
+          class NumOrLogic, class Id, class VarVec>
+double ComputeValue(
+    const CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>& ,
+    const VarVec& ) {
+  MP_RAISE(fmt::format("ComputeValue({}) not implemented.",
+                       typeid(
+                         CustomFunctionalConstraint<Args,
+                         Params, NumOrLogic, Id>).name()));
+  return 0.0;
+}
+
+/// Compute violation of functional constraint.
+/// Can only be used for proper functional constraints,
+/// should be redefined for cones.
+template <class Args, class Params,
+          class NumOrLogic, class Id, class VarVec>
+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));
+}
+
+template <class Args, class Params,
+          class NumOrLogic, class Id>
+template <class VarVec>
+double
+CustomFunctionalConstraint<Args, Params, NumOrLogic, Id>
+::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 {

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 {
 

include/mp/flat/constr_viol.h

@@ -0,0 +1,72 @@
+#ifndef CONSTR_VIOL_H
+#define CONSTR_VIOL_H
+
+/**
+  * Violations of (mainly functional) constraints
+  */
+
+#include <cmath>
+
+#include "constr_functional.h"
+#include "constr_general.h"
+
+namespace mp {
+
+/// Compute result of the max constraint.
+template <class VarVec>
+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 <class VarVec>
+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 <class VarVec>
+double ComputeValue(const AbsConstraint& con, const VarVec& x) {
+  return std::fabs(x[con.GetArguments()[0]]);
+}
+
+/// Compute result of the and constraint.
+template <class VarVec>
+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 <class VarVec>
+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 <class VarVec>
+double ComputeValue(const NotConstraint& con, const VarVec& x) {
+  return x[con.GetArguments()[0]] < 0.5;
+}
+
+
+}  // namespace mp
+
+#endif // CONSTR_VIOL_H

include/mp/flat/preprocess.h

@@ -1,7 +1,6 @@
 #ifndef PREPROCESS_H
 #define PREPROCESS_H
 
-#include <limits>
 #include <cmath>
 #include <algorithm>
 

include/mp/flat/redef/conic/cones.h

@@ -8,7 +8,6 @@
 #define MP_FLAT_CVT_CONES_H
 
 #include <cmath>
-#include <limits>
 #include <vector>
 #include <cassert>
 

include/mp/nl-reader.h

@@ -48,7 +48,7 @@
 #include <algorithm>
 #include <cctype>
 #include <cstdlib>
-#include <limits>
+#include <cmath>
 #include <string>
 
 namespace mp {
@@ -1823,7 +1823,7 @@ void NLReader<Reader, Handler>::ReadBounds() {
   double lb = 0, ub = 0;
   BoundHandler bh(*this);
   int num_bounds = bh.num_items();
-  double infinity = std::numeric_limits<double>::infinity();
+  double infinity = INFINITY;
   for (int i = 0; i < num_bounds; ++i) {
     switch (reader_.ReadChar() - '0') {
     case RANGE:

solvers/mosek/mosekmodelapi.cc

@@ -1,3 +1,5 @@
+#include <cmath>
+
 #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<double>::infinity();
-    bool freeub = v.pub()[i] == std::numeric_limits<double>::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<double>::infinity())
+  if (lb == -INFINITY)
   {
-    if (ub == std::numeric_limits<double>::infinity())
+    if (ub == INFINITY)
       key = MSK_BK_FR;
     else
       key = MSK_BK_UP;
   }
   else
   {
-    if (ub == std::numeric_limits<double>::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<double>::infinity())
+  if (lb == -INFINITY)
   {
-    if (ub == std::numeric_limits<double>::infinity())
+    if (ub == INFINITY)
       key = MSK_BK_FR;
     else
       key = MSK_BK_UP;
   }
   else
   {
-    if (ub == std::numeric_limits<double>::infinity())
+    if (ub == INFINITY)
       key = MSK_BK_LO;
     else if (lb == ub)
       key = MSK_BK_FX;

src/expr-writer.h

@@ -23,6 +23,8 @@
 #ifndef MP_EXPR_WRITER_H_
 #define MP_EXPR_WRITER_H_
 
+#include <cmath>
+
 #include "mp/basic-expr-visitor.h"
 
 namespace mp {
@@ -396,7 +398,7 @@ void WriteExpr(fmt::Writer &w, const LinearExpr &linear,
 template <typename Problem>
 void Write(fmt::Writer &w, const Problem &p) {
   // Write variables.
-  double inf = std::numeric_limits<double>::infinity();
+  double inf = INFINITY;
   int num_vars = p.num_vars();
   for (int i = 0; i < num_vars; ++i) {
     w << "var x" << (i + 1);

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<double>::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);

test/function-test.cc

@@ -21,7 +21,6 @@
  */
 
 #include <functional>
-#include <limits>
 #include <sstream>
 #include <cmath>
 #include <cstring>
@@ -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<double>::infinity(), std::log(0.0));
+  EXPECT_EQ(-INFINITY, std::log(0.0));
   EXPECT_NE(0, isnan(diff(GetDoubleFun(std::log), 0)));
 }
 

test/nl-reader-test.cc

@@ -25,6 +25,7 @@
 // Include the source file to test the implementation.
 #include "../src/nl-reader.cc"
 
+#include <cmath>
 #include <climits>
 #include <cstring>
 
@@ -488,7 +489,7 @@ class TestNLHandler {
 
   void WriteBounds(char type, int index, double lb, double ub) {
     WriteSep();
-    double infinity = std::numeric_limits<double>::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<double>::infinity();
+  double inf = INFINITY;
   EXPECT_CALL(builder, SetComplementarity(66, 77, MatchComplInfo(0, inf)));
   adapter.OnComplementarity(66, 77, mp::ComplInfo(1));
 }

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<double>::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<double>::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<double>::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());