Numerical stability issue in SFT (#94)

* WIP

* WIP

* wip

* cleanup

* cleanup

* cleanup

* cleanup

* compiler warning

src/basis_factorization/CSRMatrix.cpp

@@ -574,7 +574,7 @@ unsigned CSRMatrix::findArrayIndexForEntry( unsigned row, unsigned column ) cons
 {
     int low = _IA[row];
     int high = _IA[row + 1] - 1;
-    int mid;
+    int mid = -1;
 
     bool found = false;
     while ( !found && low <= high )

src/basis_factorization/IBasisFactorization.h

@@ -16,6 +16,7 @@
 class SparseColumnsOfBasis;
 class SparseMatrix;
 class SparseUnsortedList;
+class Statistics;
 
 class IBasisFactorization
 {
@@ -102,6 +103,11 @@ class IBasisFactorization
      */
     virtual void invertBasis( double *result ) = 0;
 
+    /*
+      Have the Basis Factoriaztion object start reporting statistics.
+    */
+    virtual void setStatistics( Statistics * ) {};
+
     /*
       For debugging
     */

src/basis_factorization/SparseFTFactorization.cpp

@@ -23,6 +23,7 @@ SparseFTFactorization::SparseFTFactorization( unsigned m, const BasisColumnOracl
 	, _m( m )
     , _sparseLUFactors( m )
     , _sparseGaussianEliminator( m )
+    , _statistics( NULL )
     , _z1( NULL )
     , _z2( NULL )
     , _z3( NULL )
@@ -161,7 +162,10 @@ void SparseFTFactorization::updateToAdjacentBasis( unsigned columnIndex,
 
     ASSERT( foundNonZeroEntry );
     if ( lastNonZeroEntryInU <= uColumnIndex )
+    {
+        ASSERT( uColumnIndex == lastNonZeroEntryInU ); // Otherwise, singular matrix
         return;
+    }
 
     /*
       Step 3:
@@ -261,6 +265,16 @@ void SparseFTFactorization::updateToAdjacentBasis( unsigned columnIndex,
         }
     }
 
+    if ( -GlobalConfiguration::SPARSE_FORREST_TOMLIN_DIAGONAL_ELEMENT_TOLERANCE <
+         _z3[columnIndex]
+         &&
+         _z3[columnIndex] <
+         GlobalConfiguration::SPARSE_FORREST_TOMLIN_DIAGONAL_ELEMENT_TOLERANCE )
+    {
+        obtainFreshBasis();
+        return;
+    }
+
     /*
       Step 5:
 
@@ -538,6 +552,12 @@ void SparseFTFactorization::dumpExplicitBasis() const
     }
 }
 
+void SparseFTFactorization::setStatistics( Statistics *statistics )
+{
+    _statistics = statistics;
+    _sparseGaussianEliminator.setStatistics( statistics );
+}
+
 //
 // Local Variables:
 // compile-command: "make -C ../.. "

src/basis_factorization/SparseFTFactorization.h

@@ -18,6 +18,7 @@
 #include "SparseEtaMatrix.h"
 #include "SparseGaussianEliminator.h"
 #include "SparseLUFactors.h"
+#include "Statistics.h"
 
 /*
   This class performs a sparse FT factorization of a given matrix.
@@ -125,6 +126,11 @@ class SparseFTFactorization : public IBasisFactorization
     */
     SparseGaussianEliminator _sparseGaussianEliminator;
 
+    /*
+      An object for reporting statistics
+    */
+    Statistics *_statistics;
+
     /*
       Work memory.
     */
@@ -165,6 +171,11 @@ class SparseFTFactorization : public IBasisFactorization
     */
 	void clearFactorization();
 
+    /*
+      Have the Basis Factoriaztion object start reporting statistics.
+    */
+    void setStatistics( Statistics *statistics );
+
     static void log( const String &message );
 };
 

src/basis_factorization/SparseGaussianEliminator.cpp

@@ -25,9 +25,9 @@ SparseGaussianEliminator::SparseGaussianEliminator( unsigned m )
     : _m( m )
     , _work( NULL )
     , _work2( NULL )
+    , _statistics( NULL )
     , _numURowElements( NULL )
     , _numUColumnElements( NULL )
-
 {
     _work = new double[_m];
     if ( !_work )
@@ -137,6 +137,34 @@ void SparseGaussianEliminator::run( const SparseColumnsOfBasis *A, SparseLUFacto
 
     // Do the work
     factorize();
+
+    // DEBUG({
+    //         // Check that the factorization is correct
+    //         double *product = new double[_m * _m];
+    //         std::fill_n( product, _m * _m, 0 );
+
+    //         for ( unsigned i = 0; i < _m; ++i )
+    //             for ( unsigned j = 0; j < _m; ++j )
+    //                 for ( unsigned k = 0; k < _m; ++k )
+    //                 {
+    //                     double fValue = ( i == k ) ? 1.0 : _sparseLUFactors->_F->get( i, k );
+    //                     double vValue = _sparseLUFactors->_V->get( k, j );
+
+    //                     ASSERT( FloatUtils::wellFormed( fValue ) );
+    //                     ASSERT( FloatUtils::wellFormed( vValue ) );
+
+    //                     product[i*_m + j] += fValue * vValue;
+    //                 }
+
+    //         for ( unsigned i = 0; i < _m; ++i )
+    //             for ( unsigned j = 0; j < _m; ++j )
+    //             {
+    //                 ASSERT( FloatUtils::areEqual( product[i*_m+j],
+    //                                               A->_columns[j]->get( i ) ) );
+    //             }
+
+    //         delete[] product;
+    //     });
 }
 
 void SparseGaussianEliminator::factorize()
@@ -147,7 +175,7 @@ void SparseGaussianEliminator::factorize()
         /*
           Step 1:
           -------
-           Choose a pivot element from the active submatrix of U. This
+          Choose a pivot element from the active submatrix of U. This
           can be any non-zero coefficient. Store the result in:
              _uPivotRow, _uPivotColumn (indices in U)
              _vPivotRow, _vPivotColumn (indices in V)
@@ -171,29 +199,6 @@ void SparseGaussianEliminator::factorize()
         */
         eliminate();
     }
-
-    // DEBUG({
-    //         // Check that the factorization is correct
-    //         double *product = new double[_m * _m];
-    //         std::fill_n( product, _m * _m, 0 );
-
-    //         for ( unsigned i = 0; i < _m; ++i )
-    //             for ( unsigned j = 0; j < _m; ++j )
-    //                 for ( unsigned k = 0; k < _m; ++k )
-    //                 {
-    //                     double fValue = ( i == k ) ? 1.0 : _sparseLUFactors->_F->get( i, k );
-    //                     double vValue = _sparseLUFactors->_V->get( k, j );
-    //                     product[i*_m + j] += fValue * vValue;
-    //                 }
-
-    //         for ( unsigned i = 0; i < _m; ++i )
-    //             for ( unsigned j = 0; j < _m; ++j )
-    //             {
-    //                 ASSERT( FloatUtils::areEqual( product[i*_m+j], A->get( i, j ) ) );
-    //             }
-
-    //         delete[] product;
-    //     });
 }
 
 void SparseGaussianEliminator::choosePivot()
@@ -283,7 +288,7 @@ void SparseGaussianEliminator::choosePivot()
     }
 
     // No singletons, apply the Markowitz rule. Find the element with acceptable
-    // magnitude that has the smallet Markowitz rule.
+    // magnitude that has the smallet Markowitz value.
     // Fail if no elements exists that are within acceptable magnitude
 
     // Todo: more clever heuristics to reduce the search space
@@ -406,7 +411,6 @@ void SparseGaussianEliminator::eliminate()
           The multiplier is: - U[row,k] / pivotElement
         */
         double rowMultiplier = - columnIt->_value / _pivotElement;
-        log( Stringf( "\tWorking on V row: %u. Multiplier: %lf", vRow, rowMultiplier ) );
 
         // Get the row being eliminated in dense format
         _sparseLUFactors->_V->getRowDense( vRow, _work2 );
@@ -473,6 +477,11 @@ void SparseGaussianEliminator::log( const String &message )
         printf( "SparseGaussianEliminator: %s\n", message.ascii() );
 }
 
+void SparseGaussianEliminator::setStatistics( Statistics *statistics )
+{
+    _statistics = statistics;
+}
+
 //
 // Local Variables:
 // compile-command: "make -C ../.. "

src/basis_factorization/SparseGaussianEliminator.h

@@ -17,6 +17,7 @@
 #include "SparseColumnsOfBasis.h"
 #include "SparseLUFactors.h"
 #include "SparseMatrix.h"
+#include "Statistics.h"
 
 class SparseGaussianEliminator
 {
@@ -30,6 +31,11 @@ class SparseGaussianEliminator
     */
     void run( const SparseColumnsOfBasis *A, SparseLUFactors *sparseLUFactors );
 
+    /*
+      Have the eliminator start reporting statistics.
+    */
+    void setStatistics( Statistics *statistics );
+
 private:
     /*
       The dimension of the (square) matrix being factorized
@@ -57,6 +63,11 @@ class SparseGaussianEliminator
     double *_work;
     double *_work2;
 
+    /*
+      An object for reporting statistics
+    */
+    Statistics *_statistics;
+
     /*
       Information on the number of non-zero elements in
       every row of the current active submatrix of U

src/basis_factorization/SparseLUFactors.cpp

@@ -289,6 +289,8 @@ void SparseLUFactors::vForwardTransformation( const double *y, double *x ) const
                 x[xBeingSolved] -= ( entry._value * x[vColumn] );
         }
 
+        ASSERT( !FloatUtils::isZero( diagonalCoefficient ) );
+
         if ( FloatUtils::isZero( x[xBeingSolved] ) )
             x[xBeingSolved] = 0.0;
         else

src/common/Sources.mk

@@ -6,6 +6,7 @@ SOURCES += \
 	HeapData.cpp \
 	MString.cpp \
 	SignalHandler.cpp \
+	Statistics.cpp \
 	TimeUtils.cpp \
 
 #

src/configuration/GlobalConfiguration.cpp

@@ -27,6 +27,7 @@ const double GlobalConfiguration::HARRIS_RATIO_CONSTRAINT_ADDITIVE_TOLERANCE = 0
 const double GlobalConfiguration::HARRIS_RATIO_CONSTRAINT_MULTIPLICATIVE_TOLERANCE = 0.001 * 0.0000001 * 0.5;
 const double GlobalConfiguration::BASIC_COSTS_ADDITIVE_TOLERANCE = 0.0000001;
 const double GlobalConfiguration::BASIC_COSTS_MULTIPLICATIVE_TOLERANCE = 0.001 * 0.0000001;
+const double GlobalConfiguration::SPARSE_FORREST_TOMLIN_DIAGONAL_ELEMENT_TOLERANCE = 0.00001;
 const unsigned GlobalConfiguration::DEGRADATION_CHECKING_FREQUENCY = 100;
 const double GlobalConfiguration::DEGRADATION_THRESHOLD = 0.1;
 const double GlobalConfiguration::ACCEPTABLE_SIMPLEX_PIVOT_THRESHOLD = 0.0001;

src/configuration/GlobalConfiguration.h

@@ -49,6 +49,9 @@ class GlobalConfiguration
     static const double BASIC_COSTS_ADDITIVE_TOLERANCE;
     static const double BASIC_COSTS_MULTIPLICATIVE_TOLERANCE;
 
+    // Sparse ForrestTomlin diagonal element tolerance constant
+    static const double SPARSE_FORREST_TOMLIN_DIAGONAL_ELEMENT_TOLERANCE;
+
     // Toggle use of Harris' two-pass ratio test for selecting the leaving variable
     static const bool USE_HARRIS_RATIO_TEST;
 

src/engine/Makefile

@@ -22,8 +22,7 @@ CFLAGS = \
 	-Werror \
 	-Wno-deprecated \
 	-std=c++0x \
-	-O3
-	\
+	-O3 \
 
 include $(BASIS_FACTORIZATION_DIR)/Sources.mk
 include $(COMMON_DIR)/Sources.mk

src/engine/ReluplexError.h

@@ -24,7 +24,7 @@ class ReluplexError : public Error
         VARIABLE_DOESNT_EXIST_IN_SOLUTION = 2,
         PARTICIPATING_VARIABLES_ABSENT = 3,
         NON_EQUALITY_INPUT_EQUATIONS_NOT_YET_SUPPORTED = 4,
-        INVALID_BOUND_TIGHTENING = 5,
+        SIMULATOR_ERROR = 5,
         MISSING_PL_CONSTRAINT_STATE = 6,
         REQUESTED_CASE_SPLITS_FROM_FIXED_CONSTRAINT = 7,
         UNBOUNDED_VARIABLES_NOT_YET_SUPPORTED = 8,
@@ -36,7 +36,6 @@ 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,
     };

src/engine/Sources.mk

@@ -20,7 +20,6 @@ SOURCES += \
 	RowBoundTightener.cpp \
 	Simulator.cpp \
 	SmtCore.cpp \
-	Statistics.cpp \
 	Tableau.cpp \
 	TableauRow.cpp \
 	TableauState.cpp \

src/engine/Tableau.cpp

@@ -294,6 +294,7 @@ void Tableau::setDimensions( unsigned m, unsigned n )
     _basisFactorization = BasisFactorizationFactory::createBasisFactorization( _m, *this );
     if ( !_basisFactorization )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::basisFactorization" );
+    _basisFactorization->setStatistics( _statistics );
 
     _workM = new double[m];
     if ( !_workM )