Simulator + numerical stability in preprocessor (#91)

* simulator

* fix some issues in the preprocessor

* bug fix

* a constant for almost-fixed vars in the preprocessor

* numerical stability in the preprocessor, some work on the simulator

src/configuration/GlobalConfiguration.cpp

@@ -42,6 +42,7 @@ const bool GlobalConfiguration::USE_HARRIS_RATIO_TEST = true;
 const bool GlobalConfiguration::PREPROCESS_INPUT_QUERY = true;
 const bool GlobalConfiguration::PREPROCESSOR_ELIMINATE_VARIABLES = true;
 const bool GlobalConfiguration::PREPROCESSOR_PL_CONSTRAINTS_ADD_AUX_EQUATIONS = false;
+const double GlobalConfiguration::PREPROCESSOR_ALMOST_FIXED_THRESHOLD = 0.00001;
 
 const unsigned GlobalConfiguration::PSE_ITERATIONS_BEFORE_RESET = 1000;
 const double GlobalConfiguration::PSE_GAMMA_ERROR_THRESHOLD = 0.001;

src/configuration/GlobalConfiguration.h

@@ -63,6 +63,10 @@ class GlobalConfiguration
     // to add auxiliary variables and equations.
     static const bool PREPROCESSOR_PL_CONSTRAINTS_ADD_AUX_EQUATIONS;
 
+    // If the difference between a variable's lower and upper bounds is smaller than this
+    // threshold, the preprocessor will treat it as fixed.
+    static const double PREPROCESSOR_ALMOST_FIXED_THRESHOLD;
+
     // How often should the main loop check the current degradation?
     static const unsigned DEGRADATION_CHECKING_FREQUENCY;
 

src/engine/Engine.cpp

@@ -677,12 +677,24 @@ void Engine::extractSolution( InputQuery &inputQuery )
     {
         if ( _preprocessingEnabled )
         {
+            // Fixed variables are easy: return the value they've been fixed to.
             if ( _preprocessor.variableIsFixed( i ) )
+            {
                 inputQuery.setSolutionValue( i, _preprocessor.getFixedValue( i ) );
-            else if ( _preprocessor.variableIsMerged( i ) )
-                inputQuery.setSolutionValue( i, _tableau->getValue( _preprocessor.getMergedIndex( i ) ) );
-            else
-                inputQuery.setSolutionValue( i, _tableau->getValue( _preprocessor.getNewIndex( i ) ) );
+                continue;
+            }
+
+            // Has the variable been merged into another?
+            unsigned variable;
+            if ( _preprocessor.variableIsMerged( i ) )
+                variable = _preprocessor.getMergedIndex( i );
+
+            // We know which variable to look for, but it may have been assigned
+            // a new index, due to variable elimination
+            unsigned finalIndex = _preprocessor.getNewIndex( variable );
+
+            // Finally, set the assigned value
+            inputQuery.setSolutionValue( i, _tableau->getValue( finalIndex ) );
         }
         else
         {

src/engine/InputQuery.cpp

@@ -263,6 +263,15 @@ unsigned InputQuery::getNumOutputVariables() const
     return _outputIndexToVariable.size();
 }
 
+List<unsigned> InputQuery::getInputVariables() const
+{
+    List<unsigned> result;
+    for ( const auto &pair : _variableToInputIndex )
+        result.append( pair.first );
+
+    return result;
+}
+
 void InputQuery::printInputOutputBounds() const
 {
     printf( "Dumping bounds of the input and output variables:\n" );

src/engine/InputQuery.h

@@ -55,6 +55,7 @@ class InputQuery
     unsigned outputVariableByIndex( unsigned index ) const;
     unsigned getNumInputVariables() const;
     unsigned getNumOutputVariables() const;
+    List<unsigned> getInputVariables() const;
 
     /*
       Methods for setting and getting the solution.

src/engine/Preprocessor.cpp

@@ -58,12 +58,16 @@ InputQuery Preprocessor::preprocess( const InputQuery &query, bool attemptVariab
     {
         continueTightening = processEquations();
         continueTightening = processConstraints() || continueTightening;
-        continueTightening = processIdenticalVariables() || continueTightening;
+
+        if ( attemptVariableElimination )
+            continueTightening = processIdenticalVariables() || continueTightening;
 
         if ( _statistics )
             _statistics->ppIncNumTighteningIterations();
     }
 
+    collectFixedValues();
+
     if ( attemptVariableElimination )
         eliminateVariables();
 
@@ -111,13 +115,16 @@ bool Preprocessor::processEquations()
 
     bool tighterBoundFound = false;
 
-    for ( const auto &equation : _preprocessed.getEquations() )
+    List<Equation> &equations( _preprocessed.getEquations() );
+    List<Equation>::iterator equation = equations.begin();
+
+    while ( equation != equations.end() )
     {
         // The equation is of the form sum (ci * xi) - b ? 0
-        Equation::EquationType type = equation._type;
+        Equation::EquationType type = equation->_type;
 
-        unsigned maxVar = equation._addends.begin()->_variable;
-        for ( const auto &addend : equation._addends )
+        unsigned maxVar = equation->_addends.begin()->_variable;
+        for ( const auto &addend : equation->_addends )
         {
             if ( addend._variable > maxVar )
                 maxVar = addend._variable;
@@ -143,9 +150,9 @@ bool Preprocessor::processEquations()
 
         // The first goal is to compute the LB and UB of: sum (ci * xi) - b
         // For this we first identify unbounded variables
-        double auxLb = -equation._scalar;
-        double auxUb = -equation._scalar;
-        for ( const auto &addend : equation._addends )
+        double auxLb = -equation->_scalar;
+        double auxUb = -equation->_scalar;
+        for ( const auto &addend : equation->_addends )
         {
             ci = addend._coefficient;
             xi = addend._variable;
@@ -158,7 +165,7 @@ bool Preprocessor::processEquations()
                 continue;
             }
 
-            ciSign[xi] = FloatUtils::isPositive( ci ) ? POSITIVE : NEGATIVE;
+            ciSign[xi] = ci > 0 ? POSITIVE : NEGATIVE;
 
             xiLB = _preprocessed.getLowerBound( xi );
             xiUB = _preprocessed.getUpperBound( xi );
@@ -173,7 +180,7 @@ bool Preprocessor::processEquations()
             }
             else
             {
-                if ( FloatUtils::isPositive( ci ) )
+                if ( ci > 0 )
                     excludedFromLB.insert( xi );
                 else
                     excludedFromUB.insert( xi );
@@ -189,15 +196,15 @@ bool Preprocessor::processEquations()
             }
             else
             {
-                if ( FloatUtils::isPositive( ci ) )
+                if ( ci > 0 )
                     excludedFromUB.insert( xi );
                 else
                     excludedFromLB.insert( xi );
             }
         }
 
         // Now, go over each addend in sum (ci * xi) - b ? 0, and see what can be done
-        for ( const auto &addend : equation._addends )
+        for ( const auto &addend : equation->_addends )
         {
             ci = addend._coefficient;
             xi = addend._variable;
@@ -268,7 +275,7 @@ bool Preprocessor::processEquations()
 
                 lowerBound /= -ci;
 
-                if ( FloatUtils::gt( lowerBound, _preprocessed.getLowerBound( xi ) ) )
+                if ( lowerBound > _preprocessed.getLowerBound( xi ) )
                 {
                     tighterBoundFound = true;
                     _preprocessed.setLowerBound( xi, lowerBound );
@@ -292,14 +299,16 @@ bool Preprocessor::processEquations()
 
                 upperBound /= -ci;
 
-                if ( FloatUtils::lt( upperBound, _preprocessed.getUpperBound( xi ) ) )
+                if ( upperBound < _preprocessed.getUpperBound( xi ) )
                 {
                     tighterBoundFound = true;
                     _preprocessed.setUpperBound( xi, upperBound );
                 }
             }
 
-            if ( FloatUtils::gt( _preprocessed.getLowerBound( xi ), _preprocessed.getUpperBound( xi ) ) )
+            if ( FloatUtils::gt( _preprocessed.getLowerBound( xi ),
+                                 _preprocessed.getUpperBound( xi ),
+                                 GlobalConfiguration::PREPROCESSOR_ALMOST_FIXED_THRESHOLD ) )
             {
                 delete[] ciTimesLb;
                 delete[] ciTimesUb;
@@ -312,6 +321,41 @@ bool Preprocessor::processEquations()
         delete[] ciTimesLb;
         delete[] ciTimesUb;
         delete[] ciSign;
+
+        /*
+          Next, do another sweep over the equation.
+          Look for almost-fixed variables and fix them, and remove the equation
+          entirely if it has nothing left to contribute.
+        */
+        bool allFixed = true;
+        for ( const auto &addend : equation->_addends )
+        {
+            unsigned var = addend._variable;
+            double lb = _preprocessed.getLowerBound( var );
+            double ub = _preprocessed.getUpperBound( var );
+
+            if ( FloatUtils::areEqual( lb, ub, GlobalConfiguration::PREPROCESSOR_ALMOST_FIXED_THRESHOLD ) )
+                _preprocessed.setUpperBound( var, _preprocessed.getLowerBound( var ) );
+            else
+                allFixed = false;
+        }
+
+        if ( !allFixed )
+        {
+            ++equation;
+        }
+        else
+        {
+            double sum = 0;
+            for ( const auto &addend : equation->_addends )
+                sum += addend._coefficient * _preprocessed.getLowerBound( addend._variable );
+
+            if ( FloatUtils::areDisequal( sum, equation->_scalar, GlobalConfiguration::PREPROCESSOR_ALMOST_FIXED_THRESHOLD ) )
+            {
+                throw InfeasibleQueryException();
+            }
+            equation = equations.erase( equation );
+        }
     }
 
     return tighterBoundFound;
@@ -335,18 +379,24 @@ bool Preprocessor::processConstraints()
         for ( const auto &tightening : tightenings )
 		{
 			if ( ( tightening._type == Tightening::LB ) &&
-                 FloatUtils::gt( tightening._value, _preprocessed.getLowerBound( tightening._variable ) ) )
+                 ( tightening._value > _preprocessed.getLowerBound( tightening._variable ) ) )
             {
                 tighterBoundFound = true;
                 _preprocessed.setLowerBound( tightening._variable, tightening._value );
             }
 
 			else if ( ( tightening._type == Tightening::UB ) &&
-				 FloatUtils::lt( tightening._value, _preprocessed.getUpperBound( tightening._variable ) ) )
+                      ( tightening._value < _preprocessed.getUpperBound( tightening._variable ) ) )
             {
                 tighterBoundFound = true;
                 _preprocessed.setUpperBound( tightening._variable, tightening._value );
             }
+
+            if ( FloatUtils::areEqual( _preprocessed.getLowerBound( tightening._variable ),
+                                       _preprocessed.getUpperBound( tightening._variable ),
+                                       GlobalConfiguration::PREPROCESSOR_ALMOST_FIXED_THRESHOLD ) )
+                _preprocessed.setUpperBound( tightening._variable,
+                                             _preprocessed.getLowerBound( tightening._variable ) );
 		}
 	}
 
