New bound tightening technique + bug fix (#38)

* new way for performing bound tightening

* cleanup. Control the 3 kinds of explicit basis bound tightenings from
the configuration file. Remove some code duplications

* bug fix

* undo change

src/configuration/GlobalConfiguration.cpp

@@ -33,7 +33,10 @@ const bool GlobalConfiguration::PREPROCESSOR_PL_CONSTRAINTS_ADD_AUX_EQUATIONS =
 const unsigned GlobalConfiguration::PSE_ITERATIONS_BEFORE_RESET = 1000;
 const double GlobalConfiguration::PSE_GAMMA_ERROR_THRESHOLD = 0.001;
 
-const bool GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS = true;
+const GlobalConfiguration::ExplicitBasisBoundTighteningType GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE =
+    GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX;
+const bool GlobalConfiguration::EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION = false;
+
 const unsigned GlobalConfiguration::REFACTORIZATION_THRESHOLD = 100;
 const GlobalConfiguration::BasisFactorizationType GlobalConfiguration::BASIS_FACTORIZATION_TYPE =
     GlobalConfiguration::LU_FACTORIZATION;
@@ -69,8 +72,30 @@ void GlobalConfiguration::print()
             PREPROCESSOR_PL_CONSTRAINTS_ADD_AUX_EQUATIONS ? "Yes" : "No" );
     printf( "  PSE_ITERATIONS_BEFORE_RESET: %u\n", PSE_ITERATIONS_BEFORE_RESET );
     printf( "  PSE_GAMMA_ERROR_THRESHOLD: %.15lf\n", PSE_GAMMA_ERROR_THRESHOLD );
-    printf( "  EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS: %s\n",
-            EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS ? "Yes" : "No" );
+
+    String basisBoundTighteningType;
+    switch ( EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE )
+    {
+    case USE_BASIS_MATRIX:
+        basisBoundTighteningType = "Use basis matrix";
+        break;
+
+    case COMPUTE_INVERTED_BASIS_MATRIX:
+        basisBoundTighteningType = "Compute inverted basis matrix";
+        break;
+
+    case USE_IMPLICIT_INVERTED_BASIS_MATRIX:
+        basisBoundTighteningType = "Use implicit inverted basis matrix";
+        break;
+
+    default:
+        basisBoundTighteningType = "Unknown";
+        break;
+    }
+
+    printf( "  EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS: %s\n", basisBoundTighteningType.ascii() );
+    printf( "  EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION: %s\n",
+            EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION ? "Yes" : "No" );
     printf( "  REFACTORIZATION_THRESHOLD: %u\n", REFACTORIZATION_THRESHOLD );
 
     String basisFactorizationType;

src/configuration/GlobalConfiguration.h

@@ -71,11 +71,27 @@ class GlobalConfiguration
     static const double PSE_GAMMA_ERROR_THRESHOLD;
 
     /*
-      Basis factorization options
+      Bound tightening options
     */
 
+    enum ExplicitBasisBoundTighteningType {
+        // Use the basis matrix without inverting it
+        USE_BASIS_MATRIX = 0,
+        // Compute the inverse basis matrix and use it
+        COMPUTE_INVERTED_BASIS_MATRIX = 1,
+        // Use the inverted basis matrix without using it, via transformations
+        USE_IMPLICIT_INVERTED_BASIS_MATRIX = 2,
+    };
+
     // When doing bound tightening using the explicit basis matrix, should the basis matrix be inverted?
-    static const bool EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS;
+    static const ExplicitBasisBoundTighteningType EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE;
+
+    // When doing explicit bound tightening, should we repeat until saturation?
+    static const bool EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION;
+
+    /*
+      Basis factorization options
+    */
 
     // The number of accumualted eta matrices, after which the basis will be refactorized
 	static const unsigned REFACTORIZATION_THRESHOLD;

src/engine/Engine.cpp

@@ -414,7 +414,7 @@ void Engine::performSimplexStep()
 
     // Tighten
     if ( !fakePivot &&
-         FloatUtils::lt( bestPivotEntry, GlobalConfiguration::ACCEPTABLE_SIMPLEX_PIVOT_THRESHOLD ) )
+         FloatUtils::gte( bestPivotEntry, GlobalConfiguration::ACCEPTABLE_SIMPLEX_PIVOT_THRESHOLD ) )
         _rowBoundTightener->examinePivotRow( _tableau );
 
     struct timespec end = TimeUtils::sampleMicro();
@@ -894,10 +894,22 @@ void Engine::explicitBasisBoundTightening()
 {
     struct timespec start = TimeUtils::sampleMicro();
 
-    if ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_INVERT_BASIS )
-        _rowBoundTightener->examineInvertedBasisMatrix( _tableau, false );
-    else
-        _rowBoundTightener->examineBasisMatrix( _tableau, false );
+    bool saturation = GlobalConfiguration::EXPLICIT_BOUND_TIGHTENING_UNTIL_SATURATION;
+
+    switch ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE )
+    {
+    case GlobalConfiguration::USE_BASIS_MATRIX:
+        _rowBoundTightener->examineBasisMatrix( _tableau, saturation );
+        break;
+
+    case GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX:
+        _rowBoundTightener->examineInvertedBasisMatrix( _tableau, saturation );
+        break;
+
+    case GlobalConfiguration::USE_IMPLICIT_INVERTED_BASIS_MATRIX:
+        _rowBoundTightener->examineImplicitInvertedBasisMatrix( _tableau, saturation );
+        break;
+    }
 
     struct timespec end = TimeUtils::sampleMicro();
     _statistics.addTimeForExplicitBasisBoundTightening( TimeUtils::timePassed( start, end ) );

