Compute viol: PL, CondCon, Cones #200

include/mp/flat/constr_algebraic.h

@@ -70,11 +70,12 @@ class AlgebraicConstraint :
     return Body::ComputeValue(x) - RhsOrRange::lb();
   }
 
-  /// Compute violation
+  /// Compute violation.
+  /// Negative if if holds with slack.
   double
   ComputeViolation(const ArrayRef<double>& x) const {
     auto bd = Body::ComputeValue(x);
-    return std::max(
+    return std::max(      // same for strict cmp?
           RhsOrRange::lb() - bd, bd - RhsOrRange::ub());
   }
 

include/mp/flat/constr_base.h

@@ -190,14 +190,6 @@ double ComputeViolation(
   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
@@ -286,6 +278,27 @@ class ConditionalConstraint :
     return GetConstraint()==cc.GetConstraint();
   }
 
+  /// Compute violation of a conditional constraint.
+  /// If the subconstr is violated but should hold,
+  /// return the exact gap, despite this is a logical constraint.
+  /// If the subconstr holds but should not,
+  /// return the opposite gap. Similar to Indicator.
+  template <class VarVec>
+  double ComputeViolation(const VarVec& x) {
+    auto viol = GetConstraint().ComputeViolation(x);
+    bool ccon_valid = viol<=0.0;
+    bool has_arg = x[GetResultVar()] >= 0.5;
+    switch (this->GetContext().GetValue()) {
+    case Context::CTX_MIX:    // Viol is non-positive if holds
+      return has_arg == ccon_valid ? 0.0 : std::fabs(viol);
+    case Context::CTX_POS:
+      return has_arg <= ccon_valid ? 0.0 : viol;
+    case Context::CTX_NEG:
+      return has_arg >= ccon_valid ? 0.0 : -viol;
+    default:
+      return INFINITY;
+    }
+  }
 };
 
 ////////////////////////////////////////////////////////////////////////

include/mp/flat/constr_eval.h

@@ -235,29 +235,117 @@ double ComputeValue(const TanhConstraint& con, const VarVec& x) {
   return std::tanh(x[con.GetArguments()[0]]);
 }
 
+/// Compute result of the asinh constraint.
+template <class VarVec>
+double ComputeValue(const AsinhConstraint& con, const VarVec& x) {
+  return std::asinh(x[con.GetArguments()[0]]);
+}
+
+/// Compute result of the acosh constraint.
+template <class VarVec>
+double ComputeValue(const AcoshConstraint& con, const VarVec& x) {
+  return std::acosh(x[con.GetArguments()[0]]);
+}
+
+/// Compute result of the atanh constraint.
+template <class VarVec>
+double ComputeValue(const AtanhConstraint& con, const VarVec& x) {
+  return std::atanh(x[con.GetArguments()[0]]);
+}
+
 /// Compute result of a conditional constraint.
-/// Actually, for violation itself, we could do this:
-/// If the subconstr is violated but should hold,
-/// return the exact gap, despite this is a logical constraint.
-/// If the subconstr holds but should not, return 0.0.
+/// This is not used to compute violation.
 template <class Con, class VarVec>
 double ComputeValue(
     const ConditionalConstraint<Con>& con, const VarVec& x) {
   auto viol = con.GetConstraint().ComputeViolation(x);
   bool ccon_valid = viol<=0.0;
   bool has_arg = x[con.GetResultVar()] >= 0.5;
-  switch (con.GetContext()) {
+  switch (con.GetContext().GetValue()) {
   case Context::CTX_MIX:
     return has_arg == ccon_valid;
   case Context::CTX_POS:
-    return has_arg < ccon_valid;
+    return has_arg <= ccon_valid;
   case Context::CTX_NEG:
-    return has_arg > ccon_valid;
+    return has_arg >= ccon_valid;
   default:
     return INFINITY;
   }
 }
 
+/// Compute violation of the QuadraticCone constraint.
+template <class VarVec>
+double ComputeViolation(
+    const QuadraticConeConstraint& con, const VarVec& x) {
+  const auto& args = con.GetArguments();
+  const auto& params = con.GetParameters();
+  assert(args.size()==params.size());
+  double sum = 0.0;
+  for (auto i=args.size(); --i; )
+    sum += std::pow( params[i]*x[args[i]], 2.0 );
+  return std::sqrt(sum) - params[0]*x[args[0]];
+}
+
+/// Compute violation of the RotatedQuadraticCone constraint.
+template <class VarVec>
+double ComputeViolation(
+    const RotatedQuadraticConeConstraint& con, const VarVec& x) {
+  const auto& args = con.GetArguments();
+  const auto& params = con.GetParameters();
+  assert(args.size()==params.size());
+  double sum = 0.0;
+  for (auto i=args.size(); --i>1; )
+    sum += std::pow( params[i]*x[args[i]], 2.0 );
+  return sum
+      - 2.0 * params[0]*x[args[0]] * params[1]*x[args[1]];
+}
+
+/// Compute violation of the ExponentialCone constraint.
+template <class VarVec>        // ax >= by exp(cz / (by))
+double ComputeViolation(       // where ax, by >= 0
+    const ExponentialConeConstraint& con, const VarVec& x) {
+  const auto& args = con.GetArguments();
+  const auto& params = con.GetParameters();
+  auto ax = params[0]*x[args[0]];
+  auto by = params[1]*x[args[1]];
+  if (0.0==std::fabs(by))
+    return -ax;
+  auto cz = params[2]*x[args[2]];
+  return by * std::exp(cz / by) - ax;
+}
+
+/// Compute result of the PL constraint.
+template <class VarVec>
+double ComputeValue(const PLConstraint& con, const VarVec& x) {
+  const auto& plp = con.GetParameters().GetPLPoints();
+  assert(!plp.empty());
+  auto x0 = x[con.GetArguments()[0]];        // position
+  if (x0<plp.x_.front())
+    return plp.y_.front()
+        - plp.PreSlope()*(plp.x_.front() - x0);
+  if (x0>plp.x_.back())
+    return plp.y_.back()
+        + plp.PostSlope()*(x0 - plp.x_.back());
+  int i0=0;
+  for ( ; x0 > plp.x_[i0]; ++i0);
+  return plp.x_[i0]==x0
+      ? plp.y_[i0]
+        : plp.y_[i0]
+        + (plp.y_[i0+1]-plp.y_[i0])
+      * (x0-plp.x_[i0]) / (plp.x_[i0+1]-plp.x_[i0]);
+}
+
+/// Should be here,
+/// after ComputeViolation() is specialized
+/// for some constraints.
+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); }
+
 
 }  // namespace mp
 