@@ -355,58 +405,70 @@ bool Preprocessor::processConstraints()
 
 bool Preprocessor::processIdenticalVariables()
 {
-    // Find distinct v1 and v2 which are exactly equal to each other
-    unsigned v1 = 0, v2 = 0;
-    Equation equToRemove;
-    for ( const auto &equation : _preprocessed.getEquations() )
+    List<Equation> &equations( _preprocessed.getEquations() );
+    List<Equation>::iterator equation = equations.begin();
+
+    bool found = false;
+    while ( equation != equations.end() )
     {
-        if ( equation._addends.size() != 2 || equation._type != Equation::EQ )
+        // We are only looking for equations of type c(v1 - v2) = 0
+        if ( equation->_addends.size() != 2 || equation->_type != Equation::EQ )
+        {
+            ++equation;
             continue;
+        }
 
-        auto term1 = equation._addends.front();
-        auto term2 = equation._addends.back();
+        Equation::Addend term1 = equation->_addends.front();
+        Equation::Addend term2 = equation->_addends.back();
 
         if ( FloatUtils::areDisequal( term1._coefficient, -term2._coefficient ) ||
-             !FloatUtils::isZero( equation._scalar ) )
+             !FloatUtils::isZero( equation->_scalar ) )
+        {
+            ++equation;
             continue;
+        }
 
         ASSERT( term1._variable != term2._variable );
-        if ( term1._variable == term2._variable )
-            continue;
 
-        v1 = term1._variable;
-        v2 = term2._variable;
-        equToRemove = equation;
-        break;
-    }
+        // The equation matches the pattern, process and remove it
+        found = true;
+        unsigned v1 = term1._variable;
+        unsigned v2 = term2._variable;
+
+        double bestLowerBound =
+            _preprocessed.getLowerBound( v1 ) > _preprocessed.getLowerBound( v2 ) ?
+            _preprocessed.getLowerBound( v1 ) :
+            _preprocessed.getLowerBound( v2 );
+
+        double bestUpperBound =
+            _preprocessed.getUpperBound( v1 ) < _preprocessed.getUpperBound( v2 ) ?
+            _preprocessed.getUpperBound( v1 ) :
+            _preprocessed.getUpperBound( v2 );
+
+        equation = equations.erase( equation );
 
