cvt:cmp:eps default 0 #102

Normalize conditional comparisons for better unification, makes positive tolerance unnecessary in current tests

Author: glebbelov

CHANGES.mp.md

@@ -1,10 +1,17 @@
 Summary of recent updates to the AMPL MP Library
 ================================================
 
+## unreleased
+- Changed default tolerance for strict comparisons
+  to 0 (option cvt:cmp:eps, #102.)
+- Fixed a bug where equivalent conditional
+  comparisons were not unified.
+
+
 ## 20230728
-- Option 'tech:writesolution' #218
-- Option 'writeprob' ('tech:writemodel') ASL-compatible
-- Hint when 'writeprob' fails: use 'writesol'
+- Option 'tech:writesolution' #218.
+- Option 'writeprob' ('tech:writemodel') ASL-compatible.
+- Hint when 'writeprob' fails: use 'writesol'.
 
 
 ## 20230726

include/mp/flat/constr_algebraic.h

@@ -70,7 +70,13 @@ class AlgebraicConstraint :
   /// Sorting and merging terms, some solvers require
   void sort_terms() { Body::sort_terms(); }
 
-  /// Negate
+  /// Is Normalized?
+  bool is_normalized() {
+    sort_terms();
+    return GetBody().is_normalized();
+  }
+
+  /// Negate all terms
   void negate() { Body::negate(); RhsOrRange::negate(); }
 
   /// Testing API

include/mp/flat/constr_base.h

@@ -35,7 +35,6 @@ class BasicConstraint {
   /// For functional constraints, result variable index
   int GetResultVar() const { return -1; }
 
-
 private:
   std::string name_;
 };
@@ -210,7 +209,8 @@ class ConditionalConstraint :
 
   /// Base class
   using Base = CustomFunctionalConstraint<
-    Con, ParamArray0, LogicalFunctionalConstraintTraits, CondConId<Con> >;
+    Con, ParamArray0,
+  LogicalFunctionalConstraintTraits, CondConId<Con> >;
 
   /// Default constructor
   ConditionalConstraint() = default;

include/mp/flat/constr_prepro.h

@@ -2,7 +2,13 @@
 #define CONSTR_PREPRO_H
 
 /**
- * Preprocess flat constraints before adding
+ * Preprocess flat constraints before adding.
+ *
+ * Possible tasks:
+ * 1. Simplify constraints
+ * 2. Replace a functional constraint by a different one,
+ *    via returning its result variable from another
+ *    (see conditional inequalities).
  */
 
 #include <cmath>
@@ -126,6 +132,8 @@ class ConstraintPreprocessors {
       CondLinConEQ& c, PreprocessInfo& prepro) {
     prepro.narrow_result_bounds(0.0, 1.0);
     prepro.set_result_type( var::INTEGER );
+    if (!IsNormalized(c))
+      c.GetConstraint().negate();   // for equality
     if (0!=MPD( IfPreproEqResBounds() ))
       if (FixEqualityResult(c, prepro))
         return;
@@ -135,6 +143,16 @@ class ConstraintPreprocessors {
         return;
   }
 
+  /// See if the argument of a conditional
+  /// algebraic constraint is normalized
+  template <class Body, int kind>
+  bool IsNormalized(
+      ConditionalConstraint<
+        AlgebraicConstraint< Body, AlgConRhs<kind> > >& cc) {
+    auto& arg = cc.GetConstraint();
+    return arg.is_normalized();
+  }
+
   /// Preprocess CondQuadConEQ
   template <class PreprocessInfo>
   void PreprocessConstraint(
@@ -218,9 +236,21 @@ class ConstraintPreprocessors {
       PreprocessInfo& prepro) {
     prepro.narrow_result_bounds(0.0, 1.0);
     prepro.set_result_type( var::INTEGER );
-    // See if we need to round the constant term
     assert(kind);
     auto& algc = cc.GetArguments();
+    if (!IsNormalized(cc)) {
+      auto arg1 = algc;
+      arg1.negate();    // Negate the terms and sense
+      prepro.set_result_var(
+            MPD( AssignResultVar2Args(
+                   ConditionalConstraint<
+                     AlgebraicConstraint< Body, AlgConRhs<
+                   -kind> > > { {
+                   std::move(arg1.GetBody()), arg1.rhs()
+                 } } ) ));
+      return;
+    }
+    // See if we need to round the constant term
     auto rhs = algc.rhs();
     auto bnt_body = MPD(
           ComputeBoundsAndType(algc.GetBody()) );

include/mp/flat/convert_functional.h

@@ -33,7 +33,7 @@ class BasicFCC {
       SetResultVar(GetConverter().
                    template GetConstraint<Constraint>(i).
                    GetResultVar());
-			GetConverter().IncrementVarUsage(GetResultVar());
+      GetConverter().IncrementVarUsage(GetResultVar()); // already here
 			if (GetConverter().DoingAutoLinking()) {  // Autolink known targets
         auto& varvn = GetConverter().GetVarValueNode();
         GetConverter().AutoLink( varvn.Select(GetResultVar()) );

include/mp/flat/expr_affine.h

@@ -110,6 +110,12 @@ class LinTerms {
       add_term(term.first, term.second);
   }
 
+  /// Is normalized? Assume terms are sorted.
+  bool is_normalized() const {
+    assert(size());
+    return coef(0) > 0.0;
+  }
+
   /// Negate
   void negate() {
     for (auto& c: coefs_)

include/mp/flat/expr_quadratic.h

@@ -73,6 +73,12 @@ class QuadTerms {
     vars2_.reserve(num_terms);
   }
 
+  /// Is normalized? Assume sorted.
+  bool is_normalized() const {
+    assert(size());
+    return coef(0) > 0.0;
+  }
+
   /// Arithmetic
   void negate() {
     for (auto& cf: coefs_)
@@ -172,6 +178,15 @@ class QuadAndLinTerms :
   /// add_term(c, v1, v2)
   using QuadTerms::add_term;
 
+  /// Is normalized? Assume sorted.
+  bool is_normalized() const {
+    assert(QuadTerms::size());
+    return
+        LinTerms::size()
+        ? LinTerms::is_normalized()
+        : QuadTerms::is_normalized();
+  }
+
   /// Negate
   void negate() {
     LinTerms::negate();

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

@@ -303,7 +303,7 @@ class MIPFlatConverter
 
 private:
   struct Options {
-    double cmpEps_ { 1e-4 };
+    double cmpEps_ { 0.0 };
     double bigM_default_ { -1 };
     double PLApproxRelTol_ { 1e-2 };
     double PLApproxDomain_ { 1e6 };
@@ -313,7 +313,9 @@ class MIPFlatConverter
   void InitOwnOptions() {
     this->GetEnv().AddOption("cvt:mip:eps cvt:cmp:eps",
                        "Tolerance for strict comparison of continuous variables for MIP. "
-                       "Ensure larger than the solver's feasibility tolerance.",
+                       "Also applies to negation of conditional comparisons: "
+                       "b==1 <==> x<=5 means that with b==0, x>=5+eps. "
+                             "Default: 0.",
                        options_.cmpEps_, 0.0, 1e100);
     this->GetEnv().AddOption("cvt:bigM cvt:bigm cvt:mip:bigM cvt:mip:bigm",
                        "Default value of big-M for linearization of logical constraints. "

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

@@ -4,7 +4,7 @@
     "files" : ["econ2.mod", "econ2.dat"],
     "tags" : ["linear", "continuous", "complementarity"],
     "options": {
-      "gcg_options": "cvt:bigm=1e5"
+      "ANYSOLVER_options": "cvt:bigm=1e5"
     },
     "values": {
       "Price['AA1']": 0.0,

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

@@ -11,7 +11,7 @@
     "tags": [ "logical" ],
     "objective": 22,
     "options": {
-      "cbc_options": "cvt:mip:bigM=100000",
+      "ANYSOLVER_options": "cvt:mip:bigM=100000",
       "gcg_options": "cvt:mip:bigM=100000 mode=2"
     }
   },
@@ -21,7 +21,7 @@
     "tags": [ "logical" ],
     "objective": 22,
     "options": {
-      "cbc_options": "cvt:mip:bigM=100000",
+      "ANYSOLVER_options": "cvt:mip:bigM=100000",
       "gcg_options": "cvt:mip:bigM=100000 mode=2"
     }
   },
@@ -80,7 +80,7 @@
     "tags" : ["logical"],
     "options": { 
       "solution_round": "6",
-      "gcg_options": "cvt:mip:bigM=100000" 
+      "ANYSOLVER_options": "cvt:mip:bigM=100000"
     },
     "comment_options": "For solution_round, see #200",
     "objective" : 218125

test/end2end/cases/categorized/fast/logical/ifthen_var.mod

@@ -1,4 +1,8 @@
-/** Test expression map as well as if-then */
+/**
+  * Test expression map as well as if-then.
+  * Requires cvt:cmp:eps > feastol
+  * for AMPL to correctly compute obj value, see #102.
+  */
 
 var x >=-100, <= 200;
 var y >=-300, <= 460;

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

@@ -15,7 +15,14 @@
   {
     "name" : "ifthen_var",
     "objective" : -5,
-    "tags" : ["logical"]
+    "tags" : ["logical"],
+    "options": {
+      "ANYSOLVER_options": "cvt:cmp:eps=1e-4"
+    },
+    "comment": [
+      "For AMPL to correctly compute obj value, ",
+      "need strict inequality for the opposite case"
+    ]
   },
   {
     "name" : "test_int_non_int",
@@ -113,17 +120,25 @@
   {
     "name" : "booleq_01",
     "options": {
-      "gcg_options": "cvt:bigm=1e5"
+      "ANYSOLVER_options": "cvt:bigm=1e5"
     },
     "objective" : 1,
     "tags" : ["logical"]
   },
   {
     "name" : "booleq_02",
     "options": {
-      "gcg_options": "cvt:bigm=1e5"
+      "ANYSOLVER_options": "cvt:bigm=1e5"
     },
     "objective" : 1,
     "tags" : ["logical"]
+  },
+  {
+    "name" : "x-multmip3_small",
+    "options": {
+      "ANYSOLVER_options": "cvt:bigm=1e5"
+    },
+    "objective" : 150,
+    "tags" : ["logical"]
   }
 ]

test/end2end/cases/categorized/fast/logical/x-multmip3_small.mod

@@ -0,0 +1,25 @@
+## Trying to replicate the cvt:cmp:eps problem from x-multmip3.mod.
+## Checks that all cases of x[i] >= minl are assigned
+## the same auxiliary variable, otherwise with cvt:cmp:eps=0
+## the optimal objective becomes 0. See also #102.
+
+param n default 4;
+param minl default 375;
+param limU default 500;
+param fcost default 50;
+param overall default 1200;
+param maxserve default n-1;
+
+var x{1..n} >= 0;
+
+minimize Total:
+   sum {i in 1..n} if x[i]>=minl then fcost;
+
+s.t. Lin01:
+   sum {i in 1..n} x[i] >= overall;
+
+s.t. Disj{i in 1..n}:
+   x[i]==0  or minl <= x[i] <= limU;
+
+s.t. Count:
+   count {i in 1..n} (x[i] >= minl) <= maxserve;

test/end2end/scripts/python/AMPLRunner.py

@@ -4,7 +4,7 @@
 from shutil import which
 
 from Solver import Solver
-from amplpy import AMPL, Kind, OutputHandler, ErrorHandler, Runnable, ampl
+from amplpy import AMPL, Kind, OutputHandler, ErrorHandler, ampl
 from Model import Model
 import time
 from TimeMe import TimeMe