include/mp/flat/constr_functional.h

@@ -359,11 +359,15 @@ class PLConParams {
   /// Construct from PLPoints
   PLConParams(PLPoints plp) : plp_(std::move(plp)) { }
   /// Produce PLPoints, either stored or from PLSlopes
-  operator PLPoints() const { return plp_.empty() ? pls_ : plp_; }
+  const PLPoints& GetPLPoints() const {
+    if (plp_.empty())
+      plp_ = pls_;
+    return plp_;
+  }
 
 private:
   PLSlopes pls_;
-  PLPoints plp_;
+  mutable PLPoints plp_;
 };
 
 

include/mp/flat/constr_general.h

@@ -223,7 +223,8 @@ DEF_STATIC_CONSTR_WITH_PRM( QuadraticConeConstraint, VarArray, DblParamArray,
                             " with factors p1..pn");
 /// Rotated quadratic cone
 DEF_STATIC_CONSTR_WITH_PRM( RotatedQuadraticConeConstraint, VarArray, DblParamArray,
-                            "Rotated quadratic cone p1*x1*p2*x2 >= sqrt((p3*x3)^2 + ...)),"
+                            "Rotated quadratic cone "
+                            "2 * p1*x1*p2*x2 >= (p3*x3)^2 + ...),"
                             " x1, x2 >= 0, with factors p1..pn");
 /// Exponential cone
 DEF_STATIC_CONSTR_WITH_PRM( ExponentialConeConstraint, VarArray3, DblParamArray3,

include/mp/flat/redef/MIP/piecewise_linear.h

@@ -29,7 +29,7 @@ class PLConverter_MIP :
 
   /// Convert in any context
   void Convert(const ItemType& cc, int ) {
-    points_ = cc.GetParameters();      // extract PLPoints
+    points_ = cc.GetParameters().GetPLPoints();
     i0=0;                              // first breakpoint
     i1=points_.x_.size()-1;            // last breakpoint
     y = cc.GetResultVar();

solvers/cplexmp/cplexmpmodelapi.cc

@@ -141,7 +141,7 @@ void CplexModelAPI::AddConstraint(const IndicatorConstraintLinGE &ic)  {
 }
 
 void CplexModelAPI::AddConstraint(const PLConstraint& plc) {
-  PLPoints plp(plc.GetParameters());
+  const auto& plp = plc.GetParameters().GetPLPoints();
   CPLEX_CALL( CPXaddpwl(env(), lp(),
               plc.GetResultVar(), plc.GetArguments()[0],
               plp.PreSlope(), plp.PostSlope(),

solvers/gurobi/gurobimodelapi.cc

@@ -225,7 +225,7 @@ void GurobiModelAPI::AddConstraint(const TanConstraint &cc)  {
 }
 
 void GurobiModelAPI::AddConstraint(const PLConstraint& plc) {
-  PLPoints plp(plc.GetParameters());
+  const auto& plp = plc.GetParameters().GetPLPoints();
   GRB_CALL( GRBaddgenconstrPWL(model(), plc.name(),
               plc.GetArguments()[0], plc.GetResultVar(),
               plp.x_.size(), plp.x_.data(), plp.y_.data()) );

test/end2end/cases/categorized/fast/nonlinear/expA_1.mod

@@ -4,14 +4,14 @@
 # expA_1.mod
 # -------------------------------------------------------------
 
-param ubx integer := 10;
-param uby integer := 20;
+param ubx integer := 100;
+param uby integer := 200;
 
-var x  >= 0, <= ubx;
-var y  >= -41, <= uby;
+var x  >= -100, <= ubx;
+var y  >= -205, <= uby;
 
 maximize TotalSum:
-    x-y;
+    -x-y;
 
 subj to ExpA: 
     y >= 0.1 ^ x;

test/end2end/cases/categorized/fast/nonlinear/modellist.json

@@ -3,15 +3,15 @@
     "name" : "expA_1",
     "tags" : ["nonlinear"],
     "values": {
-      "if abs(_obj[1]-10)<=0.01 then 1 else 0": 1
+      "if abs(_obj[1]+0.7977905982)<=0.01 then 1 else 0": 1
     }
   },
   {
     "name" : "expA_1 check_pl_approx_expA",
     "tags" : ["nonlinear", "check_pl_approx_expA"],
     "options": { "ANYSOLVER_options": "acc:expa=1" },
     "values": {
-      "if abs(_obj[1]-10)<=0.01 then 1 else 0": 1
+      "if abs(_obj[1]+0.7977905982)<=0.01 then 1 else 0": 1
     }
   },
   {