src/engine/IRowBoundTightener.h

@@ -48,6 +48,15 @@ class IRowBoundTightener : public ITableau::VariableWatcher
      */
     virtual void examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ) = 0;
 
+    /*
+      Derive and enqueue new bounds for all varaibles, implicitly using the
+      inverse of the explicit basis matrix, inv(B0), which should be available
+      through the tableau. Inv(B0) is not computed directly --- instead, the computation
+      is performed via FTRANs. Can also do this until saturation, meaning that we
+      continue until no new bounds are learned.
+    */
+    virtual void examineImplicitInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ) = 0;
+
     /*
       Derive and enqueue new bounds for all varaibles, using the
       original constraint matrix A and right hands side vector b. Can

src/engine/RowBoundTightener.cpp

@@ -21,6 +21,8 @@ RowBoundTightener::RowBoundTightener()
     , _upperBounds( NULL )
     , _tightenedLower( NULL )
     , _tightenedUpper( NULL )
+    , _rows( NULL )
+    , _z( NULL )
     , _statistics( NULL )
 {
 }
@@ -30,6 +32,7 @@ void RowBoundTightener::initialize( const ITableau &tableau )
     freeMemoryIfNeeded();
 
     _n = tableau.getN();
+    _m = tableau.getM();
 
     _lowerBounds = new double[_n];
     if ( !_lowerBounds )
@@ -55,6 +58,23 @@ void RowBoundTightener::initialize( const ITableau &tableau )
         _lowerBounds[i] = tableau.getLowerBound( i );
         _upperBounds[i] = tableau.getUpperBound( i );
     }
+
+    if ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE ==
+         GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX )
+    {
+        _rows = new TableauRow *[_m];
+        for ( unsigned i = 0; i < _m; ++i )
+            _rows[i] = new TableauRow( _n - _m );
+    }
+    else if ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE ==
+              GlobalConfiguration::USE_IMPLICIT_INVERTED_BASIS_MATRIX )
+    {
+        _rows = new TableauRow *[_m];
+        for ( unsigned i = 0; i < _m; ++i )
+            _rows[i] = new TableauRow( _n - _m );
+
+        _z = new double[_m];
+    }
 }
 
 void RowBoundTightener::clear( const ITableau &tableau )
@@ -99,6 +119,62 @@ void RowBoundTightener::freeMemoryIfNeeded()
         delete[] _tightenedUpper;
         _tightenedUpper = NULL;
     }
+
+    if ( _rows )
+    {
+        for ( unsigned i = 0; i < _m; ++i )
+            delete _rows[i];
+        delete[] _rows;
+        _rows = NULL;
+    }
+
+    if ( _z )
+    {
+        delete[] _z;
+        _z = NULL;
+    }
+}
+
+void RowBoundTightener::examineImplicitInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation )
+{
+    /*
+      Roughly (the dimensions don't add up):
+
+         xB = inv(B)*b - inv(B)*An
+    */
+
+    // Find z = inv(B) * b, by solving the forward transformation Bz = b
+    tableau.forwardTransformation( tableau.getRightHandSide(), _z );
+    for ( unsigned i = 0; i < _m; ++i )
+    {
+        _rows[i]->_scalar = _z[i];
+        _rows[i]->_lhs = tableau.basicIndexToVariable( i );
+    }
+
+    // Now, go over the columns of the constraint martrix, perform an FTRAN
+    // for each of them, and populate the rows.
+    for ( unsigned i = 0; i < _n - _m; ++i )
+    {
+        unsigned nonBasic = tableau.nonBasicIndexToVariable( i );
+        const double *ANColumn = tableau.getAColumn( nonBasic );
+        tableau.forwardTransformation( ANColumn, _z );
+
+        for ( unsigned j = 0; j < _m; ++j )
+        {
+            _rows[j]->_row[i]._var = nonBasic;
+            _rows[j]->_row[i]._coefficient = -_z[j];
+        }
+    }
+
+    // We now have all the rows, can use them for tightening.
+    // The tightening procedure may throw an exception, in which case we need
+    // to release the rows.
+    bool newBoundsLearned;
+    do
+    {
+        newBoundsLearned = onePassOverInvertedBasisRows( tableau );
+    }
+    while ( untilSaturation && newBoundsLearned );
 }
 
 void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation )
@@ -111,73 +187,62 @@ void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, boo
       We compute one row at a time.
     */
 