-    if ( v1 == v2 )
-        return false;
+        _preprocessed.setLowerBound( v2, bestLowerBound );
+        _preprocessed.setUpperBound( v2, bestUpperBound );
 
-    // Found v1 and v2 which are identical
-    _preprocessed.removeEquation( equToRemove );
-    _preprocessed.mergeIdenticalVariables( v1, v2 );
+        _preprocessed.mergeIdenticalVariables( v1, v2 );
 
-    double bestLowerBound = FloatUtils::max( _preprocessed.getLowerBound( v1 ),
-                                             _preprocessed.getLowerBound( v2 ) );
-    double bestUpperBound = FloatUtils::min( _preprocessed.getUpperBound( v1 ),
-                                             _preprocessed.getUpperBound( v2 ) );
-    _preprocessed.setLowerBound( v2, bestLowerBound );
-    _preprocessed.setUpperBound( v2, bestUpperBound );
+        _mergedVariables[v1] = v2;
+    }
 
-    _mergedVariables[v1] = v2;
-    return true;
+    return found;
 }
 
-void Preprocessor::eliminateVariables()
+void Preprocessor::collectFixedValues()
 {
-    // First, collect the variables that have become fixed, and their fixed values
 	for ( unsigned i = 0; i < _preprocessed.getNumberOfVariables(); ++i )
 	{
         if ( FloatUtils::areEqual( _preprocessed.getLowerBound( i ), _preprocessed.getUpperBound( i ) ) )
             _fixedVariables[i] = _preprocessed.getLowerBound( i );
 	}
+}
 
+void Preprocessor::eliminateVariables()
+{
     // If there's nothing to eliminate, we're done
     if ( _fixedVariables.empty() && _mergedVariables.empty() )
         return;

src/engine/Preprocessor.h

@@ -73,6 +73,11 @@ class Preprocessor
     */
     bool processIdenticalVariables();
 
+    /*
+      Collect all variables whose lower and upper bounds are equal
+    */
+    void collectFixedValues();
+
     /*
       Eliminate any variables that have become fixed or merged with an identical variable
 	*/

src/engine/ReluplexError.h

@@ -36,6 +36,7 @@ class ReluplexError : public Error
         RESTORATION_FAILED_TO_REFACTORIZE_BASIS = 14,
         ENGINE_APPLY_SPLIT_FAILED = 15,
         FILE_DOESNT_EXIST = 16,
+        SIMULATOR_ERROR = 17,
 
         DEBUGGING_ERROR = 999,
     };