Improved numerical stability for sign constraints

src/configuration/GlobalConfiguration.cpp

@@ -92,6 +92,7 @@ const bool GlobalConfiguration::ONLY_AUX_INITIAL_BASIS = false;
 const GlobalConfiguration::ExplicitBasisBoundTighteningType GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE =
     GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX;
 const bool GlobalConfiguration::EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION = false;
+const double GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE = 0.000000001;
 
 const unsigned GlobalConfiguration::REFACTORIZATION_THRESHOLD = 100;
 const GlobalConfiguration::BasisFactorizationType GlobalConfiguration::BASIS_FACTORIZATION_TYPE =

src/configuration/GlobalConfiguration.h

@@ -189,6 +189,12 @@ class GlobalConfiguration
     // When doing explicit bound tightening, should we repeat until saturation?
     static const bool EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION;
 
+    /*
+      Tolerance for numerical imprecision of bounds computed as part
+      of explicit basis tightening
+    */
+    static const double EXPLICIT_BASIS_BOUND_TOLERANCE;
+
     /*
       Symbolic bound tightening options
     */

src/engine/RowBoundTightener.cpp

@@ -296,6 +296,36 @@ unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const TableauRow &r
         }
     }
 
+    // Attempt to reduce numerical imprecision
+    double storedLowerBound = _boundManager.getLowerBound( y );
+    double storedUpperBound = _boundManager.getUpperBound( y );
+
+    // If the lower bound is almost equal to the stored upper bound,
+    // make it exactly equal - the variable becomes fixed to a
+    // constant value
+    if ( FloatUtils::areEqual( lowerBound,
+                               storedUpperBound,
+                               GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+        lowerBound = storedUpperBound;
+
+    // If the upperbound is almost equal to the stored lower bound,
+    // make it exactly equal - the variable becomes fixed to a
+    // constant value
+    if ( FloatUtils::areEqual( upperBound,
+                               storedLowerBound,
+                               GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+        upperBound = storedLowerBound;
+
+    // If the compued lower and upper bounds are nearly equal, make
+    // them exactly equal, to their average
+    if ( FloatUtils::areEqual( upperBound,
+                               lowerBound,
+                               GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+    {
+        lowerBound = ( lowerBound + upperBound ) / 2;
+        upperBound = lowerBound;
+    }
+
     result += registerTighterLowerBound( y, lowerBound, row );
     result += registerTighterUpperBound( y, upperBound, row );
 
@@ -372,8 +402,39 @@ unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const TableauRow &r
             lowerBound = temp;
         }
 
-        // If a tighter bound is found, store it
         xi = row._row[i]._var;
+
+        // Attempt to reduce numerical imprecision
+        storedLowerBound = _boundManager.getLowerBound( xi );
+        storedUpperBound = _boundManager.getUpperBound( xi );
+
+        // If the lower bound is almost equal to the stored upper bound,
+        // make it exactly equal - the variable becomes fixed to a
+        // constant value
+        if ( FloatUtils::areEqual( lowerBound,
+                                   storedUpperBound,
+                                   GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+            lowerBound = storedUpperBound;
+
+        // If the upperbound is almost equal to the stored lower bound,
+        // make it exactly equal - the variable becomes fixed to a
+        // constant value
+        if ( FloatUtils::areEqual( upperBound,
+                                   storedLowerBound,
+                                   GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+            upperBound = storedLowerBound;
+
+        // If the compued lower and upper bounds are nearly equal, make
+        // them exactly equal, to their average
+        if ( FloatUtils::areEqual( upperBound,
+                                   lowerBound,
+                                   GlobalConfiguration::EXPLICIT_BASIS_BOUND_TOLERANCE ) )
+        {
+            lowerBound = ( lowerBound + upperBound ) / 2.0;
+            upperBound = lowerBound;
+        }
+
+        // If a tighter bound is found, store it
         result += registerTighterLowerBound( xi, lowerBound, row );
         result += registerTighterUpperBound( xi, upperBound, row );
         if ( FloatUtils::gt( getLowerBound( xi ), getUpperBound( xi ) ) )

src/engine/SignConstraint.cpp

@@ -114,10 +114,10 @@ bool SignConstraint::satisfied() const
 
     // if bValue is negative, f should be -1
     if ( FloatUtils::isNegative( bValue ) )
-        return FloatUtils::areEqual( fValue, -1 );
+        return FloatUtils::areEqual( fValue, -1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE );
 
     // bValue is non-negative, f should be 1
-    return FloatUtils::areEqual( fValue, 1 );
+    return FloatUtils::areEqual( fValue, 1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE );
 }
 
 List<PiecewiseLinearCaseSplit> SignConstraint::getCaseSplits() const
@@ -390,7 +390,7 @@ void SignConstraint::notifyLowerBound( unsigned variable, double bound )
     // Otherwise - update bound
     setLowerBound( variable, bound );
 
-    if ( variable == _f && FloatUtils::gt( bound, -1 ) )
+    if ( variable == _f && FloatUtils::gt( bound, -1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE ) )
     {
         setPhaseStatus( PhaseStatus::SIGN_PHASE_POSITIVE );
         if ( _boundManager != nullptr )
@@ -420,7 +420,7 @@ void SignConstraint::notifyUpperBound( unsigned variable, double bound )
     // Otherwise - update bound
     setUpperBound( variable, bound );
 
-    if ( variable == _f && FloatUtils::lt( bound, 1 ) )
+    if ( variable == _f && FloatUtils::lt( bound, 1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE ) )
     {
         setPhaseStatus( PhaseStatus::SIGN_PHASE_NEGATIVE );
         if ( _boundManager != nullptr )
@@ -478,14 +478,14 @@ void SignConstraint::getEntailedTightenings( List<Tightening> &tightenings ) con
 
     // Additional bounds can only be propagated if we are in the POSITIVE or NEGATIVE phases
     if ( !FloatUtils::isNegative( bLowerBound ) ||
-         FloatUtils::gt( fLowerBound, -1 ) )
+         FloatUtils::gt( fLowerBound, -1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE ) )
     {
         // Positive case
         tightenings.append( Tightening( _b, 0, Tightening::LB ) );
         tightenings.append( Tightening( _f, 1, Tightening::LB ) );
     }
     else if ( FloatUtils::isNegative( bUpperBound ) ||
-              FloatUtils::lt( fUpperBound, 1 ) )
+              FloatUtils::lt( fUpperBound, 1, GlobalConfiguration::CONSTRAINT_COMPARISON_TOLERANCE ) )
     {
         // Negative case
         tightenings.append( Tightening( _b, 0, Tightening::UB ) );