-    unsigned n = tableau.getN();
-    unsigned m = tableau.getM();
-    List<TableauRow *> rows;
-
     const double *b = tableau.getRightHandSide();
     const double *invB = tableau.getInverseBasisMatrix();
 
     try
     {
-        for ( unsigned i = 0; i < m; ++i )
+        for ( unsigned i = 0; i < _m; ++i )
         {
-            TableauRow *row = new TableauRow( n - m );
+            TableauRow *row = _rows[i];
             // First, compute the scalar, using inv(B)*b
             row->_scalar = 0;
-            for ( unsigned j = 0; j < m; ++j )
-                row->_scalar += ( invB[i*m + j] * b[j] );
+            for ( unsigned j = 0; j < _m; ++j )
+                row->_scalar += ( invB[i * _m + j] * b[j] );
 
             // Now update the row's coefficients for basic variable i
-            for ( unsigned j = 0; j < n - m; ++j )
+            for ( unsigned j = 0; j < _n - _m; ++j )
             {
                 row->_row[j]._var = tableau.nonBasicIndexToVariable( j );
 
                 // Dot product of the i'th row of inv(B) with the appropriate
                 // column of An
                 const double *ANColumn = tableau.getAColumn( row->_row[j]._var );
                 row->_row[j]._coefficient = 0;
-                for ( unsigned k = 0; k < m; ++k )
-                    row->_row[j]._coefficient -= ( invB[i*m + k] * ANColumn[k] );
+                for ( unsigned k = 0; k < _m; ++k )
+                    row->_row[j]._coefficient -= ( invB[i * _m + k] * ANColumn[k] );
             }
 
             // Store the lhs variable
             row->_lhs = tableau.basicIndexToVariable( i );
-
-            // The row is ready
-            rows.append( row );
         }
 
         // We now have all the rows, can use them for tightening.
         // The tightening procedure may throw an exception, in which case we need
         // to release the rows.
+
         bool newBoundsLearned;
         do
         {
-            newBoundsLearned = onePassOverInvertedBasisRows( tableau, rows );
+            newBoundsLearned = onePassOverInvertedBasisRows( tableau );
         }
         while ( untilSaturation && newBoundsLearned );
     }
     catch ( ... )
     {
-        for ( const auto &row : rows )
-            delete row;
         delete[] invB;
-
         throw;
     }
 
-    for ( const auto &row : rows )
-        delete row;
     delete[] invB;
 }
 
-bool RowBoundTightener::onePassOverInvertedBasisRows( const ITableau &tableau, List<TableauRow *> &rows )
+bool RowBoundTightener::onePassOverInvertedBasisRows( const ITableau &tableau )
 {
     bool result = false;
 
-    for ( const auto &row : rows )
-        if ( tightenOnSingleInvertedBasisRow( tableau, *row ) )
+    for ( unsigned i = 0; i < _m; ++i )
+        if ( tightenOnSingleInvertedBasisRow( tableau, *( _rows[i] ) ) )
             result = true;
 
     return result;

src/engine/RowBoundTightener.h

@@ -60,6 +60,15 @@ class RowBoundTightener : public IRowBoundTightener
      */
     void examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation );
 
+    /*
+      Derive and enqueue new bounds for all varaibles, implicitly using the
+      inverse of the explicit basis matrix, inv(B0), which should be available
+      through the tableau. Inv(B0) is not computed directly --- instead, the computation
+      is performed via FTRANs. Can also do this until saturation, meaning that we
+      continue until no new bounds are learned.
+     */
+    void examineImplicitInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation );
+
     /*
       Derive and enqueue new bounds for all varaibles, using the
       original constraint matrix A and right hands side vector b. Can
@@ -87,6 +96,7 @@ class RowBoundTightener : public IRowBoundTightener
 
 private:
     unsigned _n;
+    unsigned _m;
 
     /*
       Work space for the tightener to derive tighter bounds. These
@@ -99,6 +109,12 @@ class RowBoundTightener : public IRowBoundTightener
     bool *_tightenedLower;
     bool *_tightenedUpper;
 
+    /*
+      Work space for the inverted basis matrix tighteners
+    */
+    TableauRow **_rows;
+    double *_z;
+
     /*
       Statistics collection
     */
@@ -145,7 +161,7 @@ class RowBoundTightener : public IRowBoundTightener
       tighter bounds. Return true if new bounds are learned, false
       otherwise.
     */
-    bool onePassOverInvertedBasisRows( const ITableau &tableau, List<TableauRow *> &rows );
+    bool onePassOverInvertedBasisRows( const ITableau &tableau );
 
     /*
       Process the inverted basis row and attempt to derive tighter

src/engine/tests/MockRowBoundTightener.h

@@ -61,6 +61,7 @@ class MockRowBoundTightener : public IRowBoundTightener
     void examinePivotRow( ITableau &/* tableau */ ) {}
     void getRowTightenings( List<Tightening> &/* tightenings */ ) const {}
     void setStatistics( Statistics */* statistics */ ) {}
+    void examineImplicitInvertedBasisMatrix( const ITableau &/* tableau */, bool /* untilSaturation */ ) {}
 };
 
 #endif // __MockRowBoundTightener_h__