From f318b4f06416358b36d9c3460252ee1550794da7 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 8 May 2018 11:11:10 +0300 Subject: [PATCH 01/14] a more efficient computation of the explicit basis: 1. Compute from left to write instead of from right to left. This makes multipliciation by eta matrices easier 2. Compute invLP once-and-for-all, when the basis is factorized --- .../ForrestTomlinFactorization.cpp | 158 ++++++++++++------ .../ForrestTomlinFactorization.h | 8 + .../tests/Test_CompareFactorizations.h | 11 ++ .../tests/Test_ForrestTomlinFactorization.h | 42 +++++ 4 files changed, 170 insertions(+), 49 deletions(-) diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp index ddf760131..eae387d00 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.cpp +++ b/src/basis_factorization/ForrestTomlinFactorization.cpp @@ -24,6 +24,7 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m ) , _Q( m ) , _invQ( m ) , _U( NULL ) + , _invLP( NULL ) , _explicitBasisAvailable( false ) , _workMatrix( NULL ) , _workVector( NULL ) @@ -42,6 +43,11 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m ) throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, "ForrestTomlinFactorization::U" ); + _invLP = new double[m * m]; + if ( !_invLP ) + throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, + "ForrestTomlinFactorization::invLPx" ); + for ( unsigned i = 0; i < _m; ++i ) { _U[i] = new EtaMatrix( _m, i ); @@ -69,6 +75,18 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m ) if ( !_workOrdering ) throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, "ForrestTomlinFactorization::workOrdering" ); + + // Initialization: assume B = I + std::fill_n( _B, _m * _m, 0.0 ); + for ( unsigned row = 0; row < _m; ++row ) + _B[row * _m + row] = 1.0; + + // Also, invLP = I + std::fill_n( _invLP, _m * _m, 0.0 ); + for ( unsigned row = 0; row < _m; ++row ) + _invLP[row * _m + row] = 1.0; + + // Us, Q and invQ are all initialized to the identity matrix anyway. } ForrestTomlinFactorization::~ForrestTomlinFactorization() @@ -94,7 +112,13 @@ ForrestTomlinFactorization::~ForrestTomlinFactorization() _U = NULL; } - if ( _workMatrix ) + if ( _invLP ) + { + delete[] _invLP; + _invLP = NULL; + } + + if ( _workMatrix ) { delete[] _workMatrix; _workMatrix = NULL; @@ -663,6 +687,56 @@ void ForrestTomlinFactorization::initialLUFactorization() for ( unsigned j = 0; j <= i; ++j ) u->_column[j] = _workMatrix[j * _m + i]; } + + // Compute inv(P1)inv(L1) ... inv(Ps)inv(LS) for later. + + // Initialize to the identity matrix + std::fill_n( _invLP, _m * _m, 0 ); + for ( unsigned i = 0; i < _m; ++i ) + _invLP[i*_m + i] = 1.0; + + // Multiply by inverted Ps and Ls + for ( auto lp = _LP.rbegin(); lp != _LP.rend(); ++lp ) + { + if ( (*lp)->_pair ) + { + // The inverse of a general permutation matrix is its transpose. + // However, since these are permutations of just two rows, the + // transpose is equal to the original - so no need to invert. + // Multiplying on the right permutes columns. + double temp; + unsigned first = (*lp)->_pair->first; + unsigned second = (*lp)->_pair->second; + + for ( unsigned i = 0; i < _m; ++i ) + { + temp = _invLP[i * _m + first]; + _invLP[i * _m + first] = _invLP[i * _m + second]; + _invLP[i * _m + second] = temp; + } + } + else + { + // The inverse of an eta matrix is an eta matrix with the special + // column negated and divided by the diagonal entry. The only + // exception is the diagonal entry itself, which is just the + // inverse of the original diagonal entry. + EtaMatrix *eta = (*lp)->_eta; + unsigned columnIndex = eta->_columnIndex; + double etaDiagonalEntry = 1 / eta->_column[columnIndex]; + + for ( unsigned row = 0; row < _m; ++row ) + { + for ( unsigned col = columnIndex + 1; col < _m; ++col ) + { + _invLP[row * _m + columnIndex] -= + ( eta->_column[col] * _invLP[row * _m + col] ); + } + + _invLP[row * _m + columnIndex] *= etaDiagonalEntry; + } + } + } } bool ForrestTomlinFactorization::explicitBasisAvailable() const @@ -686,76 +760,57 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable() * Q * Um...U1 * invQ */ - // Start by computing: Um...U1 * invQ. - std::fill_n( _workMatrix, _m * _m, 0.0 ); - for ( unsigned i = 0; i < _m; ++i ) - _workMatrix[i * _m + _invQ._ordering[i]] = 1.0; - - for ( unsigned i = 0; i < _m; ++i ) - { - for ( unsigned row = 0; row < _U[i]->_columnIndex; ++row ) - for ( unsigned col = 0; col < _m; ++col ) - _workMatrix[row * _m + col] += _U[i]->_column[row] * _workMatrix[_U[i]->_columnIndex * _m + col]; - } - - // Permute the rows according to Q - for ( unsigned i = 0; i < _m; ++i ) - memcpy( _B + ( i * _m ), _workMatrix + ( _Q._ordering[i] * _m ), sizeof(double) * _m ); + // Start with the already-computed inv(P1)inv(L1) ... inv(Ps)inv(LS) + memcpy( _B, _invLP, sizeof(double) * _m * _m ); // Multiply by the inverse As. An inverse of an almost identity matrix // is just that matrix with its special value negated (unless that // matrix is diagonal. - for ( auto a = _A.rbegin(); a != _A.rend(); ++a ) + for ( const auto &a : _A ) { - if ( (*a)->_row == (*a)->_column ) + if ( a->_row == a->_column ) { for ( unsigned j = 0; j < _m; ++j ) - _B[(*a)->_row * _m + j] *= ( 1 / (*a)->_value ); + _B[j * _m + a->_column] *= ( 1 / a->_value ); } else { for ( unsigned j = 0; j < _m; ++j ) - _B[(*a)->_row * _m + j] -= _B[(*a)->_column * _m + j] * (*a)->_value; + _B[j * _m + a->_column] -= _B[j * _m + a->_row] * a->_value; } } - // Multiply by inverted Ps and Ls - for ( const auto &lp : _LP ) + // Permute columns according to Q + for ( unsigned i = 0; i < _m; ++i ) { - // The inverse of a general permutation matrix is its transpose. - // However, since these are permutations of just two rows, the - // transpose is equal to the original - so no need to invert. - if ( lp->_pair ) - { - memcpy( _workVector, _B + (lp->_pair->first * _m), sizeof(double) * _m ); - memcpy( _B + (lp->_pair->first * _m), _B + (lp->_pair->second * _m), sizeof(double) * _m ); - memcpy( _B + (lp->_pair->second * _m), _workVector, sizeof(double) * _m ); - } - else - { - // The inverse of an eta matrix is an eta matrix with the special - // column negated and divided by the diagonal entry. The only - // exception is the diagonal entry itself, which is just the - // inverse of the original diagonal entry. - EtaMatrix *eta = lp->_eta; - double etaDiagonalEntry = 1 / eta->_column[eta->_columnIndex]; + unsigned columnToCopy = _invQ._ordering[i]; + for ( unsigned j = 0; j < _m; ++j ) + _workMatrix[j * _m + i] = _B[j * _m + columnToCopy]; + } - for ( unsigned row = eta->_columnIndex + 1; row < _m; ++row ) + // Multiply by Um...U1 + for ( int i = _m - 1; i >= 0; --i ) + { + for ( unsigned row = 0; row < _m; ++row ) + { + for ( unsigned j = 0; j < _U[i]->_columnIndex; ++j ) { - for ( unsigned col = 0; col < _m; ++col ) - { - _B[row * _m + col] -= - eta->_column[row] * etaDiagonalEntry * _B[eta->_columnIndex * _m + col]; - } + _workMatrix[row * _m + _U[i]->_columnIndex] += + ( _U[i]->_column[j] * _workMatrix[row * _m + j] ); } - - for ( unsigned col = 0; col < _m; ++col ) - _B[eta->_columnIndex * _m + col] *= etaDiagonalEntry; } } + // Permute columns according to invQ + for ( unsigned i = 0; i < _m; ++i ) + { + unsigned columnToCopy = _Q._ordering[i]; + for ( unsigned j = 0; j < _m; ++j ) + _B[j * _m + i] = _workMatrix[j * _m + columnToCopy]; + } + clearFactorization(); - initialLUFactorization(); + initialLUFactorization(); _explicitBasisAvailable = true; } @@ -932,6 +987,11 @@ void ForrestTomlinFactorization::dump() const printf( "*** Done dumping FT factorization ***\n\no" ); } +const double *ForrestTomlinFactorization::getInvLP() const +{ + return _invLP; +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/basis_factorization/ForrestTomlinFactorization.h b/src/basis_factorization/ForrestTomlinFactorization.h index 9f04d5e2f..149e13175 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.h +++ b/src/basis_factorization/ForrestTomlinFactorization.h @@ -97,6 +97,7 @@ class ForrestTomlinFactorization : public IBasisFactorization const EtaMatrix **getU() const; const List *getLP() const; const List *getA() const; + const double *getInvLP() const; void pushA( const AlmostIdentityMatrix &matrix ); void setQ( const PermutationMatrix &Q ); @@ -123,6 +124,13 @@ class ForrestTomlinFactorization : public IBasisFactorization // A1 is the first element of the list List_A; + /* + A helper matrix to facilitate computing the explicit basis. + This matrix stores: + inv(P1)inv(L1) ... inv(Ps)inv(LS) + */ + double *_invLP; + /* Indicates whether the explicit basis matrix is available. */ diff --git a/src/basis_factorization/tests/Test_CompareFactorizations.h b/src/basis_factorization/tests/Test_CompareFactorizations.h index 090986c33..4ada62db3 100644 --- a/src/basis_factorization/tests/Test_CompareFactorizations.h +++ b/src/basis_factorization/tests/Test_CompareFactorizations.h @@ -46,6 +46,17 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); TS_ASSERT( lu = new LUFactorization( 4 ) ); + double B[] = + { + 2, 0, 3, -4, + 0, 1, 10, 0, + -3, 4.5, 1, 1, + 0, 0, 2, 2 + }; + + TS_ASSERT_THROWS_NOTHING( ft->setBasis( B ) ); + TS_ASSERT_THROWS_NOTHING( lu->setBasis( B ) ); + double y[4] = { 9, 15, 10, -12 }; double x1[4]; double x2[4]; diff --git a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h index 0ef71dc31..5391e717b 100644 --- a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h +++ b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h @@ -188,6 +188,48 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite ++lIt; TS_ASSERT_EQUALS( *((*lIt)->_eta), expectedL1 ); + /* + L4 = | 1 0 0 0 | + | 0 1 0 0 | + | 0 0 1 0 | + | 0 0 0 -0.5 | + + L3 = | 1 0 0 0 | + | 0 1 0 0 | + | 0 0 -1 0 | + | 0 0 1 1 | + + L2 = | 1 0 0 0 | + | 0 0.5 0 0 | + | 0 0 1 0 | + | 0 0 0 1 | + + L1 = | 1 0 0 0 | + | -1 1 0 0 | + | -1 0 1 0 | + | 1 0 0 1 | + + InvLP = inv( L4L3L2L1 ) + + = inv( | 1 0 0 0 | ) + | -0.5 0.5 0 0 | + | 1 0 -1 0 | + | 0 0 -0.5 -0.5 | + + = | 1 0 0 0 | + | 1 2 0 0 | + | 1 0 -1 0 | + | -1 0 1 -2 | + */ + double expectedInvLP[] = + { 1, 0, 0, 0, + 1, 2, 0, 0, + 1, 0, -1, 0, + -1, 0, 1, - 2}; + + for ( unsigned i = 0; i < 16; ++i ) + TS_ASSERT( FloatUtils::areEqual( expectedInvLP[i], ft->getInvLP()[i] ) ); + TS_ASSERT_THROWS_NOTHING( delete ft ); } From 11796abb783a8f67e4710986660c286a8ece8f3c Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 8 May 2018 12:31:38 +0300 Subject: [PATCH 02/14] pivoting stats --- src/engine/Statistics.cpp | 11 ++++++++++- src/engine/Statistics.h | 4 ++++ src/engine/Tableau.cpp | 19 +++++++++++++++++++ 3 files changed, 33 insertions(+), 1 deletion(-) diff --git a/src/engine/Statistics.cpp b/src/engine/Statistics.cpp index 0ac77f663..239cad9b5 100644 --- a/src/engine/Statistics.cpp +++ b/src/engine/Statistics.cpp @@ -36,6 +36,7 @@ Statistics::Statistics() , _numTableauPivots( 0 ) , _numTableauDegeneratePivots( 0 ) , _numTableauDegeneratePivotsByRequest( 0 ) + , _timePivotsMicro( 0 ) , _numSimplexPivotSelectionsIgnoredForStability( 0 ) , _numSimplexUnstablePivots( 0 ) , _numTableauBoundHopping( 0 ) @@ -214,7 +215,10 @@ void Statistics::print() , _numTableauDegeneratePivotsByRequest , printPercents( _numTableauDegeneratePivotsByRequest, _numTableauDegeneratePivots ) ); - printf( "\t\tFake pivots: %llu\n", _numTableauBoundHopping ); + printf( "\t\tAverage time per pivot: %.2lf milli\n", + printAverage( _timePivotsMicro / 1000, _numTableauPivots ) ); + + printf( "\tTotal number of fake pivots performed: %llu\n", _numTableauBoundHopping ); printf( "\t--- SMT Core Statistics ---\n" ); printf( "\tTotal depth is %u. Total visited states: %u. Number of splits: %u. Number of pops: %u\n" @@ -368,6 +372,11 @@ void Statistics::incNumTableauDegeneratePivotsByRequest() ++_numTableauDegeneratePivotsByRequest; } +void Statistics::addTimePivots( unsigned long long time ) +{ + _timePivotsMicro += time; +} + void Statistics::incNumSimplexPivotSelectionsIgnoredForStability() { ++_numSimplexPivotSelectionsIgnoredForStability; diff --git a/src/engine/Statistics.h b/src/engine/Statistics.h index 9fd7ca350..6ebd1a35c 100644 --- a/src/engine/Statistics.h +++ b/src/engine/Statistics.h @@ -71,6 +71,7 @@ class Statistics void incNumTableauDegeneratePivotsByRequest(); void incNumSimplexPivotSelectionsIgnoredForStability(); void incNumSimplexUnstablePivots(); + void addTimePivots( unsigned long long time ); unsigned long long getNumTableauPivots() const; unsigned long long getNumSimplexPivotSelectionsIgnoredForStability() const; unsigned long long getNumSimplexUnstablePivots() const; @@ -175,6 +176,9 @@ class Statistics // by explicit request unsigned long long _numTableauDegeneratePivotsByRequest; + // Total time for performing pivots (both real and degenrate), in microseconds + unsigned long long _timePivotsMicro; + // Total number of entering/leaving variable pairs ignored because their pivot // element was too small unsigned long long _numSimplexPivotSelectionsIgnoredForStability; diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 7ce1820f5..d34cdd566 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -606,8 +606,13 @@ void Tableau::performPivot() return; } + struct timespec pivotStart; + if ( _statistics ) + { + pivotStart = TimeUtils::sampleMicro(); _statistics->incNumTableauPivots(); + } unsigned currentBasic = _basicIndexToVariable[_leavingVariable]; unsigned currentNonBasic = _nonBasicIndexToVariable[_enteringVariable]; @@ -649,12 +654,20 @@ void Tableau::performPivot() // Update the basis factorization. The column corresponding to the // leaving variable is the one that has changed _basisFactorization->pushEtaMatrix( _leavingVariable, _changeColumn ); + + if ( _statistics ) + { + struct timespec pivotEnd = TimeUtils::sampleMicro(); + _statistics->addTimePivots( TimeUtils::timePassed( pivotStart, pivotEnd ) ); + } } void Tableau::performDegeneratePivot() { + struct timespec pivotStart; if ( _statistics ) { + pivotStart = TimeUtils::sampleMicro(); _statistics->incNumTableauPivots(); _statistics->incNumTableauDegeneratePivots(); _statistics->incNumTableauDegeneratePivotsByRequest(); @@ -688,6 +701,12 @@ void Tableau::performDegeneratePivot() _basicAssignment[_leavingVariable] = _nonBasicAssignment[_enteringVariable]; _nonBasicAssignment[_enteringVariable] = temp; computeBasicStatus( _leavingVariable ); + + if ( _statistics ) + { + struct timespec pivotEnd = TimeUtils::sampleMicro(); + _statistics->addTimePivots( TimeUtils::timePassed( pivotStart, pivotEnd ) ); + } } double Tableau::ratioConstraintPerBasic( unsigned basicIndex, double coefficient, bool decrease ) From a32d35f18d96161b4c29c535ef1bde99dfa63270 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 8 May 2018 12:31:52 +0300 Subject: [PATCH 03/14] bug fix in store/restore --- src/basis_factorization/ForrestTomlinFactorization.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp index eae387d00..3c0d98635 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.cpp +++ b/src/basis_factorization/ForrestTomlinFactorization.cpp @@ -533,6 +533,8 @@ void ForrestTomlinFactorization::storeFactorization( IBasisFactorization *other for ( const auto &lp : _LP ) otherFTFactorization->_LP.append( lp->duplicate() ); + memcpy( otherFTFactorization->_invLP, _invLP, sizeof(double) * _m * _m ); + List::iterator aIt; for ( aIt = otherFTFactorization->_A.begin(); aIt != otherFTFactorization->_A.end(); ++aIt ) delete *aIt; @@ -570,6 +572,8 @@ void ForrestTomlinFactorization::restoreFactorization( const IBasisFactorization for ( const auto &lp : otherFTFactorization->_LP ) _LP.append( lp->duplicate() ); + memcpy( _invLP, otherFTFactorization->_invLP, sizeof(double) * _m * _m ); + List::iterator aIt; for ( aIt = _A.begin(); aIt != _A.end(); ++aIt ) delete *aIt; From e53123c966ce02b3e1acd8aa62f82a0d6d40037c Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 9 May 2018 09:56:08 +0300 Subject: [PATCH 04/14] supprt for a "basis column oracle", which allows the basis factorization classes to query the tableau for constraint matrix columns --- .../BasisFactorizationFactory.cpp | 6 +-- .../BasisFactorizationFactory.h | 2 +- .../ForrestTomlinFactorization.cpp | 27 +++++++++++-- .../ForrestTomlinFactorization.h | 7 +++- src/basis_factorization/IBasisFactorization.h | 16 +++++++- src/basis_factorization/LUFactorization.cpp | 5 ++- src/basis_factorization/LUFactorization.h | 2 +- .../tests/MockColumnOracle.h | 38 +++++++++++++++++++ .../tests/Test_CompareFactorizations.h | 8 +++- .../tests/Test_ForrestTomlinFactorization.h | 28 ++++++++------ .../tests/Test_LUFactorization.h | 38 ++++++++++--------- src/engine/Tableau.cpp | 13 +++++-- src/engine/Tableau.h | 6 ++- src/engine/TableauState.cpp | 4 +- src/engine/TableauState.h | 5 +-- 15 files changed, 152 insertions(+), 53 deletions(-) create mode 100644 src/basis_factorization/tests/MockColumnOracle.h diff --git a/src/basis_factorization/BasisFactorizationFactory.cpp b/src/basis_factorization/BasisFactorizationFactory.cpp index 2352863c9..f87db4367 100644 --- a/src/basis_factorization/BasisFactorizationFactory.cpp +++ b/src/basis_factorization/BasisFactorizationFactory.cpp @@ -16,14 +16,14 @@ #include "GlobalConfiguration.h" #include "LUFactorization.h" -IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize ) +IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle ) { if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE == GlobalConfiguration::LU_FACTORIZATION ) - return new LUFactorization( basisSize ); + return new LUFactorization( basisSize, basisColumnOracle ); else if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE == GlobalConfiguration::FORREST_TOMLIN_FACTORIZATION ) - return new ForrestTomlinFactorization( basisSize ); + return new ForrestTomlinFactorization( basisSize, basisColumnOracle ); throw BasisFactorizationError( BasisFactorizationError::UNKNOWN_BASIS_FACTORIZATION_TYPE ); } diff --git a/src/basis_factorization/BasisFactorizationFactory.h b/src/basis_factorization/BasisFactorizationFactory.h index 534a5be33..65c699b88 100644 --- a/src/basis_factorization/BasisFactorizationFactory.h +++ b/src/basis_factorization/BasisFactorizationFactory.h @@ -18,7 +18,7 @@ class BasisFactorizationFactory { public: - static IBasisFactorization *createBasisFactorization( unsigned basisSize ); + static IBasisFactorization *createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle ); }; #endif // __BasisFactorizationFactory_h__ diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp index 3c0d98635..0bc950526 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.cpp +++ b/src/basis_factorization/ForrestTomlinFactorization.cpp @@ -18,8 +18,9 @@ #include #include -ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m ) - : _m( m ) +ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle ) + : IBasisFactorization( basisColumnOracle ) + , _m( m ) , _B( NULL ) , _Q( m ) , _invQ( m ) @@ -309,7 +310,9 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub // If the number of A matrices is too great, condense them. if ( _A.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD ) - makeExplicitBasisAvailable(); + refactorizeBasis(); + + // makeExplicitBasisAvailable(); } void ForrestTomlinFactorization::forwardTransformation( const double *y, double *x ) const @@ -752,6 +755,10 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable() { if ( explicitBasisAvailable() ) return; + + refactorizeBasis(); + return; + /* We know that the following equation holds: @@ -996,6 +1003,20 @@ const double *ForrestTomlinFactorization::getInvLP() const return _invLP; } +void ForrestTomlinFactorization::refactorizeBasis() +{ + for ( unsigned column = 0; column < _m; ++column ) + { + const double *basisColumn = _basisColumnOracle->getColumnOfBasis( column ); + for ( unsigned row = 0; row < _m; ++row ) + _B[row * _m + column] = basisColumn[row]; + } + + clearFactorization(); + initialLUFactorization(); + _explicitBasisAvailable = true; +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/basis_factorization/ForrestTomlinFactorization.h b/src/basis_factorization/ForrestTomlinFactorization.h index 149e13175..1f22d21ef 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.h +++ b/src/basis_factorization/ForrestTomlinFactorization.h @@ -33,7 +33,7 @@ class ForrestTomlinFactorization : public IBasisFactorization { public: - ForrestTomlinFactorization( unsigned m ); + ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle ); ~ForrestTomlinFactorization(); /* @@ -88,6 +88,11 @@ class ForrestTomlinFactorization : public IBasisFactorization */ void invertBasis( double *result ); + /* + Obtain the basis matrix from the oracle and compute a fresh factorization + */ + void refactorizeBasis(); + public: /* For testing purposes only diff --git a/src/basis_factorization/IBasisFactorization.h b/src/basis_factorization/IBasisFactorization.h index fc0f90ca8..d6d7b1d66 100644 --- a/src/basis_factorization/IBasisFactorization.h +++ b/src/basis_factorization/IBasisFactorization.h @@ -19,8 +19,19 @@ class IBasisFactorization This is the interfact class for a basis factorization. */ public: - IBasisFactorization() + /* + A callback for obtaining columns of the basis matrix + */ + class BasisColumnOracle + { + public: + virtual ~BasisColumnOracle() {} + virtual const double *getColumnOfBasis( unsigned column ) const = 0; + }; + + IBasisFactorization( const BasisColumnOracle &basisColumnOracle ) : _factorizationEnabled( true ) + , _basisColumnOracle( &basisColumnOracle ) { } @@ -96,6 +107,9 @@ class IBasisFactorization disabled. */ bool _factorizationEnabled; + +protected: + const BasisColumnOracle *_basisColumnOracle; }; #endif // __IBasisFactorization_h__ diff --git a/src/basis_factorization/LUFactorization.cpp b/src/basis_factorization/LUFactorization.cpp index 507dd448c..d6c0e1153 100644 --- a/src/basis_factorization/LUFactorization.cpp +++ b/src/basis_factorization/LUFactorization.cpp @@ -21,8 +21,9 @@ #include "MStringf.h" #include "MalformedBasisException.h" -LUFactorization::LUFactorization( unsigned m ) - : _B0( NULL ) +LUFactorization::LUFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle ) + : IBasisFactorization( basisColumnOracle ) + , _B0( NULL ) , _m( m ) , _U( NULL ) , _tempY( NULL ) diff --git a/src/basis_factorization/LUFactorization.h b/src/basis_factorization/LUFactorization.h index cb2d05a0a..6f506fa6a 100644 --- a/src/basis_factorization/LUFactorization.h +++ b/src/basis_factorization/LUFactorization.h @@ -25,7 +25,7 @@ class LPElement; class LUFactorization : public IBasisFactorization { public: - LUFactorization( unsigned m ); + LUFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle ); ~LUFactorization(); /* diff --git a/src/basis_factorization/tests/MockColumnOracle.h b/src/basis_factorization/tests/MockColumnOracle.h new file mode 100644 index 000000000..fce3e9613 --- /dev/null +++ b/src/basis_factorization/tests/MockColumnOracle.h @@ -0,0 +1,38 @@ +/********************* */ +/*! \file MockColumnOracle.h +** \verbatim +** Top contributors (to current version): +** Guy Katz +** This file is part of the Marabou project. +** Copyright (c) 2016-2017 by the authors listed in the file AUTHORS +** in the top-level source directory) and their institutional affiliations. +** All rights reserved. See the file COPYING in the top-level source +** directory for licensing information.\endverbatim +**/ + +#ifndef __MockColumnOracle_h__ +#define __MockColumnOracle_h__ + +#include "IBasisFactorization.h" + +class MockColumnOracle : public IBasisFactorization::BasisColumnOracle +{ +public: + const double *_basis; + unsigned _m; + + const double *getColumnOfBasis( unsigned column ) const + { + return _basis + ( _m * column ); + } +}; + +#endif // __MockColumnOracle_h__ + +// +// Local Variables: +// compile-command: "make -C ../../.. " +// tags-file-name: "../../../TAGS" +// c-basic-offset: 4 +// End: +// diff --git a/src/basis_factorization/tests/Test_CompareFactorizations.h b/src/basis_factorization/tests/Test_CompareFactorizations.h index 4ada62db3..0b4c0ab63 100644 --- a/src/basis_factorization/tests/Test_CompareFactorizations.h +++ b/src/basis_factorization/tests/Test_CompareFactorizations.h @@ -16,6 +16,7 @@ #include "FloatUtils.h" #include "ForrestTomlinFactorization.h" #include "LUFactorization.h" +#include "MockColumnOracle.h" #include "MockErrno.h" class MockForCompareFactorizations @@ -27,14 +28,17 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite { public: MockForCompareFactorizations *mock; + MockColumnOracle *oracle; void setUp() { TS_ASSERT( mock = new MockForCompareFactorizations ); + TS_ASSERT( oracle = new MockColumnOracle ); } void tearDown() { + TS_ASSERT_THROWS_NOTHING( delete oracle ); TS_ASSERT_THROWS_NOTHING( delete mock ); } @@ -43,8 +47,8 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite ForrestTomlinFactorization *ft; LUFactorization *lu; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); - TS_ASSERT( lu = new LUFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); + TS_ASSERT( lu = new LUFactorization( 4, *oracle ) ); double B[] = { diff --git a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h index 5391e717b..839520969 100644 --- a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h +++ b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h @@ -16,9 +16,10 @@ #include "BasisFactorizationError.h" #include "EtaMatrix.h" #include "FloatUtils.h" -#include "GlobalConfiguration.h" #include "ForrestTomlinFactorization.h" +#include "GlobalConfiguration.h" #include "List.h" +#include "MockColumnOracle.h" #include "MockErrno.h" class MockForForrestTomlinFactorization @@ -30,6 +31,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { public: MockForForrestTomlinFactorization *mock; + MockColumnOracle *oracle; bool isIdentityPermutation( const PermutationMatrix *matrix ) { @@ -43,10 +45,12 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite void setUp() { TS_ASSERT( mock = new MockForForrestTomlinFactorization ); + TS_ASSERT( oracle = new MockColumnOracle ); } void tearDown() { + TS_ASSERT_THROWS_NOTHING( delete oracle ); TS_ASSERT_THROWS_NOTHING( delete mock ); } @@ -54,7 +58,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 3 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 3, *oracle ) ); TS_ASSERT( ft->factorizationEnabled() ); @@ -73,7 +77,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); double basisMatrix[16] = { 1, 3, -2, 4, @@ -237,7 +241,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); double basisMatrix[16] = { 1, 4, -2, 4, @@ -373,7 +377,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); double basisMatrix[16] = { 1, 3, -2, 4, @@ -501,7 +505,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); double basisMatrix[16] = { 1, 3, -2, 4, @@ -640,7 +644,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); // B = | 1 3 -2 4 | // | 1 5 -1 5 | @@ -854,7 +858,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); double basisMatrix[16] = { 1, 3, -2, 4, @@ -868,8 +872,8 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite double a1[] = { -4, 2, 0, 3 }; ft->pushEtaMatrix( 1, a1 ); - ForrestTomlinFactorization *ft2 = new ForrestTomlinFactorization( 4 ); - ForrestTomlinFactorization *ft3 = new ForrestTomlinFactorization( 4 ); + ForrestTomlinFactorization *ft2 = new ForrestTomlinFactorization( 4, *oracle ); + ForrestTomlinFactorization *ft3 = new ForrestTomlinFactorization( 4, *oracle ); ft->storeFactorization( ft2 ); ft3->restoreFactorization( ft2 ); @@ -894,7 +898,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); // B = | 1 3 -2 4 | // | 1 5 -1 5 | @@ -951,7 +955,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite { ForrestTomlinFactorization *ft; - TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) ); + TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) ); // B = | 1 3 -2 4 | // | 1 5 -1 5 | diff --git a/src/basis_factorization/tests/Test_LUFactorization.h b/src/basis_factorization/tests/Test_LUFactorization.h index 182cdd910..d63bc9e82 100644 --- a/src/basis_factorization/tests/Test_LUFactorization.h +++ b/src/basis_factorization/tests/Test_LUFactorization.h @@ -19,6 +19,7 @@ #include "GlobalConfiguration.h" #include "LUFactorization.h" #include "List.h" +#include "MockColumnOracle.h" #include "MockErrno.h" void matrixMultiply( unsigned dimension, const double *left, const double *right, double *result ) @@ -45,14 +46,17 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite { public: MockForLUFactorization *mock; + MockColumnOracle *oracle; void setUp() { TS_ASSERT( mock = new MockForLUFactorization ); + TS_ASSERT( oracle = new MockColumnOracle ); } void tearDown() { + TS_ASSERT_THROWS_NOTHING( delete oracle ); TS_ASSERT_THROWS_NOTHING( delete mock ); } @@ -60,7 +64,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite { LUFactorization *basis; - TS_ASSERT( basis = new LUFactorization( 3 ) ); + TS_ASSERT( basis = new LUFactorization( 3, *oracle ) ); TS_ASSERT( basis->factorizationEnabled() ); @@ -77,7 +81,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_forward_transformation() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); // If no eta matrices are provided, d = a double a1[] = { 1, 1, 3 }; @@ -127,7 +131,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_forward_transformation_with_B0() { // Same etas as test_backward_transformation() - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); double e1[] = {1., 1., 3.}; basis.pushEtaMatrix( 1, e1 ); double e2[] = {2., 1., 1.}; @@ -152,7 +156,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_backward_transformation() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); // If no eta matrices are provided, x = y double y1[] = { 1, 2, 3 }; @@ -222,7 +226,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_backward_transformation_2() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); // E1 = | -1 | // | 0 1 | @@ -247,7 +251,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_backward_transformation_with_B0() { // Same etas as test_backward_transformation() - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); double e1[] = {1., 1., 3.}; basis.pushEtaMatrix( 1, e1 ); double e2[] = {2., 1., 1.}; @@ -274,8 +278,8 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_store_and_restore() { - LUFactorization basis( 3 ); - LUFactorization otherBasis( 3 ); + LUFactorization basis( 3, *oracle ); + LUFactorization otherBasis( 3, *oracle ); double a1[] = { 1, 1, 3 }; double d1[] = { 0, 0, 0 }; @@ -327,7 +331,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_factorization_pivot() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); const int nsq = 9; double A[nsq]= {0., 1. , 0., @@ -362,7 +366,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_factorization_textbook() { // Textbook example - LUFactorization basis( 4 ); + LUFactorization basis( 4, *oracle ); const int nsq = 16; double A[nsq]= {1., 3., -2., 4., 1., 5., -1., 5., @@ -438,7 +442,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_factorization_as_black_box() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); const int nsq = 9; double A[nsq]= { 1., 2., 4., 4., 5., 7., @@ -472,7 +476,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_factorization_numerical_stability() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); double A[] = { @@ -555,8 +559,8 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite // TODO: this test fails when the REFACTORIZATION_THRESHOLD is too great (> 10 or so). // Disabling for now. - LUFactorization basis( 3 ); - LUFactorization basis2( 3 ); + LUFactorization basis( 3, *oracle ); + LUFactorization basis2( 3, *oracle ); basis.toggleFactorization( false ); int d = 3; unsigned etaCount = GlobalConfiguration::REFACTORIZATION_THRESHOLD + 2; @@ -633,7 +637,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_invert_B0_fail() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); double B0[] = { 1, 0, 0, 0, 1, 0, @@ -654,7 +658,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite void test_invert_B0() { - LUFactorization basis( 3 ); + LUFactorization basis( 3, *oracle ); { double B0[] = { 1, 0, 0, @@ -761,7 +765,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite } { - LUFactorization basis( 4 ); + LUFactorization basis( 4, *oracle ); double B0[] = { 1, 1, 1, 0, 0, 3, 1, 2, diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index d34cdd566..91ba6988a 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -218,7 +218,7 @@ void Tableau::setDimensions( unsigned m, unsigned n ) if ( !_basicStatus ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::basicStatus" ); - _basisFactorization = BasisFactorizationFactory::createBasisFactorization( _m ); + _basisFactorization = BasisFactorizationFactory::createBasisFactorization( _m, *this ); if ( !_basisFactorization ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::basisFactorization" ); @@ -1089,7 +1089,7 @@ void Tableau::dumpEquations() void Tableau::storeState( TableauState &state ) const { // Set the dimensions - state.setDimensions( _m, _n ); + state.setDimensions( _m, _n, *this ); // Store matrix A memcpy( state._A, _A, sizeof(double) * _n * _m ); @@ -1431,7 +1431,7 @@ void Tableau::addRow() // Allocate a larger basis factorization IBasisFactorization *newBasisFactorization = - BasisFactorizationFactory::createBasisFactorization( newM ); + BasisFactorizationFactory::createBasisFactorization( newM, *this ); if ( !newBasisFactorization ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::newBasisFactorization" ); delete _basisFactorization; @@ -1846,6 +1846,13 @@ void Tableau::registerCostFunctionManager( ICostFunctionManager *costFunctionMan _costFunctionManager = costFunctionManager; } +const double *Tableau::getColumnOfBasis( unsigned column ) const +{ + ASSERT( column < _m ); + unsigned variable = _basicIndexToVariable[column]; + return _A + ( variable * _m ); +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/engine/Tableau.h b/src/engine/Tableau.h index 906d1ab86..127651ebf 100644 --- a/src/engine/Tableau.h +++ b/src/engine/Tableau.h @@ -13,19 +13,19 @@ #ifndef __Tableau_h__ #define __Tableau_h__ +#include "IBasisFactorization.h" #include "ITableau.h" #include "MString.h" #include "Map.h" #include "Set.h" #include "Statistics.h" -class IBasisFactorization; class Equation; class ICostFunctionManager; class PiecewiseLinearCaseSplit; class TableauState; -class Tableau : public ITableau +class Tableau : public ITableau, public IBasisFactorization::BasisColumnOracle { public: Tableau(); @@ -385,6 +385,8 @@ class Tableau : public ITableau Equation *getBasisEquation( unsigned row ) const; double *getInverseBasisMatrix() const; + const double *getColumnOfBasis( unsigned column ) const; + private: typedef List VariableWatchers; Map _variableToWatchers; diff --git a/src/engine/TableauState.cpp b/src/engine/TableauState.cpp index fdeb77b14..a4c378939 100644 --- a/src/engine/TableauState.cpp +++ b/src/engine/TableauState.cpp @@ -91,7 +91,7 @@ TableauState::~TableauState() } } -void TableauState::setDimensions( unsigned m, unsigned n ) +void TableauState::setDimensions( unsigned m, unsigned n, const IBasisFactorization::BasisColumnOracle &oracle ) { _m = m; _n = n; @@ -132,7 +132,7 @@ void TableauState::setDimensions( unsigned m, unsigned n ) if ( !_variableToIndex ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "TableauState::variableToIndex" ); - _basisFactorization = BasisFactorizationFactory::createBasisFactorization( m ); + _basisFactorization = BasisFactorizationFactory::createBasisFactorization( m, oracle ); if ( !_basisFactorization ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "TableauState::basisFactorization" ); } diff --git a/src/engine/TableauState.h b/src/engine/TableauState.h index 114cd4af7..b117518e7 100644 --- a/src/engine/TableauState.h +++ b/src/engine/TableauState.h @@ -13,11 +13,10 @@ #ifndef __TableauState_h__ #define __TableauState_h__ +#include "IBasisFactorization.h" #include "ITableau.h" #include "Set.h" -class IBasisFactorization; - class TableauState { /* @@ -37,7 +36,7 @@ class TableauState TableauState(); ~TableauState(); - void setDimensions( unsigned m, unsigned n ); + void setDimensions( unsigned m, unsigned n, const IBasisFactorization::BasisColumnOracle &oracle ); /* The dimensions of matrix A From 6c7e6e36639f434b937b61781bbd8ba698dccb61 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 9 May 2018 11:02:01 +0300 Subject: [PATCH 05/14] cleanup in FT factorization --- .../ForrestTomlinFactorization.cpp | 146 +----------------- .../ForrestTomlinFactorization.h | 8 - src/basis_factorization/LUFactorization.cpp | 28 ++++ src/basis_factorization/LUFactorization.h | 5 + .../tests/MockColumnOracle.h | 30 +++- .../tests/Test_ForrestTomlinFactorization.h | 53 ++----- 6 files changed, 74 insertions(+), 196 deletions(-) diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp index 0bc950526..c0026168a 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.cpp +++ b/src/basis_factorization/ForrestTomlinFactorization.cpp @@ -25,7 +25,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisC , _Q( m ) , _invQ( m ) , _U( NULL ) - , _invLP( NULL ) , _explicitBasisAvailable( false ) , _workMatrix( NULL ) , _workVector( NULL ) @@ -44,11 +43,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisC throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, "ForrestTomlinFactorization::U" ); - _invLP = new double[m * m]; - if ( !_invLP ) - throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, - "ForrestTomlinFactorization::invLPx" ); - for ( unsigned i = 0; i < _m; ++i ) { _U[i] = new EtaMatrix( _m, i ); @@ -82,11 +76,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisC for ( unsigned row = 0; row < _m; ++row ) _B[row * _m + row] = 1.0; - // Also, invLP = I - std::fill_n( _invLP, _m * _m, 0.0 ); - for ( unsigned row = 0; row < _m; ++row ) - _invLP[row * _m + row] = 1.0; - // Us, Q and invQ are all initialized to the identity matrix anyway. } @@ -113,12 +102,6 @@ ForrestTomlinFactorization::~ForrestTomlinFactorization() _U = NULL; } - if ( _invLP ) - { - delete[] _invLP; - _invLP = NULL; - } - if ( _workMatrix ) { delete[] _workMatrix; @@ -311,8 +294,6 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub // If the number of A matrices is too great, condense them. if ( _A.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD ) refactorizeBasis(); - - // makeExplicitBasisAvailable(); } void ForrestTomlinFactorization::forwardTransformation( const double *y, double *x ) const @@ -536,8 +517,6 @@ void ForrestTomlinFactorization::storeFactorization( IBasisFactorization *other for ( const auto &lp : _LP ) otherFTFactorization->_LP.append( lp->duplicate() ); - memcpy( otherFTFactorization->_invLP, _invLP, sizeof(double) * _m * _m ); - List::iterator aIt; for ( aIt = otherFTFactorization->_A.begin(); aIt != otherFTFactorization->_A.end(); ++aIt ) delete *aIt; @@ -575,8 +554,6 @@ void ForrestTomlinFactorization::restoreFactorization( const IBasisFactorization for ( const auto &lp : otherFTFactorization->_LP ) _LP.append( lp->duplicate() ); - memcpy( _invLP, otherFTFactorization->_invLP, sizeof(double) * _m * _m ); - List::iterator aIt; for ( aIt = _A.begin(); aIt != _A.end(); ++aIt ) delete *aIt; @@ -694,56 +671,6 @@ void ForrestTomlinFactorization::initialLUFactorization() for ( unsigned j = 0; j <= i; ++j ) u->_column[j] = _workMatrix[j * _m + i]; } - - // Compute inv(P1)inv(L1) ... inv(Ps)inv(LS) for later. - - // Initialize to the identity matrix - std::fill_n( _invLP, _m * _m, 0 ); - for ( unsigned i = 0; i < _m; ++i ) - _invLP[i*_m + i] = 1.0; - - // Multiply by inverted Ps and Ls - for ( auto lp = _LP.rbegin(); lp != _LP.rend(); ++lp ) - { - if ( (*lp)->_pair ) - { - // The inverse of a general permutation matrix is its transpose. - // However, since these are permutations of just two rows, the - // transpose is equal to the original - so no need to invert. - // Multiplying on the right permutes columns. - double temp; - unsigned first = (*lp)->_pair->first; - unsigned second = (*lp)->_pair->second; - - for ( unsigned i = 0; i < _m; ++i ) - { - temp = _invLP[i * _m + first]; - _invLP[i * _m + first] = _invLP[i * _m + second]; - _invLP[i * _m + second] = temp; - } - } - else - { - // The inverse of an eta matrix is an eta matrix with the special - // column negated and divided by the diagonal entry. The only - // exception is the diagonal entry itself, which is just the - // inverse of the original diagonal entry. - EtaMatrix *eta = (*lp)->_eta; - unsigned columnIndex = eta->_columnIndex; - double etaDiagonalEntry = 1 / eta->_column[columnIndex]; - - for ( unsigned row = 0; row < _m; ++row ) - { - for ( unsigned col = columnIndex + 1; col < _m; ++col ) - { - _invLP[row * _m + columnIndex] -= - ( eta->_column[col] * _invLP[row * _m + col] ); - } - - _invLP[row * _m + columnIndex] *= etaDiagonalEntry; - } - } - } } bool ForrestTomlinFactorization::explicitBasisAvailable() const @@ -757,72 +684,6 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable() return; refactorizeBasis(); - return; - - /* - We know that the following equation holds: - - Ak...A1 * LsPs...L1P1 * B = Q * Um...U1 * invQ - - Or: - - B = inv( Ak...A1 * LsPs...L1P1 ) * Q * Um...U1 * invQ - = inv(P1)inv(L1) ... inv(Ps)inv(LS) * inv(A1)...inv(Ak) - * Q * Um...U1 * invQ - */ - - // Start with the already-computed inv(P1)inv(L1) ... inv(Ps)inv(LS) - memcpy( _B, _invLP, sizeof(double) * _m * _m ); - - // Multiply by the inverse As. An inverse of an almost identity matrix - // is just that matrix with its special value negated (unless that - // matrix is diagonal. - for ( const auto &a : _A ) - { - if ( a->_row == a->_column ) - { - for ( unsigned j = 0; j < _m; ++j ) - _B[j * _m + a->_column] *= ( 1 / a->_value ); - } - else - { - for ( unsigned j = 0; j < _m; ++j ) - _B[j * _m + a->_column] -= _B[j * _m + a->_row] * a->_value; - } - } - - // Permute columns according to Q - for ( unsigned i = 0; i < _m; ++i ) - { - unsigned columnToCopy = _invQ._ordering[i]; - for ( unsigned j = 0; j < _m; ++j ) - _workMatrix[j * _m + i] = _B[j * _m + columnToCopy]; - } - - // Multiply by Um...U1 - for ( int i = _m - 1; i >= 0; --i ) - { - for ( unsigned row = 0; row < _m; ++row ) - { - for ( unsigned j = 0; j < _U[i]->_columnIndex; ++j ) - { - _workMatrix[row * _m + _U[i]->_columnIndex] += - ( _U[i]->_column[j] * _workMatrix[row * _m + j] ); - } - } - } - - // Permute columns according to invQ - for ( unsigned i = 0; i < _m; ++i ) - { - unsigned columnToCopy = _Q._ordering[i]; - for ( unsigned j = 0; j < _m; ++j ) - _B[j * _m + i] = _workMatrix[j * _m + columnToCopy]; - } - - clearFactorization(); - initialLUFactorization(); - _explicitBasisAvailable = true; } const double *ForrestTomlinFactorization::getBasis() const @@ -995,12 +856,7 @@ void ForrestTomlinFactorization::dump() const printf( "\nDumping invQ:\n" ); _invQ.dump(); - printf( "*** Done dumping FT factorization ***\n\no" ); -} - -const double *ForrestTomlinFactorization::getInvLP() const -{ - return _invLP; + printf( "*** Done dumping FT factorization ***\n\n" ); } void ForrestTomlinFactorization::refactorizeBasis() diff --git a/src/basis_factorization/ForrestTomlinFactorization.h b/src/basis_factorization/ForrestTomlinFactorization.h index 1f22d21ef..f58fd3e87 100644 --- a/src/basis_factorization/ForrestTomlinFactorization.h +++ b/src/basis_factorization/ForrestTomlinFactorization.h @@ -102,7 +102,6 @@ class ForrestTomlinFactorization : public IBasisFactorization const EtaMatrix **getU() const; const List *getLP() const; const List *getA() const; - const double *getInvLP() const; void pushA( const AlmostIdentityMatrix &matrix ); void setQ( const PermutationMatrix &Q ); @@ -129,13 +128,6 @@ class ForrestTomlinFactorization : public IBasisFactorization // A1 is the first element of the list List_A; - /* - A helper matrix to facilitate computing the explicit basis. - This matrix stores: - inv(P1)inv(L1) ... inv(Ps)inv(LS) - */ - double *_invLP; - /* Indicates whether the explicit basis matrix is available. */ diff --git a/src/basis_factorization/LUFactorization.cpp b/src/basis_factorization/LUFactorization.cpp index d6c0e1153..f7955352d 100644 --- a/src/basis_factorization/LUFactorization.cpp +++ b/src/basis_factorization/LUFactorization.cpp @@ -557,6 +557,34 @@ void LUFactorization::makeExplicitBasisAvailable() condenseEtas(); } +void LUFactorization::dump() const +{ + printf( "*** Dumping LU factorization ***\n\n" ); + + printf( "\nDumping LPs:\n" ); + unsigned count = 0; + for ( const auto &lp : _LP ) + { + printf( "LP[%i]:\n", _LP.size() - 1 - count ); + ++count; + lp->dump(); + } + printf( "\n\n" ); + + printf( "Dumping Us:\n" ); + for ( unsigned col = 0; col < _m; ++col ) + { + printf( "U[%u]:\n", col ); + for ( unsigned i = 0; i < _m; ++i ) + { + printf( "\t%lf\n", _U[i * _m + col] ); + } + printf( "\n" ); + } + + printf( "*** Done dumping LU factorization ***\n\n" ); +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/basis_factorization/LUFactorization.h b/src/basis_factorization/LUFactorization.h index 6f506fa6a..9b31a6485 100644 --- a/src/basis_factorization/LUFactorization.h +++ b/src/basis_factorization/LUFactorization.h @@ -147,6 +147,11 @@ class LUFactorization : public IBasisFactorization const List getLP() const; const List getEtas() const; + /* + Debug + */ + void dump() const; + private: /* The Basis matrix. diff --git a/src/basis_factorization/tests/MockColumnOracle.h b/src/basis_factorization/tests/MockColumnOracle.h index fce3e9613..7fec144fa 100644 --- a/src/basis_factorization/tests/MockColumnOracle.h +++ b/src/basis_factorization/tests/MockColumnOracle.h @@ -18,7 +18,35 @@ class MockColumnOracle : public IBasisFactorization::BasisColumnOracle { public: - const double *_basis; + MockColumnOracle() + : _basis( NULL ) + { + } + + ~MockColumnOracle() + { + if ( _basis ) + { + delete[] _basis; + _basis = NULL; + } + } + + void storeBasis( unsigned m, const double *basis ) + { + _m = m; + _basis = new double[_m * _m]; + + for ( unsigned row = 0; row < _m; ++row ) + { + for ( unsigned column = 0; column < _m; ++column ) + { + _basis[column * _m + row] = basis[row * _m + column]; + } + } + } + + double *_basis; unsigned _m; const double *getColumnOfBasis( unsigned column ) const diff --git a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h index 839520969..d68242dd2 100644 --- a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h +++ b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h @@ -192,48 +192,6 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite ++lIt; TS_ASSERT_EQUALS( *((*lIt)->_eta), expectedL1 ); - /* - L4 = | 1 0 0 0 | - | 0 1 0 0 | - | 0 0 1 0 | - | 0 0 0 -0.5 | - - L3 = | 1 0 0 0 | - | 0 1 0 0 | - | 0 0 -1 0 | - | 0 0 1 1 | - - L2 = | 1 0 0 0 | - | 0 0.5 0 0 | - | 0 0 1 0 | - | 0 0 0 1 | - - L1 = | 1 0 0 0 | - | -1 1 0 0 | - | -1 0 1 0 | - | 1 0 0 1 | - - InvLP = inv( L4L3L2L1 ) - - = inv( | 1 0 0 0 | ) - | -0.5 0.5 0 0 | - | 1 0 -1 0 | - | 0 0 -0.5 -0.5 | - - = | 1 0 0 0 | - | 1 2 0 0 | - | 1 0 -1 0 | - | -1 0 1 -2 | - */ - double expectedInvLP[] = - { 1, 0, 0, 0, - 1, 2, 0, 0, - 1, 0, -1, 0, - -1, 0, 1, - 2}; - - for ( unsigned i = 0; i < 16; ++i ) - TS_ASSERT( FloatUtils::areEqual( expectedInvLP[i], ft->getInvLP()[i] ) ); - TS_ASSERT_THROWS_NOTHING( delete ft ); } @@ -938,6 +896,8 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite -1, -26, 3, -8, }; + oracle->storeBasis( 4, expectedB ); + TS_ASSERT( !ft->explicitBasisAvailable() ); TS_ASSERT_THROWS_NOTHING( ft->makeExplicitBasisAvailable() ); TS_ASSERT( ft->explicitBasisAvailable() ); @@ -999,6 +959,15 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite // | 1 20 -3 6 | // | -1 -26 3 -8 | + double expectedB[16] = { + 1, 14, -2, 4, + 1, 21, -1, 5, + 1, 20, -3, 6, + -1, -26, 3, -8, + }; + + oracle->storeBasis( 4, expectedB ); + // invB = | 4 -1/2 -13/4 -3/4 | // | -1/2 1/4 5/8 3/8 | // | 1 0 -2 -1 | From 450a30219eabcaa83551a8492a9a8946c15b0462 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 9 May 2018 11:58:45 +0300 Subject: [PATCH 06/14] grab a fresh basis from the oracle instead of condensing th etas --- src/basis_factorization/LUFactorization.cpp | 94 +++++++++---------- src/basis_factorization/LUFactorization.h | 9 +- .../tests/Test_CompareFactorizations.h | 11 +++ .../tests/Test_LUFactorization.h | 47 +++++----- 4 files changed, 79 insertions(+), 82 deletions(-) diff --git a/src/basis_factorization/LUFactorization.cpp b/src/basis_factorization/LUFactorization.cpp index f7955352d..108d7fafe 100644 --- a/src/basis_factorization/LUFactorization.cpp +++ b/src/basis_factorization/LUFactorization.cpp @@ -122,9 +122,8 @@ void LUFactorization::pushEtaMatrix( unsigned columnIndex, const double *column if ( ( _etas.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD ) && factorizationEnabled() ) { - log( "Number of etas exceeds threshold. Condensing and refactoring\n" ); - condenseEtas(); - factorizeMatrix( _B0 ); + log( "Number of etas exceeds threshold. Refactoring basis\n" ); + refactorizeBasis(); } } @@ -159,41 +158,10 @@ void LUFactorization::LMultiplyLeft( const EtaMatrix *L, double *X ) const void LUFactorization::setBasis( const double *B ) { memcpy( _B0, B, sizeof(double) * _m * _m ); + clearFactorization(); factorizeMatrix( _B0 ); } -void LUFactorization::condenseEtas() -{ - // Multiplication by an eta matrix on the right only changes one - // column of B0. The new column is a linear combination of the - // existing columns of B0, according to the eta column. We perform - // the computation in place. - for ( const auto &eta : _etas ) - { - unsigned col = eta->_columnIndex; - - // Iterate over the rows of B0 - for ( unsigned i = 0; i < _m; ++i ) - { - // Compute the linear combinations of the columns - double sum = 0.0; - for ( unsigned j = 0; j < _m; ++j ) - sum += _B0[i * _m + j] * eta->_column[j]; - - if ( FloatUtils::isZero( sum ) ) - sum = 0.0; - - _B0[i * _m + col] = sum; - } - - delete eta; - } - - _etas.clear(); - clearLPU(); - factorizeMatrix( _B0 ); -} - void LUFactorization::forwardTransformation( const double *y, double *x ) const { // If there's no LP factorization, it is implied that B0 = I. @@ -357,7 +325,7 @@ void LUFactorization::rowSwap( unsigned rowOne, unsigned rowTwo, double *matrix } } -void LUFactorization::clearLPU() +void LUFactorization::clearFactorization() { List::iterator element; for ( element = _LP.begin(); element != _LP.end(); ++element ) @@ -365,12 +333,13 @@ void LUFactorization::clearLPU() _LP.clear(); std::fill_n( _U, _m*_m, 0 ); + + _etas.clear(); } void LUFactorization::factorizeMatrix( double *matrix ) { - // Clear any previous factorization, initialize U - clearLPU(); + // Initialize U memcpy( _U, matrix, sizeof(double) * _m * _m ); for ( unsigned i = 0; i < _m; ++i ) @@ -448,12 +417,15 @@ void LUFactorization::storeFactorization( IBasisFactorization *other ) ASSERT( _m == otherLUFactorization->_m ); ASSERT( otherLUFactorization->_etas.size() == 0 ); - // In order to reduce space requirements, condense the etas before storing a factorization - condenseEtas(); - factorizeMatrix( _B0 ); + // In order to reduce space requirements (for the etas), compute a fresh refactorization + if ( !_etas.empty() ) + refactorizeBasis(); - // Now we simply store _B0 - otherLUFactorization->setBasis( _B0 ); + // Store the new basis and factorization + memcpy( otherLUFactorization->_B0, _B0, sizeof(double) * _m * _m ); + memcpy( otherLUFactorization->_U, _U, sizeof(double) * _m * _m ); + for ( const auto &it : _LP ) + otherLUFactorization->_LP.append( it->duplicate() ); } void LUFactorization::restoreFactorization( const IBasisFactorization *other ) @@ -464,14 +436,13 @@ void LUFactorization::restoreFactorization( const IBasisFactorization *other ) ASSERT( otherLUFactorization->_etas.size() == 0 ); // Clear any existing data - for ( const auto &it : _etas ) - delete it; + clearFactorization(); - _etas.clear(); - clearLPU(); - - // Store the new B0 and LU-factorize it - setBasis( otherLUFactorization->_B0 ); + // Store the new basis and factorization + memcpy( _B0, otherLUFactorization->_B0, sizeof(double) * _m * _m ); + memcpy( _U, otherLUFactorization->_U, sizeof(double) * _m * _m ); + for ( const auto &it : otherLUFactorization->_LP ) + _LP.append( it->duplicate() ); } void LUFactorization::invertBasis( double *result ) @@ -554,7 +525,7 @@ bool LUFactorization::explicitBasisAvailable() const void LUFactorization::makeExplicitBasisAvailable() { - condenseEtas(); + refactorizeBasis(); } void LUFactorization::dump() const @@ -582,9 +553,30 @@ void LUFactorization::dump() const printf( "\n" ); } + printf( "Dumping etas:\n" ); + for ( const auto &eta : _etas ) + { + eta->dump(); + printf( "\n" ); + } + printf( "*** Done dumping LU factorization ***\n\n" ); } +void LUFactorization::refactorizeBasis() +{ + clearFactorization(); + + for ( unsigned column = 0; column < _m; ++column ) + { + const double *basisColumn = _basisColumnOracle->getColumnOfBasis( column ); + for ( unsigned row = 0; row < _m; ++row ) + _B0[row * _m + column] = basisColumn[row]; + } + + factorizeMatrix( _B0 ); +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/basis_factorization/LUFactorization.h b/src/basis_factorization/LUFactorization.h index 9b31a6485..53f8e9df0 100644 --- a/src/basis_factorization/LUFactorization.h +++ b/src/basis_factorization/LUFactorization.h @@ -107,12 +107,6 @@ class LUFactorization : public IBasisFactorization */ void rowSwap( unsigned rowOne, unsigned rowTwo, double *matrix ); - /* - Compute B0 * E1 ... *En for all stored eta matrices, and place - the result in B0. - */ - void condenseEtas(); - /* Return true iff the basis matrix B0 is explicitly available. */ @@ -122,6 +116,7 @@ class LUFactorization : public IBasisFactorization Make the basis explicitly available */ void makeExplicitBasisAvailable(); + void refactorizeBasis(); /* Get the explicit basis matrix @@ -185,7 +180,7 @@ class LUFactorization : public IBasisFactorization /* Clear a previous factorization. */ - void clearLPU(); + void clearFactorization(); /* Helper functions for backward- and forward-transformations. diff --git a/src/basis_factorization/tests/Test_CompareFactorizations.h b/src/basis_factorization/tests/Test_CompareFactorizations.h index 0b4c0ab63..34fa2dd0e 100644 --- a/src/basis_factorization/tests/Test_CompareFactorizations.h +++ b/src/basis_factorization/tests/Test_CompareFactorizations.h @@ -91,7 +91,18 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite for ( unsigned i = 0; i < 4; ++i ) TS_ASSERT( FloatUtils::areEqual( x1[i], x2[i] ) ); + + double basisAtThisPoint[] = { + 2, -20, 64, -4, + 0, 2, 96, 0, + -3, 24, -45, 1, + 0, 6, 14, 2, + }; + + oracle->storeBasis( 4, basisAtThisPoint ); + ft->makeExplicitBasisAvailable(); + lu->makeExplicitBasisAvailable(); TS_ASSERT_THROWS_NOTHING( ft->forwardTransformation( y, x1 ) ); TS_ASSERT_THROWS_NOTHING( lu->forwardTransformation( y, x2 ) ); diff --git a/src/basis_factorization/tests/Test_LUFactorization.h b/src/basis_factorization/tests/Test_LUFactorization.h index d63bc9e82..6a2cab76f 100644 --- a/src/basis_factorization/tests/Test_LUFactorization.h +++ b/src/basis_factorization/tests/Test_LUFactorization.h @@ -146,7 +146,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite double a[] = { 2., -1., 4. }; double d[] = { 0., 0., 0. }; - double expected[] = { 42, 116, -131 }; + double expected[] = { -20, 27, -8 }; basis.forwardTransformation( a, d ); @@ -264,12 +264,13 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite double y[] = {19., 12., 17.}; double x[] = {0., 0., 0.}; - double expected[] = {-6, 9, -4}; - // | 1 2 4 | | 1 1 | | 2 | | 1 0.5 | - // x * | 4 5 7 | * | 1 | * | 1 1 | * | 1 0.5 | = | 19 12 17 | - // | 7 8 9 | | 3 1 | | 1 1 | | 0.5 | + double expected[] = { -104.0/3, 140.0/3, -19 }; + + // | 1 2 4 | + // x * | 4 5 7 | = | 19 12 17 | + // | 7 8 9 | // - // --> x = [ -6 9 -4 ] + // --> x = [ -104/3, 140/3, -19 ] basis.backwardTransformation( y, x ); for ( unsigned i = 0; i < 3; ++i ) @@ -287,6 +288,14 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( basis.forwardTransformation( a1, d1 ) ); basis.pushEtaMatrix( 1, a1 ); + // Save the expected basis after this push + double currentBasis[] = { + 1, 1, 0, + 0, 1, 0, + 0, 3, 1 + }; + oracle->storeBasis( 3, currentBasis ); + // Do a computation using both basis, see that we get the same result. double a2[] = { 3, 1, 4 }; @@ -554,11 +563,8 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite TS_ASSERT( FloatUtils::areEqual( eta->_column[i], expectedCol3[i] ) ); } - void xtest_refactor() + void test_refactor() { - // TODO: this test fails when the REFACTORIZATION_THRESHOLD is too great (> 10 or so). - // Disabling for now. - LUFactorization basis( 3, *oracle ); LUFactorization basis2( 3, *oracle ); basis.toggleFactorization( false ); @@ -578,6 +584,13 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite etaPool[i] = (float)(rand()) / (float)(RAND_MAX); } + double dummyBasis[] = { + 1, 2, 0, + 0, 5, 0, + 0, 8, 9, + }; + oracle->storeBasis( 3, dummyBasis ); + // Generate random etas for ( unsigned i = 0; i < etaCount; ++i ) { @@ -590,20 +603,6 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite // Check if etas have disappeared TS_ASSERT_EQUALS( basis2.getEtas().size(), etaCount - GlobalConfiguration::REFACTORIZATION_THRESHOLD - 1 ); - double a[] = {2., -1., 4.}; - double x1[] = {0., 0., 0.}; - double y1[] = {0., 0., 0.}; - double x2[] = {0., 0., 0.}; - double y2[] = {0., 0., 0.}; - basis2.forwardTransformation( a, x1 ); - basis.forwardTransformation( a, y1 ); - basis2.backwardTransformation( a, x2 ); - basis.backwardTransformation( a, y2 ); - for ( unsigned i = 0; i < 3; ++i ) - { - TS_ASSERT( FloatUtils::areEqual( x1[i], y1[i], 0.001 ) ); - TS_ASSERT( FloatUtils::areEqual( x2[i], y2[i], 0.001 ) ); - } } void test_matrix_multiply() From 23ba33a55807551de91415b4604655e80e9fea96 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Thu, 10 May 2018 15:00:09 +0300 Subject: [PATCH 07/14] bug fix in test --- src/input_parsers/acas_example/main.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/input_parsers/acas_example/main.cpp b/src/input_parsers/acas_example/main.cpp index 0186e54b5..ba99d6021 100644 --- a/src/input_parsers/acas_example/main.cpp +++ b/src/input_parsers/acas_example/main.cpp @@ -60,9 +60,7 @@ int main() // Feed the query to the engine Engine engine; - engine.processInputQuery( inputQuery ); - - if ( !engine.solve() ) + if ( !engine.processInputQuery( inputQuery ) || !engine.solve() ) { printf( "\n\nQuery is unsat\n" ); return 0; From 701a3c2bb23bd01ae2de4f31734d371b79d83984 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Mon, 21 May 2018 17:04:08 +0300 Subject: [PATCH 08/14] started separating out the gaussian elimination functionality from the basis factorization classes --- .../GaussianEliminator.cpp | 167 +++++++++++++ src/basis_factorization/GaussianEliminator.h | 73 ++++++ src/basis_factorization/Sources.mk | 1 + .../tests/Test_GaussianEliminator.h | 222 ++++++++++++++++++ 4 files changed, 463 insertions(+) create mode 100644 src/basis_factorization/GaussianEliminator.cpp create mode 100644 src/basis_factorization/GaussianEliminator.h create mode 100644 src/basis_factorization/tests/Test_GaussianEliminator.h diff --git a/src/basis_factorization/GaussianEliminator.cpp b/src/basis_factorization/GaussianEliminator.cpp new file mode 100644 index 000000000..613053ac2 --- /dev/null +++ b/src/basis_factorization/GaussianEliminator.cpp @@ -0,0 +1,167 @@ +/******************** */ +/*! \file GaussianEliminator.cpp + ** \verbatim + ** Top contributors (to current version): + ** Guy Katz + ** This file is part of the Marabou project. + ** Copyright (c) 2016-2017 by the authors listed in the file AUTHORS + ** in the top-level source directory) and their institutional affiliations. + ** All rights reserved. See the file COPYING in the top-level source + ** directory for licensing information.\endverbatim + **/ + +#include "BasisFactorizationError.h" +#include "EtaMatrix.h" +#include "FloatUtils.h" +#include "GaussianEliminator.h" +#include "MalformedBasisException.h" + +#include + +GaussianEliminator::GaussianEliminator( const double *A, unsigned m ) + : _A( A ) + , _m( m ) + , _pivotColumn( NULL ) + , _LCol( NULL ) +{ + _pivotColumn = new double[_m]; + if ( !_pivotColumn ) + throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, + "GaussianEliminator::pivotColumn" ); + + _LCol = new double[_m]; + if ( !_LCol ) + throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED, "LUFactorization::LCol" ); +} + +GaussianEliminator::~GaussianEliminator() +{ + if ( _pivotColumn ) + { + delete[] _pivotColumn; + _pivotColumn = NULL; + } + + if ( _LCol ) + { + delete[] _LCol; + _LCol = NULL; + } +} + +unsigned GaussianEliminator::choosePivotElement( unsigned columnIndex, FactorizationStrategy factorizationStrategy ) +{ + if ( factorizationStrategy != PARTIAL_PIVOTING ) + { + printf( "Error! Only strategy currently supported is partial pivoting!\n" ); + exit( 1 ); + } + + double largestElement = FloatUtils::abs( _pivotColumn[columnIndex] ); + unsigned bestRowIndex = columnIndex; + + for ( unsigned i = columnIndex + 1; i < _m; ++i ) + { + double contender = FloatUtils::abs( _pivotColumn[i] ); + if ( FloatUtils::gt( contender, largestElement ) ) + { + largestElement = contender; + bestRowIndex = i; + } + } + + // No non-zero pivot has been found, matrix cannot be factorized + if ( FloatUtils::isZero( largestElement ) ) + throw MalformedBasisException(); + + return bestRowIndex; +} + +void GaussianEliminator::factorize( List *lp, + double *U, + unsigned *rowHeaders, + FactorizationStrategy factorizationStrategy ) +{ + // Initialize U to the stored A + memcpy( U, _A, sizeof(double) * _m * _m ); + + // Reset the row headers. The i'th row is stored as the rowHeaders[i]'th row in memory + for ( unsigned i = 0; i < _m; ++i ) + rowHeaders[i] = i; + + // Iterate over the columns of U + for ( unsigned pivotColumnIndex = 0; pivotColumnIndex < _m; ++pivotColumnIndex ) + { + // Extract the current column + for ( unsigned i = pivotColumnIndex; i < _m; ++i ) + _pivotColumn[i] = U[rowHeaders[i] * _m + pivotColumnIndex]; + + // Select the new pivot element according to the factorization strategy + unsigned bestRowIndex = choosePivotElement( pivotColumnIndex, factorizationStrategy ); + + // Swap the pivot row to the top of the column if needed + if ( bestRowIndex != pivotColumnIndex ) + { + // Swap the rows through the headers + unsigned tempIndex = rowHeaders[pivotColumnIndex]; + rowHeaders[pivotColumnIndex] = rowHeaders[bestRowIndex]; + rowHeaders[bestRowIndex] = tempIndex; + + // Also swap in the work vector + double tempVal = _pivotColumn[pivotColumnIndex]; + _pivotColumn[pivotColumnIndex] = _pivotColumn[bestRowIndex]; + _pivotColumn[bestRowIndex] = tempVal; + + std::pair *P = new std::pair( pivotColumnIndex, bestRowIndex ); + lp->appendHead( new LPElement( NULL, P ) ); + } + + // The matrix now has a non-zero value at entry (pivotColumnIndex, pivotColumnIndex), + // so we can perform Gaussian elimination for the subsequent rows + + // TODO: don't need the full LCol, can switch to a dedicated lower-triangular matrix datatype + bool eliminationRequired = !FloatUtils::areEqual( _pivotColumn[pivotColumnIndex], 1.0 ); + std::fill_n( _LCol, _m, 0 ); + _LCol[pivotColumnIndex] = 1 / _pivotColumn[pivotColumnIndex]; + for ( unsigned j = pivotColumnIndex + 1; j < _m; ++j ) + { + if ( !FloatUtils::isZero( _pivotColumn[j] ) ) + { + eliminationRequired = true; + _LCol[j] = -_pivotColumn[j] * _LCol[pivotColumnIndex]; + } + } + + if ( eliminationRequired ) + { + // Store the resulting lower-triangular eta matrix + EtaMatrix *L = new EtaMatrix( _m, pivotColumnIndex, _LCol ); + lp->appendHead( new LPElement( L, NULL ) ); + + // Perform the actual elimination step on U. First, perform in-place multiplication + // for all rows below the pivot row + for ( unsigned row = pivotColumnIndex + 1; row < _m; ++row ) + { + U[rowHeaders[row] * _m + pivotColumnIndex] = 0.0; + for ( unsigned col = pivotColumnIndex + 1; col < _m; ++col ) + U[rowHeaders[row] * _m + col] += L->_column[row] * U[rowHeaders[pivotColumnIndex] * _m + col]; + } + + // Finally, perform the multiplication for the pivot row + // itself. We change this row last because it is required for all + // previous multiplication operations. + for ( unsigned i = pivotColumnIndex + 1; i < _m; ++i ) + U[rowHeaders[pivotColumnIndex] * _m + i] *= L->_column[pivotColumnIndex]; + + U[rowHeaders[pivotColumnIndex] * _m + pivotColumnIndex] = 1.0; + } + } +} + +// +// Local Variables: +// compile-command: "make -C ../.. " +// tags-file-name: "../../TAGS" +// c-basic-offset: 4 +// End: +// diff --git a/src/basis_factorization/GaussianEliminator.h b/src/basis_factorization/GaussianEliminator.h new file mode 100644 index 000000000..63444f0c3 --- /dev/null +++ b/src/basis_factorization/GaussianEliminator.h @@ -0,0 +1,73 @@ +/********************* */ +/*! \file GaussianEliminator.h + ** \verbatim + ** Top contributors (to current version): + ** Guy Katz + ** This file is part of the Marabou project. + ** Copyright (c) 2016-2017 by the authors listed in the file AUTHORS + ** in the top-level source directory) and their institutional affiliations. + ** All rights reserved. See the file COPYING in the top-level source + ** directory for licensing information.\endverbatim + **/ + +#ifndef __GaussianEliminator_h__ +#define __GaussianEliminator_h__ + +#include "LPElement.h" +#include "List.h" + +class GaussianEliminator +{ +public: + enum FactorizationStrategy { + PARTIAL_PIVOTING = 0, + MARKOWITZ, + }; + + GaussianEliminator( const double *A, unsigned m ); + ~GaussianEliminator(); + + /* + The classes main method: factorize matrix A, given in row-major format, + into a matrix U and a sequence of L and P matrices, such that: + + - U is upper triangular with its diagonal set to 1s + - The Ls are lower triangular + - The Ps are permutation matrices + - The rowHeaders array indicates the orders of the rows of U, where the i'th + row is stored in the rowHeaders[i] location in memory + */ + void factorize( List *lp, + double *U, + unsigned *rowHeaders, + FactorizationStrategy factorizationStrategy = PARTIAL_PIVOTING ); + +private: + /* + The (square) matrix being factorized and its dimension + */ + const double *_A; + unsigned _m; + + /* + Given a column with element in indices [columnIndex, .., _m], choose the next pivot according to + the factorization strategy. + */ + unsigned choosePivotElement( unsigned columnIndex, FactorizationStrategy factorizationStrategy ); + + /* + Work memory + */ + double *_pivotColumn; + double *_LCol; +}; + +#endif // __GaussianEliminator_h__ + +// +// Local Variables: +// compile-command: "make -C ../.. " +// tags-file-name: "../../TAGS" +// c-basic-offset: 4 +// End: +// diff --git a/src/basis_factorization/Sources.mk b/src/basis_factorization/Sources.mk index 32d60e40e..edf4d7ec7 100644 --- a/src/basis_factorization/Sources.mk +++ b/src/basis_factorization/Sources.mk @@ -2,6 +2,7 @@ SOURCES += \ BasisFactorizationFactory.cpp \ EtaMatrix.cpp \ ForrestTomlinFactorization.cpp \ + GaussianEliminator.cpp \ LPElement.cpp \ LUFactorization.cpp \ PermutationMatrix.cpp \ diff --git a/src/basis_factorization/tests/Test_GaussianEliminator.h b/src/basis_factorization/tests/Test_GaussianEliminator.h new file mode 100644 index 000000000..6ca99dfd9 --- /dev/null +++ b/src/basis_factorization/tests/Test_GaussianEliminator.h @@ -0,0 +1,222 @@ +/********************* */ +/*! \file Test_GaussianEliminator.h +** \verbatim +** Top contributors (to current version): +** Derek Huang +** Guy Katz +** This file is part of the Marabou project. +** Copyright (c) 2016-2017 by the authors listed in the file AUTHORS +** in the top-level source directory) and their institutional affiliations. +** All rights reserved. See the file COPYING in the top-level source +** directory for licensing information.\endverbatim +**/ + +#include + +#include "EtaMatrix.h" +#include "FloatUtils.h" +#include "GaussianEliminator.h" +#include + +class MockForGaussianEliminator +{ +public: +}; + +class GaussianEliminatorTestSuite : public CxxTest::TestSuite +{ +public: + MockForGaussianEliminator *mock; + + void setUp() + { + TS_ASSERT( mock = new MockForGaussianEliminator ); + } + + void tearDown() + { + TS_ASSERT_THROWS_NOTHING( delete mock ); + } + + void test_identity() + { + double A[] = + { + 1, 0, 0, + 0, 1, 0, + 0, 0, 1 + }; + + GaussianEliminator eliminator( A, 3 ); + + List lp; + double U[9]; + unsigned rowHeaders[3]; + + TS_ASSERT_THROWS_NOTHING( eliminator.factorize( &lp, U, rowHeaders ) ); + + double expectedU[] = + { + 1, 0, 0, + 0, 1, 0, + 0, 0, 1 + }; + + TS_ASSERT_SAME_DATA( U, expectedU, sizeof(U) ); + TS_ASSERT( lp.empty() ); + } + + void applyLpOnA( List &lp, double *A, unsigned m ) + { + for ( auto op = lp.rbegin(); op != lp.rend(); ++op ) + { + LPElement &element( *(*op) ); + if ( element._pair ) + { + double *temp = new double[m]; + memcpy( temp, A + ( m * element._pair->first ), sizeof(double) * m ); + memcpy( A + ( m * element._pair->first ), A + ( m * element._pair->second ), sizeof(double) * m ); + memcpy( A + ( m * element._pair->second ), temp, sizeof(double) * m ); + delete[] temp; + } + else + { + EtaMatrix &eta( *( element._eta ) ); + for ( unsigned i = 0; i < m; ++i ) + { + if ( i == eta._columnIndex ) + continue; + + for ( unsigned j = 0; j < m; ++j ) + { + A[i*m + j] += eta._column[i] * A[eta._columnIndex*m + j]; + } + } + + for ( unsigned j = 0; j < m; ++j ) + A[eta._columnIndex*m + j] *= eta._column[eta._columnIndex]; + } + } + } + + void dumpMatrix( double *A, unsigned m ) + { + for ( unsigned i = 0; i < m; ++i ) + { + for ( unsigned j = 0; j < m; ++j ) + { + printf( "%.2lf ", A[i*m + j] ); + } + printf( "\n" ); + } + + printf( "\n\n" ); + } + + void dumpMatrix( double *A, unsigned m, unsigned *rowHeaders ) + { + for ( unsigned i = 0; i < m; ++i ) + { + for ( unsigned j = 0; j < m; ++j ) + { + printf( "%.2lf ", A[rowHeaders[i]*m + j] ); + } + printf( "\n" ); + } + + printf( "\n\n" ); + } + + void test_sanity() + { + { + double A[] = + { + 1, 1, 1, + 0, 1, 0, + 0, 1, 1 + }; + + unsigned m = 3; + GaussianEliminator eliminator( A, m ); + List lp; + double U[m * m]; + unsigned rowHeaders[m]; + + TS_ASSERT_THROWS_NOTHING( eliminator.factorize( &lp, U, rowHeaders ) ); + + applyLpOnA( lp, A, m ); + + for ( unsigned i = 0; i < m; ++i ) + { + for ( unsigned j = 0; j < m; ++j ) + { + TS_ASSERT( FloatUtils::areEqual( A[i*m + j], U[rowHeaders[i]*m + j] ) ); + } + } + } + + { + double A[] = + { + 1, 2, 3, + 0, 4, 5, + 2, 1, 3 + }; + + unsigned m = 3; + GaussianEliminator eliminator( A, m ); + List lp; + double U[m * m]; + unsigned rowHeaders[m]; + + TS_ASSERT_THROWS_NOTHING( eliminator.factorize( &lp, U, rowHeaders ) ); + + applyLpOnA( lp, A, m ); + + for ( unsigned i = 0; i < m; ++i ) + { + for ( unsigned j = 0; j < m; ++j ) + { + TS_ASSERT( FloatUtils::areEqual( A[i*m + j], U[rowHeaders[i]*m + j] ) ); + } + } + } + + { + double A[] = + { + 1, 2, 3, -3, + -2, 4, 5, 1, + 2, 1, 3, 0, + 0, 0, 1, 0, + }; + + unsigned m = 4; + GaussianEliminator eliminator( A, m ); + List lp; + double U[m * m]; + unsigned rowHeaders[m]; + + TS_ASSERT_THROWS_NOTHING( eliminator.factorize( &lp, U, rowHeaders ) ); + + applyLpOnA( lp, A, m ); + + for ( unsigned i = 0; i < m; ++i ) + { + for ( unsigned j = 0; j < m; ++j ) + { + TS_ASSERT( FloatUtils::areEqual( A[i*m + j], U[rowHeaders[i]*m + j] ) ); + } + } + } + } +}; + +// +// Local Variables: +// compile-command: "make -C ../../.. " +// tags-file-name: "../../../TAGS" +// c-basic-offset: 4 +// End: +// From ab0a98c011ad82691b34915b466fd735f1c7258a Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 22 May 2018 09:59:54 +0300 Subject: [PATCH 09/14] error when adding an equation --- src/engine/Engine.cpp | 10 +++++----- src/engine/ReluplexError.h | 1 + src/engine/Tableau.cpp | 11 ++++++++++- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/engine/Engine.cpp b/src/engine/Engine.cpp index ad9c36e70..05d8d7b5e 100644 --- a/src/engine/Engine.cpp +++ b/src/engine/Engine.cpp @@ -832,16 +832,16 @@ void Engine::applyValidConstraintCaseSplit( PiecewiseLinearConstraint *constrain { if ( constraint->isActive() && constraint->phaseFixed() ) { + String constraintString; + constraint->dump( constraintString ); + log( Stringf( "A constraint has become valid. Dumping constraint: %s", + constraintString.ascii() ) ); + constraint->setActiveConstraint( false ); PiecewiseLinearCaseSplit validSplit = constraint->getValidCaseSplit(); _smtCore.recordImpliedValidSplit( validSplit ); applySplit( validSplit ); ++_numPlConstraintsDisabledByValidSplits; - - String constraintString; - constraint->dump( constraintString ); - log( Stringf( "A constraint has become valid. Dumping constraint: %s", - constraintString.ascii() ) ); } } diff --git a/src/engine/ReluplexError.h b/src/engine/ReluplexError.h index 07bce341a..081347103 100644 --- a/src/engine/ReluplexError.h +++ b/src/engine/ReluplexError.h @@ -33,6 +33,7 @@ class ReluplexError : public Error PREPROCESSOR_CANT_FIND_NEW_AUXILIARY_VAR = 11, RESTORATION_FAILED_TO_RESTORE_PRECISION = 12, CANNOT_RESTORE_TABLEAU = 13, + FAILURE_TO_ADD_NEW_EQUATION = 14, DEBUGGING_ERROR = 999, }; diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 88e3e4e2b..464f0e129 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -17,6 +17,7 @@ #include "FloatUtils.h" #include "ICostFunctionManager.h" #include "MStringf.h" +#include "MalformedBasisException.h" #include "PiecewiseLinearCaseSplit.h" #include "ReluplexError.h" #include "Tableau.h" @@ -1259,7 +1260,15 @@ void Tableau::addEquation( const Equation &equation ) _basicAssignment[_m - 1] = 0.0; // Refactorize the basis - _basisFactorization->refactorizeBasis(); + try + { + _basisFactorization->refactorizeBasis(); + } + catch ( MalformedBasisException & ) + { + log( "addEquation failed - could not refactorize basis" ); + throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); + } // Notify about the new variable's assignment and compute its status notifyVariableValue( _basicIndexToVariable[_m - 1], _basicAssignment[_m - 1] ); From 5a04155d5d98676c8f69160257896360ee8178d1 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 23 May 2018 09:08:02 +0300 Subject: [PATCH 10/14] handle the initial assignment of basic variables as part of Tableau::initializeTableau(0 --- src/engine/Engine.cpp | 38 +++++-------- src/engine/ITableau.h | 2 +- src/engine/Tableau.cpp | 95 ++++++++++++++++++++++----------- src/engine/Tableau.h | 2 +- src/engine/tests/MockTableau.h | 5 +- src/engine/tests/Test_Tableau.h | 79 +++++++++------------------ 6 files changed, 111 insertions(+), 110 deletions(-) diff --git a/src/engine/Engine.cpp b/src/engine/Engine.cpp index 05d8d7b5e..f195294c8 100644 --- a/src/engine/Engine.cpp +++ b/src/engine/Engine.cpp @@ -562,28 +562,6 @@ bool Engine::processInputQuery( InputQuery &inputQuery, bool preprocess ) ++equationIndex; } - // Placeholder: better constraint matrix analysis as part - // of the preprocessing phase. - - AutoConstraintMatrixAnalyzer analyzer; - analyzer->analyze( _tableau->getA(), _tableau->getM(), _tableau->getN() ); - - if ( analyzer->getRank() != _tableau->getM() ) - { - printf( "Warning!! Contraint matrix rank is %u (out of %u)\n", - analyzer->getRank(), _tableau->getM() ); - } - - List independentColumns = analyzer->getIndependentColumns(); - - unsigned assigned = 0; - for( unsigned basicVar : independentColumns ) - { - _tableau->markAsBasic( basicVar ); - _tableau->assignIndexToBasicVariable( basicVar, assigned ); - assigned++; - } - for ( unsigned i = 0; i < n; ++i ) { _tableau->setLowerBound( i, _preprocessedQuery.getLowerBound( i ) ); @@ -591,7 +569,6 @@ bool Engine::processInputQuery( InputQuery &inputQuery, bool preprocess ) } _tableau->registerToWatchAllVariables( _rowBoundTightener ); - _rowBoundTightener->initialize( _tableau ); _plConstraints = _preprocessedQuery.getPiecewiseLinearConstraints(); @@ -601,7 +578,20 @@ bool Engine::processInputQuery( InputQuery &inputQuery, bool preprocess ) constraint->setStatistics( &_statistics ); } - _tableau->initializeTableau(); + // Placeholder: better constraint matrix analysis as part + // of the preprocessing phase. + + AutoConstraintMatrixAnalyzer analyzer; + analyzer->analyze( _tableau->getA(), _tableau->getM(), _tableau->getN() ); + + if ( analyzer->getRank() != _tableau->getM() ) + { + printf( "Warning!! Contraint matrix rank is %u (out of %u)\n", + analyzer->getRank(), _tableau->getM() ); + } + List independentColumns = analyzer->getIndependentColumns(); + _tableau->initializeTableau( independentColumns ); + _costFunctionManager->initialize(); _tableau->registerCostFunctionManager( _costFunctionManager ); _activeEntryStrategy->initialize( _tableau ); diff --git a/src/engine/ITableau.h b/src/engine/ITableau.h index 9dd221878..127a1de38 100644 --- a/src/engine/ITableau.h +++ b/src/engine/ITableau.h @@ -82,7 +82,7 @@ class ITableau virtual void setRightHandSide( const double *b ) = 0; virtual void setRightHandSide( unsigned index, double value ) = 0; virtual void markAsBasic( unsigned variable ) = 0; - virtual void initializeTableau() = 0; + virtual void initializeTableau( const List &initialBasicVariables ) = 0; virtual double getValue( unsigned variable ) = 0; virtual bool allBoundsValid() const = 0; virtual double getLowerBound( unsigned variable ) const = 0; diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 464f0e129..4180d5f94 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -244,11 +244,21 @@ void Tableau::assignIndexToBasicVariable( unsigned variable, unsigned index ) _variableToIndex[variable] = index; } -void Tableau::initializeTableau() +void Tableau::initializeTableau( const List &initialBasicVariables ) { - unsigned nonBasicIndex = 0; + _basicVariables.clear(); - // Assign variable indices + // Assign the basic indices + unsigned basicIndex = 0; + for( unsigned basicVar : initialBasicVariables ) + { + markAsBasic( basicVar ); + assignIndexToBasicVariable( basicVar, basicIndex ); + ++basicIndex; + } + + // Assign the non-basic indices + unsigned nonBasicIndex = 0; for ( unsigned i = 0; i < _n; ++i ) { if ( !_basicVariables.exists( i ) ) @@ -1233,46 +1243,71 @@ void Tableau::addEquation( const Equation &equation ) // Invalidate the cost function, so that it is recomputed in the next iteration. _costFunctionManager->invalidateCostFunction(); - // Mark the auxiliary variable as basic, add to indices - _basicVariables.insert( equation._auxVariable ); - _basicIndexToVariable[_m - 1] = equation._auxVariable; - _variableToIndex[equation._auxVariable] = _m - 1; - - // Populate the new row of A, compute the assignment for the new basic variable + // Populate the new row of A and b _b[_m - 1] = equation._scalar; - _basicAssignment[_m - 1] = equation._scalar; - double auxCoefficient = 0.0; for ( const auto &addend : equation._addends ) - { setEntryValue( _m - 1, addend._variable, addend._coefficient ); - if ( addend._variable == equation._auxVariable ) - auxCoefficient = addend._coefficient; - else - _basicAssignment[_m - 1] -= addend._coefficient * getValue( addend._variable ); - } - - ASSERT( !FloatUtils::isZero( auxCoefficient ) ); - _basicAssignment[_m - 1] = _basicAssignment[_m - 1] / auxCoefficient; - ASSERT( FloatUtils::wellFormed( _basicAssignment[_m - 1] ) ); - - if ( FloatUtils::isZero( _basicAssignment[_m - 1] ) ) - _basicAssignment[_m - 1] = 0.0; + /* + Attempt to make the auxiliary variable the new basic variable. + This usually works. + If it doesn't, compute a new set of basic variables (which is more + computationally expensive) + */ + _basicIndexToVariable[_m - 1] = equation._auxVariable; + _variableToIndex[equation._auxVariable] = _m - 1; + _basicVariables.insert( equation._auxVariable ); - // Refactorize the basis + // Attempt to refactorize the basis + bool factorizationSuccessful = true; try { _basisFactorization->refactorizeBasis(); } catch ( MalformedBasisException & ) { - log( "addEquation failed - could not refactorize basis" ); - throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); + factorizationSuccessful = false; + } + + if ( factorizationSuccessful ) + { + _basicAssignment[_m - 1] = equation._scalar; + double auxCoefficient = 0.0; + for ( const auto &addend : equation._addends ) + { + if ( addend._variable == equation._auxVariable ) + auxCoefficient = addend._coefficient; + else + _basicAssignment[_m - 1] -= addend._coefficient * getValue( addend._variable ); + } + + ASSERT( !FloatUtils::isZero( auxCoefficient ) ); + _basicAssignment[_m - 1] = _basicAssignment[_m - 1] / auxCoefficient; + ASSERT( FloatUtils::wellFormed( _basicAssignment[_m - 1] ) ); + + if ( FloatUtils::isZero( _basicAssignment[_m - 1] ) ) + _basicAssignment[_m - 1] = 0.0; + + // Notify about the new variable's assignment and compute its status + notifyVariableValue( _basicIndexToVariable[_m - 1], _basicAssignment[_m - 1] ); + computeBasicStatus( _m - 1 ); + return; + } + else + { + // ConstraintMatrixAnalyzer analyzer; + // analyzer.analyze( _A, _m, _n ); + // List independentColumns = analyzer->getIndependentColumns(); + + // Need to assign basic variables and initialize a whole bunch of stuff... } - // Notify about the new variable's assignment and compute its status - notifyVariableValue( _basicIndexToVariable[_m - 1], _basicAssignment[_m - 1] ); - computeBasicStatus( _m - 1 ); + // if ( !factorizationSuccessful ) + // { + // printf( "Out of candidates, terminating\n" ); + // log( "addEquation failed - could not refactorize basis" ); + // throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); + // } } void Tableau::addRow() diff --git a/src/engine/Tableau.h b/src/engine/Tableau.h index 127651ebf..ef45d13fd 100644 --- a/src/engine/Tableau.h +++ b/src/engine/Tableau.h @@ -85,7 +85,7 @@ class Tableau : public ITableau, public IBasisFactorization::BasisColumnOracle to lower bounds and computes the assignment. Assign the initial basic indices according to the equations. */ - void initializeTableau(); + void initializeTableau( const List &initialBasicVariables ); void assignIndexToBasicVariable( unsigned variable, unsigned index ); /* diff --git a/src/engine/tests/MockTableau.h b/src/engine/tests/MockTableau.h index c0fdfbbac..fd8a05f0c 100644 --- a/src/engine/tests/MockTableau.h +++ b/src/engine/tests/MockTableau.h @@ -152,10 +152,13 @@ class MockTableau : public ITableau } bool initializeTableauCalled; - void initializeTableau() + void initializeTableau( const List &initialBasicVariables ) { TS_ASSERT( !initializeTableauCalled ); initializeTableauCalled = true; + + for ( const auto &basic : initialBasicVariables ) + lastBasicVariables.insert( basic ); } Map nextValues; diff --git a/src/engine/tests/Test_Tableau.h b/src/engine/tests/Test_Tableau.h index 43f61f149..1cc89c86d 100644 --- a/src/engine/tests/Test_Tableau.h +++ b/src/engine/tests/Test_Tableau.h @@ -162,10 +162,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); // All non-basics are set to lower bounds. TS_ASSERT_EQUALS( tableau->getValue( 0 ), 1.0 ); @@ -209,10 +207,6 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - MockVariableWatcher watcher1; MockVariableWatcher watcher2; @@ -220,7 +214,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->registerToWatchVariable( &watcher1, 5 ) ); TS_ASSERT_THROWS_NOTHING( tableau->registerToWatchVariable( &watcher2, 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); // The basic values get computed, so the watchers should be called @@ -270,10 +265,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 4 ), Tableau::BELOW_LB ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 5 ), Tableau::BETWEEN ); @@ -324,10 +317,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 412 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 4 ), Tableau::BETWEEN ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 5 ), Tableau::BELOW_LB ); @@ -378,10 +369,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 4 ), Tableau::BELOW_LB ); TS_ASSERT_EQUALS( tableau->getBasicStatus( 5 ), Tableau::BETWEEN ); @@ -487,10 +476,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_THROWS_NOTHING( tableau->computeCostFunction() ); costFunctionManager.nextCostFunction = new double[4]; @@ -547,10 +534,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_THROWS_NOTHING( tableau->computeCostFunction() ); costFunctionManager.nextCostFunction = new double[4]; @@ -621,10 +606,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TableauRow row( 4 ); @@ -825,10 +808,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_THROWS_NOTHING( tableau->computeAssignment() ); @@ -1028,10 +1009,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 200 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 202 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_THROWS_NOTHING( tableau->computeCostFunction() ); costFunctionManager.nextCostFunction = new double[4]; @@ -1149,10 +1128,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); // Do a pivot to shuffle the basis tableau->setEnteringVariableIndex( 2u ); @@ -1304,10 +1281,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); for ( unsigned i = 0; i < 4; ++i ) { @@ -1384,10 +1359,8 @@ class TableauTestSuite : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( tableau->setLowerBound( 6, 400 ) ); TS_ASSERT_THROWS_NOTHING( tableau->setUpperBound( 6, 402 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 4 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 5 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->markAsBasic( 6 ) ); - TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau() ); + List basics = { 4, 5, 6 }; + TS_ASSERT_THROWS_NOTHING( tableau->initializeTableau( basics ) ); TS_ASSERT_THROWS_NOTHING( tableau->computeCostFunction() ); costFunctionManager.nextCostFunction = new double[4]; From 08153234c25fcc6b85a13e7b013eec05e7dcd1c5 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 23 May 2018 09:49:24 +0300 Subject: [PATCH 11/14] make RowBoundTightener remember the tableau, instead of passing the tableau again and again every time --- src/engine/AutoRowBoundTightener.h | 4 +- src/engine/Engine.cpp | 23 ++--- src/engine/IRowBoundTightener.h | 14 +-- src/engine/ITableau.h | 6 ++ src/engine/RowBoundTightener.cpp | 97 ++++++++++--------- src/engine/RowBoundTightener.h | 32 +++--- src/engine/T/RowBoundTightenerFactory.h | 8 +- src/engine/real/RowBoundTightenerFactory.cpp | 4 +- src/engine/tests/MockRowBoundTightener.h | 18 ++-- .../tests/MockRowBoundTightenerFactory.h | 4 +- src/engine/tests/Test_RowBoundTightener.h | 36 +++---- 11 files changed, 133 insertions(+), 113 deletions(-) diff --git a/src/engine/AutoRowBoundTightener.h b/src/engine/AutoRowBoundTightener.h index 100c42b38..3bd26bf08 100644 --- a/src/engine/AutoRowBoundTightener.h +++ b/src/engine/AutoRowBoundTightener.h @@ -19,9 +19,9 @@ class AutoRowBoundTightener { public: - AutoRowBoundTightener() + AutoRowBoundTightener( const ITableau &tableau ) { - _rowBoundTightener = T::createRowBoundTightener(); + _rowBoundTightener = T::createRowBoundTightener( tableau ); } ~AutoRowBoundTightener() diff --git a/src/engine/Engine.cpp b/src/engine/Engine.cpp index f195294c8..839225269 100644 --- a/src/engine/Engine.cpp +++ b/src/engine/Engine.cpp @@ -25,7 +25,8 @@ #include "TimeUtils.h" Engine::Engine() - : _smtCore( this ) + : _rowBoundTightener( *_tableau ) + , _smtCore( this ) , _numPlConstraintsDisabledByValidSplits( 0 ) , _preprocessingEnabled( false ) , _work( NULL ) @@ -407,7 +408,7 @@ void Engine::performSimplexStep() if ( !fakePivot ) { _tableau->computePivotRow(); - _rowBoundTightener->examinePivotRow( _tableau ); + _rowBoundTightener->examinePivotRow(); } // Perform the actual pivot @@ -569,7 +570,7 @@ bool Engine::processInputQuery( InputQuery &inputQuery, bool preprocess ) } _tableau->registerToWatchAllVariables( _rowBoundTightener ); - _rowBoundTightener->initialize( _tableau ); + _rowBoundTightener->initialize(); _plConstraints = _preprocessedQuery.getPiecewiseLinearConstraints(); for ( const auto &constraint : _plConstraints ) @@ -701,7 +702,7 @@ void Engine::restoreState( const EngineState &state ) _numPlConstraintsDisabledByValidSplits = state._numPlConstraintsDisabledByValidSplits; // Make sure the data structures are initialized to the correct size - _rowBoundTightener->initialize( _tableau ); + _rowBoundTightener->initialize(); adjustWorkMemorySize(); _activeEntryStrategy->resizeHook( _tableau ); _costFunctionManager->initialize(); @@ -745,7 +746,7 @@ void Engine::applySplit( const PiecewiseLinearCaseSplit &split ) adjustWorkMemorySize(); - _rowBoundTightener->initialize( _tableau ); + _rowBoundTightener->initialize(); for ( auto &bound : bounds ) { @@ -868,7 +869,7 @@ void Engine::tightenBoundsOnConstraintMatrix() if ( _statistics.getNumMainLoopIterations() % GlobalConfiguration::BOUND_TIGHTING_ON_CONSTRAINT_MATRIX_FREQUENCY == 0 ) { - _rowBoundTightener->examineConstraintMatrix( _tableau, true ); + _rowBoundTightener->examineConstraintMatrix( true ); _statistics.incNumBoundTighteningOnConstraintMatrix(); } @@ -887,15 +888,15 @@ void Engine::explicitBasisBoundTightening() switch ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE ) { case GlobalConfiguration::USE_BASIS_MATRIX: - _rowBoundTightener->examineBasisMatrix( _tableau, saturation ); + _rowBoundTightener->examineBasisMatrix( saturation ); break; case GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX: - _rowBoundTightener->examineInvertedBasisMatrix( _tableau, saturation ); + _rowBoundTightener->examineInvertedBasisMatrix( saturation ); break; case GlobalConfiguration::USE_IMPLICIT_INVERTED_BASIS_MATRIX: - _rowBoundTightener->examineImplicitInvertedBasisMatrix( _tableau, saturation ); + _rowBoundTightener->examineImplicitInvertedBasisMatrix( saturation ); break; } @@ -916,7 +917,7 @@ void Engine::performPrecisionRestoration( PrecisionRestorer::RestoreBasics resto _statistics.addTimeForPrecisionRestoration( TimeUtils::timePassed( start, end ) ); _statistics.incNumPrecisionRestorations(); - _rowBoundTightener->clear( _tableau ); + _rowBoundTightener->clear(); // debug double after = _degradationChecker.computeDegradation( *_tableau ); @@ -936,7 +937,7 @@ void Engine::performPrecisionRestoration( PrecisionRestorer::RestoreBasics resto _statistics.addTimeForPrecisionRestoration( TimeUtils::timePassed( start, end ) ); _statistics.incNumPrecisionRestorations(); - _rowBoundTightener->clear( _tableau ); + _rowBoundTightener->clear(); // debug double afterSecond = _degradationChecker.computeDegradation( *_tableau ); diff --git a/src/engine/IRowBoundTightener.h b/src/engine/IRowBoundTightener.h index a78d2af8b..fb19cfec9 100644 --- a/src/engine/IRowBoundTightener.h +++ b/src/engine/IRowBoundTightener.h @@ -25,12 +25,12 @@ class IRowBoundTightener : public ITableau::VariableWatcher Allocate internal work memory according to the tableau size and initialize tightest lower/upper bounds using the talbeau. */ - virtual void initialize( const ITableau &tableau ) = 0; + virtual void initialize() = 0; /* Clear all learned bounds, without reallocating memory. */ - virtual void clear( const ITableau &tableau ) = 0; + virtual void clear() = 0; /* Derive and enqueue new bounds for all varaibles, using the @@ -38,7 +38,7 @@ class IRowBoundTightener : public ITableau::VariableWatcher tableau. Can also do this until saturation, meaning that we continue until no new bounds are learned. */ - virtual void examineBasisMatrix( const ITableau &tableau, bool untilSaturation ) = 0; + virtual void examineBasisMatrix( bool untilSaturation ) = 0; /* Derive and enqueue new bounds for all varaibles, using the @@ -46,7 +46,7 @@ class IRowBoundTightener : public ITableau::VariableWatcher through the tableau. Can also do this until saturation, meaning that we continue until no new bounds are learned. */ - virtual void examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ) = 0; + virtual void examineInvertedBasisMatrix( bool untilSaturation ) = 0; /* Derive and enqueue new bounds for all varaibles, implicitly using the @@ -55,7 +55,7 @@ class IRowBoundTightener : public ITableau::VariableWatcher 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; + virtual void examineImplicitInvertedBasisMatrix( bool untilSaturation ) = 0; /* Derive and enqueue new bounds for all varaibles, using the @@ -63,14 +63,14 @@ class IRowBoundTightener : public ITableau::VariableWatcher also do this until saturation, meaning that we continue until no new bounds are learned. */ - virtual void examineConstraintMatrix( const ITableau &tableau, bool untilSaturation ) = 0; + virtual void examineConstraintMatrix( bool untilSaturation ) = 0; /* Derive and enqueue new bounds immedaitely following a pivot operation in the given tableau. The tightening is performed for the entering variable (which is now basic). */ - virtual void examinePivotRow( ITableau &tableau ) = 0; + virtual void examinePivotRow() = 0; /* Get the tightenings entailed by the constraint. diff --git a/src/engine/ITableau.h b/src/engine/ITableau.h index 127a1de38..e6fd0a342 100644 --- a/src/engine/ITableau.h +++ b/src/engine/ITableau.h @@ -67,6 +67,12 @@ class ITableau */ virtual void notifyLowerBound( unsigned /* variable */, double /* bound */ ) {} virtual void notifyUpperBound( unsigned /* variable */, double /* bound */ ) {} + + /* + This callback will be invoked when the tableau size changes, + typically when new variables are added. + */ + virtual void notifyDimensionChange( unsigned /* m */, unsigned /* n */ ) {} }; virtual void registerToWatchAllVariables( VariableWatcher *watcher ) = 0; diff --git a/src/engine/RowBoundTightener.cpp b/src/engine/RowBoundTightener.cpp index 85f6ca35c..528a169c3 100644 --- a/src/engine/RowBoundTightener.cpp +++ b/src/engine/RowBoundTightener.cpp @@ -16,8 +16,9 @@ #include "RowBoundTightener.h" #include "Statistics.h" -RowBoundTightener::RowBoundTightener() - : _lowerBounds( NULL ) +RowBoundTightener::RowBoundTightener( const ITableau &tableau ) + : _tableau( tableau ) + , _lowerBounds( NULL ) , _upperBounds( NULL ) , _tightenedLower( NULL ) , _tightenedUpper( NULL ) @@ -30,12 +31,12 @@ RowBoundTightener::RowBoundTightener() { } -void RowBoundTightener::initialize( const ITableau &tableau ) +void RowBoundTightener::initialize() { freeMemoryIfNeeded(); - _n = tableau.getN(); - _m = tableau.getM(); + _n = _tableau.getN(); + _m = _tableau.getM(); _lowerBounds = new double[_n]; if ( !_lowerBounds ) @@ -58,8 +59,8 @@ void RowBoundTightener::initialize( const ITableau &tableau ) for ( unsigned i = 0; i < _n; ++i ) { - _lowerBounds[i] = tableau.getLowerBound( i ); - _upperBounds[i] = tableau.getUpperBound( i ); + _lowerBounds[i] = _tableau.getLowerBound( i ); + _upperBounds[i] = _tableau.getUpperBound( i ); } if ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE == @@ -84,15 +85,15 @@ void RowBoundTightener::initialize( const ITableau &tableau ) _ciSign = new char[_n]; } -void RowBoundTightener::clear( const ITableau &tableau ) +void RowBoundTightener::clear() { std::fill( _tightenedLower, _tightenedLower + _n, false ); std::fill( _tightenedUpper, _tightenedUpper + _n, false ); for ( unsigned i = 0; i < _n; ++i ) { - _lowerBounds[i] = tableau.getLowerBound( i ); - _upperBounds[i] = tableau.getUpperBound( i ); + _lowerBounds[i] = _tableau.getLowerBound( i ); + _upperBounds[i] = _tableau.getUpperBound( i ); } } @@ -160,7 +161,7 @@ void RowBoundTightener::freeMemoryIfNeeded() } } -void RowBoundTightener::examineImplicitInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ) +void RowBoundTightener::examineImplicitInvertedBasisMatrix( bool untilSaturation ) { /* Roughly (the dimensions don't add up): @@ -169,20 +170,20 @@ void RowBoundTightener::examineImplicitInvertedBasisMatrix( const ITableau &tabl */ // Find z = inv(B) * b, by solving the forward transformation Bz = b - tableau.forwardTransformation( tableau.getRightHandSide(), _z ); + _tableau.forwardTransformation( _tableau.getRightHandSide(), _z ); for ( unsigned i = 0; i < _m; ++i ) { _rows[i]->_scalar = _z[i]; - _rows[i]->_lhs = tableau.basicIndexToVariable( 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 ); + unsigned nonBasic = _tableau.nonBasicIndexToVariable( i ); + const double *ANColumn = _tableau.getAColumn( nonBasic ); + _tableau.forwardTransformation( ANColumn, _z ); for ( unsigned j = 0; j < _m; ++j ) { @@ -197,7 +198,7 @@ void RowBoundTightener::examineImplicitInvertedBasisMatrix( const ITableau &tabl unsigned newBoundsLearned; do { - newBoundsLearned = onePassOverInvertedBasisRows( tableau ); + newBoundsLearned = onePassOverInvertedBasisRows(); if ( _statistics && ( newBoundsLearned > 0 ) ) _statistics->incNumTighteningsFromExplicitBasis( newBoundsLearned ); @@ -205,7 +206,7 @@ void RowBoundTightener::examineImplicitInvertedBasisMatrix( const ITableau &tabl while ( untilSaturation && ( newBoundsLearned > 0 ) ); } -void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ) +void RowBoundTightener::examineInvertedBasisMatrix( bool untilSaturation ) { /* Roughly (the dimensions don't add up): @@ -215,8 +216,8 @@ void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, boo We compute one row at a time. */ - const double *b = tableau.getRightHandSide(); - const double *invB = tableau.getInverseBasisMatrix(); + const double *b = _tableau.getRightHandSide(); + const double *invB = _tableau.getInverseBasisMatrix(); try { @@ -231,18 +232,18 @@ void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, boo // Now update the row's coefficients for basic variable i for ( unsigned j = 0; j < _n - _m; ++j ) { - row->_row[j]._var = tableau.nonBasicIndexToVariable( 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 ); + 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] ); } // Store the lhs variable - row->_lhs = tableau.basicIndexToVariable( i ); + row->_lhs = _tableau.basicIndexToVariable( i ); } // We now have all the rows, can use them for tightening. @@ -252,7 +253,7 @@ void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, boo unsigned newBoundsLearned; do { - newBoundsLearned = onePassOverInvertedBasisRows( tableau ); + newBoundsLearned = onePassOverInvertedBasisRows(); if ( _statistics && ( newBoundsLearned > 0 ) ) _statistics->incNumTighteningsFromExplicitBasis( newBoundsLearned ); @@ -268,17 +269,17 @@ void RowBoundTightener::examineInvertedBasisMatrix( const ITableau &tableau, boo delete[] invB; } -unsigned RowBoundTightener::onePassOverInvertedBasisRows( const ITableau &tableau ) +unsigned RowBoundTightener::onePassOverInvertedBasisRows() { unsigned newBounds = 0; for ( unsigned i = 0; i < _m; ++i ) - newBounds += tightenOnSingleInvertedBasisRow( tableau, *( _rows[i] ) ); + newBounds += tightenOnSingleInvertedBasisRow( *( _rows[i] ) ); return newBounds; } -unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const ITableau &tableau, const TableauRow &row ) +unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const TableauRow &row ) { /* A row is of the form @@ -287,8 +288,8 @@ unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const ITableau &tab We wish to tighten once for y, but also once for every x. */ - unsigned n = tableau.getN(); - unsigned m = tableau.getM(); + unsigned n = _tableau.getN(); + unsigned m = _tableau.getM(); unsigned result = 0; @@ -447,7 +448,7 @@ unsigned RowBoundTightener::tightenOnSingleInvertedBasisRow( const ITableau &tab return result; } -void RowBoundTightener::examineBasisMatrix( const ITableau &tableau, bool untilSaturation ) +void RowBoundTightener::examineBasisMatrix( bool untilSaturation ) { unsigned newBoundsLearned; @@ -457,7 +458,7 @@ void RowBoundTightener::examineBasisMatrix( const ITableau &tableau, bool untilS */ do { - newBoundsLearned = onePassOverBasisMatrix( tableau ); + newBoundsLearned = onePassOverBasisMatrix(); if ( _statistics && ( newBoundsLearned > 0 ) ) _statistics->incNumTighteningsFromExplicitBasis( newBoundsLearned ); @@ -465,12 +466,12 @@ void RowBoundTightener::examineBasisMatrix( const ITableau &tableau, bool untilS while ( untilSaturation && ( newBoundsLearned > 0 ) ); } -unsigned RowBoundTightener::onePassOverBasisMatrix( const ITableau &tableau ) +unsigned RowBoundTightener::onePassOverBasisMatrix() { unsigned newBounds = 0; List basisEquations; - tableau.getBasisEquations( basisEquations ); + _tableau.getBasisEquations( basisEquations ); for ( const auto &equation : basisEquations ) for ( const auto &addend : equation->_addends ) newBounds += tightenOnSingleEquation( *equation, addend ); @@ -552,7 +553,7 @@ unsigned RowBoundTightener::tightenOnSingleEquation( Equation &equation, return result; } -void RowBoundTightener::examineConstraintMatrix( const ITableau &tableau, bool untilSaturation ) +void RowBoundTightener::examineConstraintMatrix( bool untilSaturation ) { unsigned newBoundsLearned; @@ -562,7 +563,7 @@ void RowBoundTightener::examineConstraintMatrix( const ITableau &tableau, bool u */ do { - newBoundsLearned = onePassOverConstraintMatrix( tableau ); + newBoundsLearned = onePassOverConstraintMatrix(); if ( _statistics && ( newBoundsLearned > 0 ) ) _statistics->incNumTighteningsFromConstraintMatrix( newBoundsLearned ); @@ -570,19 +571,19 @@ void RowBoundTightener::examineConstraintMatrix( const ITableau &tableau, bool u while ( untilSaturation && ( newBoundsLearned > 0 ) ); } -unsigned RowBoundTightener::onePassOverConstraintMatrix( const ITableau &tableau ) +unsigned RowBoundTightener::onePassOverConstraintMatrix() { unsigned result = 0; - unsigned m = tableau.getM(); + unsigned m = _tableau.getM(); for ( unsigned i = 0; i < m; ++i ) - result += tightenOnSingleConstraintRow( tableau, i ); + result += tightenOnSingleConstraintRow( i ); return result; } -unsigned RowBoundTightener::tightenOnSingleConstraintRow( const ITableau &tableau, unsigned row ) +unsigned RowBoundTightener::tightenOnSingleConstraintRow( unsigned row ) { /* The cosntraint matrix A satisfies Ax = b. @@ -594,13 +595,13 @@ unsigned RowBoundTightener::tightenOnSingleConstraintRow( const ITableau &tablea sum ci xi - b */ - unsigned n = tableau.getN(); - unsigned m = tableau.getM(); + unsigned n = _tableau.getN(); + unsigned m = _tableau.getM(); unsigned result = 0; - const double *A = tableau.getA(); - const double *b = tableau.getRightHandSide(); + const double *A = _tableau.getA(); + const double *b = _tableau.getRightHandSide(); double ci; @@ -723,13 +724,13 @@ unsigned RowBoundTightener::tightenOnSingleConstraintRow( const ITableau &tablea return result; } -void RowBoundTightener::examinePivotRow( ITableau &tableau ) +void RowBoundTightener::examinePivotRow() { if ( _statistics ) _statistics->incNumRowsExaminedByRowTightener(); - const TableauRow &row = *tableau.getPivotRow(); - unsigned newBoundsLearned = tightenOnSingleInvertedBasisRow( tableau, row ); + const TableauRow &row( *_tableau.getPivotRow() ); + unsigned newBoundsLearned = tightenOnSingleInvertedBasisRow( row ); if ( _statistics && ( newBoundsLearned > 0 ) ) _statistics->incNumTighteningsFromRows( newBoundsLearned ); @@ -776,6 +777,10 @@ void RowBoundTightener::notifyUpperBound( unsigned variable, double bound ) } } +void RowBoundTightener::notifyDimensionChange( unsigned /* m */ , unsigned /* n */ ) +{ +} + // // Local Variables: // compile-command: "make -C ../.. " diff --git a/src/engine/RowBoundTightener.h b/src/engine/RowBoundTightener.h index 153c526e0..78093e54d 100644 --- a/src/engine/RowBoundTightener.h +++ b/src/engine/RowBoundTightener.h @@ -23,19 +23,19 @@ class RowBoundTightener : public IRowBoundTightener { public: - RowBoundTightener(); + RowBoundTightener( const ITableau &tableau ); ~RowBoundTightener(); /* Allocate internal work memory according to the tableau size and initialize tightest lower/upper bounds using the talbeau. */ - void initialize( const ITableau &tableau ); + void initialize(); /* Clear all learned bounds, without reallocating memory. */ - void clear( const ITableau &tableau ); + void clear(); /* Callbacks from the Tableau, to inform of bounds tightened by, @@ -44,13 +44,18 @@ class RowBoundTightener : public IRowBoundTightener void notifyLowerBound( unsigned variable, double bound ); void notifyUpperBound( unsigned variable, double bound ); + /* + Callbacks from the Tableau, to inform of a change in dimesions + */ + void notifyDimensionChange( unsigned m, unsigned n ); + /* Derive and enqueue new bounds for all varaibles, using the explicit basis matrix B0 that should be available through the tableau. Can also do this until saturation, meaning that we continue until no new bounds are learned. */ - void examineBasisMatrix( const ITableau &tableau, bool untilSaturation ); + void examineBasisMatrix( bool untilSaturation ); /* Derive and enqueue new bounds for all varaibles, using the @@ -58,7 +63,7 @@ class RowBoundTightener : public IRowBoundTightener through the tableau. Can also do this until saturation, meaning that we continue until no new bounds are learned. */ - void examineInvertedBasisMatrix( const ITableau &tableau, bool untilSaturation ); + void examineInvertedBasisMatrix( bool untilSaturation ); /* Derive and enqueue new bounds for all varaibles, implicitly using the @@ -67,7 +72,7 @@ class RowBoundTightener : public IRowBoundTightener 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 ); + void examineImplicitInvertedBasisMatrix( bool untilSaturation ); /* Derive and enqueue new bounds for all varaibles, using the @@ -75,14 +80,14 @@ class RowBoundTightener : public IRowBoundTightener also do this until saturation, meaning that we continue until no new bounds are learned. */ - void examineConstraintMatrix( const ITableau &tableau, bool untilSaturation ); + void examineConstraintMatrix( bool untilSaturation ); /* Derive and enqueue new bounds immedaitely following a pivot operation in the given tableau. The tightening is performed for the entering variable (which is now basic). */ - void examinePivotRow( ITableau &tableau ); + void examinePivotRow(); /* Get the tightenings entailed by the constraint. @@ -95,6 +100,7 @@ class RowBoundTightener : public IRowBoundTightener void setStatistics( Statistics *statistics ); private: + const ITableau &_tableau; unsigned _n; unsigned _m; @@ -132,7 +138,7 @@ class RowBoundTightener : public IRowBoundTightener Do a single pass over the basis matrix and derive any tighter bounds. Return the number of new bounds are learned. */ - unsigned onePassOverBasisMatrix( const ITableau &tableau ); + unsigned onePassOverBasisMatrix(); /* Process the basis row and attempt to derive tighter @@ -146,27 +152,27 @@ class RowBoundTightener : public IRowBoundTightener Do a single pass over the constraint matrix and derive any tighter bounds. Return the number of new bounds learned. */ - unsigned onePassOverConstraintMatrix( const ITableau &tableau ); + unsigned onePassOverConstraintMatrix(); /* Process the tableau row and attempt to derive tighter lower/upper bounds for the specified variable. Return the number of tighter bounds found. */ - unsigned tightenOnSingleConstraintRow( const ITableau &tableau, unsigned row ); + unsigned tightenOnSingleConstraintRow( unsigned row ); /* Do a single pass over the inverted basis rows and derive any tighter bounds. Return the number of new bounds learned. */ - unsigned onePassOverInvertedBasisRows( const ITableau &tableau ); + unsigned onePassOverInvertedBasisRows(); /* Process the inverted basis row and attempt to derive tighter lower/upper bounds for the specified variable. Return the number of tighter bounds found. */ - unsigned tightenOnSingleInvertedBasisRow( const ITableau &tableau, const TableauRow &row ); + unsigned tightenOnSingleInvertedBasisRow( const TableauRow &row ); }; #endif // __RowBoundTightener_h__ diff --git a/src/engine/T/RowBoundTightenerFactory.h b/src/engine/T/RowBoundTightenerFactory.h index 013b03757..00b3c7111 100644 --- a/src/engine/T/RowBoundTightenerFactory.h +++ b/src/engine/T/RowBoundTightenerFactory.h @@ -16,20 +16,20 @@ #include "cxxtest/Mock.h" class IRowBoundTightener; -class ISelector; +class ITableau; namespace T { - IRowBoundTightener *createRowBoundTightener(); + IRowBoundTightener *createRowBoundTightener( const ITableau &tableau ); void discardRowBoundTightener( IRowBoundTightener *rowBoundTightener ); } CXXTEST_SUPPLY( createRowBoundTightener, IRowBoundTightener *, createRowBoundTightener, - (), + ( const ITableau &tableau ), T::createRowBoundTightener, - () ); + ( tableau ) ); CXXTEST_SUPPLY_VOID( discardRowBoundTightener, discardRowBoundTightener, diff --git a/src/engine/real/RowBoundTightenerFactory.cpp b/src/engine/real/RowBoundTightenerFactory.cpp index d7c1cf04a..5f06ca585 100644 --- a/src/engine/real/RowBoundTightenerFactory.cpp +++ b/src/engine/real/RowBoundTightenerFactory.cpp @@ -14,9 +14,9 @@ namespace T { - IRowBoundTightener *createRowBoundTightener() + IRowBoundTightener *createRowBoundTightener( const ITableau &tableau ) { - return new RowBoundTightener(); + return new RowBoundTightener( tableau ); } void discardRowBoundTightener( IRowBoundTightener *rowBoundTightener ) diff --git a/src/engine/tests/MockRowBoundTightener.h b/src/engine/tests/MockRowBoundTightener.h index 0b524aec3..840b17723 100644 --- a/src/engine/tests/MockRowBoundTightener.h +++ b/src/engine/tests/MockRowBoundTightener.h @@ -32,11 +32,13 @@ class MockRowBoundTightener : public IRowBoundTightener bool wasCreated; bool wasDiscarded; + const ITableau *lastTableau; - void mockConstructor() + void mockConstructor( const ITableau &tableau ) { TS_ASSERT( !wasCreated ); wasCreated = true; + lastTableau = &tableau; } void mockDestructor() @@ -47,21 +49,21 @@ class MockRowBoundTightener : public IRowBoundTightener } bool initializeWasCalled; - void initialize( const ITableau &/* tableau */ ) + void initialize() { initializeWasCalled = true; } - void clear( const ITableau &/* tableau */ ) {} + void clear() {} void notifyLowerBound( unsigned /* variable */, double /* bound */ ) {} void notifyUpperBound( unsigned /* variable */, double /* bound */ ) {} - void examineBasisMatrix( const ITableau &/* tableau */, bool /* untilSaturation */ ) {} - void examineInvertedBasisMatrix( const ITableau &/* tableau */, bool /* untilSaturation */ ) {} - void examineConstraintMatrix( const ITableau &/* tableau */, bool /* untilSaturation */ ) {} - void examinePivotRow( ITableau &/* tableau */ ) {} + void examineBasisMatrix( bool /* untilSaturation */ ) {} + void examineInvertedBasisMatrix( bool /* untilSaturation */ ) {} + void examineConstraintMatrix( bool /* untilSaturation */ ) {} + void examinePivotRow() {} void getRowTightenings( List &/* tightenings */ ) const {} void setStatistics( Statistics */* statistics */ ) {} - void examineImplicitInvertedBasisMatrix( const ITableau &/* tableau */, bool /* untilSaturation */ ) {} + void examineImplicitInvertedBasisMatrix( bool /* untilSaturation */ ) {} }; #endif // __MockRowBoundTightener_h__ diff --git a/src/engine/tests/MockRowBoundTightenerFactory.h b/src/engine/tests/MockRowBoundTightenerFactory.h index dc3d5065a..39e6c4b92 100644 --- a/src/engine/tests/MockRowBoundTightenerFactory.h +++ b/src/engine/tests/MockRowBoundTightenerFactory.h @@ -31,9 +31,9 @@ class MockRowBoundTightenerFactory : } } - IRowBoundTightener *createRowBoundTightener() + IRowBoundTightener *createRowBoundTightener( const ITableau &tableau ) { - mockRowBoundTightener.mockConstructor(); + mockRowBoundTightener.mockConstructor( tableau ); return &mockRowBoundTightener; } diff --git a/src/engine/tests/Test_RowBoundTightener.h b/src/engine/tests/Test_RowBoundTightener.h index b8ebdeb48..81d494f83 100644 --- a/src/engine/tests/Test_RowBoundTightener.h +++ b/src/engine/tests/Test_RowBoundTightener.h @@ -40,7 +40,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_pivot_row__both_bounds_tightened() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 2, 5 ); @@ -60,7 +60,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, -100 ); tableau->setUpperBound( 4, 100 ); - tightener.initialize( *tableau ); + tightener.initialize(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -76,7 +76,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite // x0 entering, x4 leaving // x0 = 1 - x1 + 2 x2 - x4 - TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow( *tableau ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow() ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); @@ -102,7 +102,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_pivot_row__just_upper_tightend() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 2, 5 ); @@ -122,7 +122,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, 0 ); tableau->setUpperBound( 4, 0 ); - tightener.initialize( *tableau ); + tightener.initialize(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -134,7 +134,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->nextPivotRow = &row; - TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow( *tableau ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow() ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); @@ -153,7 +153,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_pivot_row__just_lower_tightend() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 2, 5 ); @@ -168,7 +168,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, 0 ); tableau->setUpperBound( 4, 0 ); - tightener.initialize( *tableau ); + tightener.initialize(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -180,7 +180,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->nextPivotRow = &row; - TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow( *tableau ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow() ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); @@ -198,7 +198,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_pivot_row__nothing_tightened() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 2, 5 ); @@ -213,7 +213,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, -100 ); tableau->setUpperBound( 4, 100 ); - tightener.initialize( *tableau ); + tightener.initialize(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -225,7 +225,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->nextPivotRow = &row; - TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow( *tableau ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examinePivotRow() ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); @@ -235,7 +235,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_examine_constraint_matrix_single_equation() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 1, 5 ); @@ -267,7 +267,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite x1 >= 1.5 */ - tightener.initialize( *tableau ); + tightener.initialize(); double A[] = { 1, -2, 0, 1, 2 }; double b[] = { 1 }; @@ -275,7 +275,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->A = A; tableau->b = b; - TS_ASSERT_THROWS_NOTHING( tightener.examineConstraintMatrix( *tableau, false ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examineConstraintMatrix( false ) ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); @@ -296,7 +296,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite void test_examine_constraint_matrix_multiple_equations() { - RowBoundTightener tightener; + RowBoundTightener tightener( *tableau ); tableau->setDimensions( 2, 5 ); @@ -332,7 +332,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite Equation 2 gives us that x2 <= 2, >= 1 */ - tightener.initialize( *tableau ); + tightener.initialize(); double A[] = { 1, 0, -2, -2, @@ -345,7 +345,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->A = A; tableau->b = b; - TS_ASSERT_THROWS_NOTHING( tightener.examineConstraintMatrix( *tableau, false ) ); + TS_ASSERT_THROWS_NOTHING( tightener.examineConstraintMatrix( false ) ); List tightenings; TS_ASSERT_THROWS_NOTHING( tightener.getRowTightenings( tightenings ) ); From e4e50895dcc8b9fedf81d40c7e311be63ab6145d Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 23 May 2018 12:36:45 +0300 Subject: [PATCH 12/14] Tableau's addRow() informs objects registered as "resizeWatchers" that the tableau dimensions have changed. Tableau handles the case where the basis factorization fails after a row is added, by choosing a new set of basic variables using gaussian elimination of the constraint matrix --- src/engine/Engine.cpp | 8 ++- src/engine/IRowBoundTightener.h | 12 ++-- src/engine/ITableau.h | 6 ++ src/engine/RowBoundTightener.cpp | 24 ++++--- src/engine/RowBoundTightener.h | 12 ++-- src/engine/Tableau.cpp | 77 +++++++++++++++++++---- src/engine/Tableau.h | 13 ++++ src/engine/tests/MockRowBoundTightener.h | 12 ++-- src/engine/tests/MockTableau.h | 6 ++ src/engine/tests/Test_Engine.h | 5 +- src/engine/tests/Test_RowBoundTightener.h | 12 ++-- 11 files changed, 143 insertions(+), 44 deletions(-) diff --git a/src/engine/Engine.cpp b/src/engine/Engine.cpp index 839225269..c74df43f6 100644 --- a/src/engine/Engine.cpp +++ b/src/engine/Engine.cpp @@ -570,7 +570,9 @@ bool Engine::processInputQuery( InputQuery &inputQuery, bool preprocess ) } _tableau->registerToWatchAllVariables( _rowBoundTightener ); - _rowBoundTightener->initialize(); + _tableau->registerResizeWatcher( _rowBoundTightener ); + + _rowBoundTightener->setDimensions(); _plConstraints = _preprocessedQuery.getPiecewiseLinearConstraints(); for ( const auto &constraint : _plConstraints ) @@ -702,7 +704,7 @@ void Engine::restoreState( const EngineState &state ) _numPlConstraintsDisabledByValidSplits = state._numPlConstraintsDisabledByValidSplits; // Make sure the data structures are initialized to the correct size - _rowBoundTightener->initialize(); + _rowBoundTightener->setDimensions(); adjustWorkMemorySize(); _activeEntryStrategy->resizeHook( _tableau ); _costFunctionManager->initialize(); @@ -746,7 +748,7 @@ void Engine::applySplit( const PiecewiseLinearCaseSplit &split ) adjustWorkMemorySize(); - _rowBoundTightener->initialize(); + _rowBoundTightener->resetBounds(); for ( auto &bound : bounds ) { diff --git a/src/engine/IRowBoundTightener.h b/src/engine/IRowBoundTightener.h index fb19cfec9..dd6f1fceb 100644 --- a/src/engine/IRowBoundTightener.h +++ b/src/engine/IRowBoundTightener.h @@ -16,16 +16,20 @@ #include "ITableau.h" #include "Tightening.h" -class IRowBoundTightener : public ITableau::VariableWatcher +class IRowBoundTightener : public ITableau::VariableWatcher, public ITableau::ResizeWatcher { public: virtual ~IRowBoundTightener() {}; /* - Allocate internal work memory according to the tableau size and - initialize tightest lower/upper bounds using the talbeau. + Allocate internal work memory according to the tableau size. */ - virtual void initialize() = 0; + virtual void setDimensions() = 0; + + /* + Initialize tightest lower/upper bounds using the talbeau. + */ + virtual void resetBounds() = 0; /* Clear all learned bounds, without reallocating memory. diff --git a/src/engine/ITableau.h b/src/engine/ITableau.h index e6fd0a342..6fbbc4f48 100644 --- a/src/engine/ITableau.h +++ b/src/engine/ITableau.h @@ -67,7 +67,11 @@ class ITableau */ virtual void notifyLowerBound( unsigned /* variable */, double /* bound */ ) {} virtual void notifyUpperBound( unsigned /* variable */, double /* bound */ ) {} + }; + class ResizeWatcher + { + public: /* This callback will be invoked when the tableau size changes, typically when new variables are added. @@ -79,6 +83,8 @@ class ITableau virtual void registerToWatchVariable( VariableWatcher *watcher, unsigned variable ) = 0; virtual void unregisterToWatchVariable( VariableWatcher *watcher, unsigned variable ) = 0; + virtual void registerResizeWatcher( ResizeWatcher *watcher ) = 0; + virtual void registerCostFunctionManager( ICostFunctionManager *costFunctionManager ) = 0; virtual ~ITableau() {}; diff --git a/src/engine/RowBoundTightener.cpp b/src/engine/RowBoundTightener.cpp index 528a169c3..bb29edfc1 100644 --- a/src/engine/RowBoundTightener.cpp +++ b/src/engine/RowBoundTightener.cpp @@ -31,7 +31,7 @@ RowBoundTightener::RowBoundTightener( const ITableau &tableau ) { } -void RowBoundTightener::initialize() +void RowBoundTightener::setDimensions() { freeMemoryIfNeeded(); @@ -54,14 +54,7 @@ void RowBoundTightener::initialize() if ( !_tightenedUpper ) throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "RowBoundTightener::tightenedUpper" ); - std::fill( _tightenedLower, _tightenedLower + _n, false ); - std::fill( _tightenedUpper, _tightenedUpper + _n, false ); - - for ( unsigned i = 0; i < _n; ++i ) - { - _lowerBounds[i] = _tableau.getLowerBound( i ); - _upperBounds[i] = _tableau.getUpperBound( i ); - } + resetBounds(); if ( GlobalConfiguration::EXPLICIT_BASIS_BOUND_TIGHTENING_TYPE == GlobalConfiguration::COMPUTE_INVERTED_BASIS_MATRIX ) @@ -85,6 +78,18 @@ void RowBoundTightener::initialize() _ciSign = new char[_n]; } +void RowBoundTightener::resetBounds() +{ + std::fill( _tightenedLower, _tightenedLower + _n, false ); + std::fill( _tightenedUpper, _tightenedUpper + _n, false ); + + for ( unsigned i = 0; i < _n; ++i ) + { + _lowerBounds[i] = _tableau.getLowerBound( i ); + _upperBounds[i] = _tableau.getUpperBound( i ); + } +} + void RowBoundTightener::clear() { std::fill( _tightenedLower, _tightenedLower + _n, false ); @@ -779,6 +784,7 @@ void RowBoundTightener::notifyUpperBound( unsigned variable, double bound ) void RowBoundTightener::notifyDimensionChange( unsigned /* m */ , unsigned /* n */ ) { + setDimensions(); } // diff --git a/src/engine/RowBoundTightener.h b/src/engine/RowBoundTightener.h index 78093e54d..0c2d5c593 100644 --- a/src/engine/RowBoundTightener.h +++ b/src/engine/RowBoundTightener.h @@ -27,10 +27,14 @@ class RowBoundTightener : public IRowBoundTightener ~RowBoundTightener(); /* - Allocate internal work memory according to the tableau size and - initialize tightest lower/upper bounds using the talbeau. + Allocate internal work memory according to the tableau size */ - void initialize(); + void setDimensions(); + + /* + Initialize tightest lower/upper bounds using the talbeau. + */ + void resetBounds(); /* Clear all learned bounds, without reallocating memory. @@ -45,7 +49,7 @@ class RowBoundTightener : public IRowBoundTightener void notifyUpperBound( unsigned variable, double bound ); /* - Callbacks from the Tableau, to inform of a change in dimesions + Callback from the Tableau, to inform of a change in dimensions */ void notifyDimensionChange( unsigned m, unsigned n ); diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 4180d5f94..edd688a05 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -11,6 +11,7 @@ **/ #include "BasisFactorizationFactory.h" +#include "ConstraintMatrixAnalyzer.h" #include "Debug.h" #include "EntrySelectionStrategy.h" #include "Equation.h" @@ -1243,6 +1244,47 @@ void Tableau::addEquation( const Equation &equation ) // Invalidate the cost function, so that it is recomputed in the next iteration. _costFunctionManager->invalidateCostFunction(); + // All variables except the new one have finite bounds. Use this to compute + // finite bounds for the new variable. + double lb = equation._scalar; + double ub = equation._scalar; + double auxCoefficient; + double coefficient; + unsigned variable; + for ( const auto &addend : equation._addends ) + { + coefficient = addend._coefficient; + variable = addend._variable; + if ( variable == equation._auxVariable ) + { + auxCoefficient = coefficient; + } + else + { + if ( FloatUtils::isPositive( coefficient ) ) + { + lb -= coefficient * _upperBounds[variable]; + ub -= coefficient * _lowerBounds[variable]; + } + else + { + lb -= coefficient * _lowerBounds[variable]; + ub -= coefficient * _upperBounds[variable]; + } + } + + if ( FloatUtils::isPositive( auxCoefficient ) ) + { + setLowerBound( equation._auxVariable, lb / auxCoefficient ); + setUpperBound( equation._auxVariable, ub / auxCoefficient ); + } + else + { + setLowerBound( equation._auxVariable, ub / auxCoefficient ); + setUpperBound( equation._auxVariable, lb / auxCoefficient ); + } + } + // Populate the new row of A and b _b[_m - 1] = equation._scalar; for ( const auto &addend : equation._addends ) @@ -1251,8 +1293,8 @@ void Tableau::addEquation( const Equation &equation ) /* Attempt to make the auxiliary variable the new basic variable. This usually works. - If it doesn't, compute a new set of basic variables (which is more - computationally expensive) + If it doesn't, compute a new set of basic variables and re-initialize + the tableau (which is more computationally expensive) */ _basicIndexToVariable[_m - 1] = equation._auxVariable; _variableToIndex[equation._auxVariable] = _m - 1; @@ -1295,19 +1337,20 @@ void Tableau::addEquation( const Equation &equation ) } else { - // ConstraintMatrixAnalyzer analyzer; - // analyzer.analyze( _A, _m, _n ); - // List independentColumns = analyzer->getIndependentColumns(); + ConstraintMatrixAnalyzer analyzer; + analyzer.analyze( _A, _m, _n ); + List independentColumns = analyzer.getIndependentColumns(); - // Need to assign basic variables and initialize a whole bunch of stuff... + try + { + initializeTableau( independentColumns ); + } + catch ( MalformedBasisException & ) + { + log( "addEquation failed - could not refactorize basis" ); + throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); + } } - - // if ( !factorizationSuccessful ) - // { - // printf( "Out of candidates, terminating\n" ); - // log( "addEquation failed - could not refactorize basis" ); - // throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); - // } } void Tableau::addRow() @@ -1436,6 +1479,9 @@ void Tableau::addRow() _m = newM; _n = newN; _costFunctionManager->initialize(); + + for ( const auto &watcher : _resizeWatchers ) + watcher->notifyDimensionChange( _m, _n ); } void Tableau::registerToWatchVariable( VariableWatcher *watcher, unsigned variable ) @@ -1453,6 +1499,11 @@ void Tableau::registerToWatchAllVariables( VariableWatcher *watcher ) _globalWatchers.append( watcher ); } +void Tableau::registerResizeWatcher( ResizeWatcher *watcher ) +{ + _resizeWatchers.append( watcher ); +} + void Tableau::notifyVariableValue( unsigned variable, double value ) { for ( auto &watcher : _globalWatchers ) diff --git a/src/engine/Tableau.h b/src/engine/Tableau.h index ef45d13fd..fc286c67c 100644 --- a/src/engine/Tableau.h +++ b/src/engine/Tableau.h @@ -333,6 +333,11 @@ class Tableau : public ITableau, public IBasisFactorization::BasisColumnOracle void registerToWatchVariable( VariableWatcher *watcher, unsigned variable ); void unregisterToWatchVariable( VariableWatcher *watcher, unsigned variable ); + /* + Register to watch for tableau dimension changes. + */ + void registerResizeWatcher( ResizeWatcher *watcher ); + /* Register the cost function manager. */ @@ -388,10 +393,18 @@ class Tableau : public ITableau, public IBasisFactorization::BasisColumnOracle const double *getColumnOfBasis( unsigned column ) const; private: + /* + Variable watchers + */ typedef List VariableWatchers; Map _variableToWatchers; List _globalWatchers; + /* + Resize watchers + */ + List _resizeWatchers; + /* The dimensions of matrix A */ diff --git a/src/engine/tests/MockRowBoundTightener.h b/src/engine/tests/MockRowBoundTightener.h index 840b17723..67065dc11 100644 --- a/src/engine/tests/MockRowBoundTightener.h +++ b/src/engine/tests/MockRowBoundTightener.h @@ -23,7 +23,7 @@ class MockRowBoundTightener : public IRowBoundTightener wasCreated = false; wasDiscarded = false; - initializeWasCalled = false; + setDimensionsWasCalled = false; } ~MockRowBoundTightener() @@ -48,10 +48,14 @@ class MockRowBoundTightener : public IRowBoundTightener wasDiscarded = true; } - bool initializeWasCalled; - void initialize() + bool setDimensionsWasCalled; + void setDimensions() + { + setDimensionsWasCalled = true; + } + + void resetBounds() { - initializeWasCalled = true; } void clear() {} diff --git a/src/engine/tests/MockTableau.h b/src/engine/tests/MockTableau.h index fd8a05f0c..1915b8d42 100644 --- a/src/engine/tests/MockTableau.h +++ b/src/engine/tests/MockTableau.h @@ -432,6 +432,12 @@ class MockTableau : public ITableau { } + Set lastResizeWatchers; + void registerResizeWatcher( ResizeWatcher *watcher ) + { + lastResizeWatchers.insert( watcher ); + } + ICostFunctionManager *lastCostFunctionManager; void registerCostFunctionManager( ICostFunctionManager *costFunctionManager ) { diff --git a/src/engine/tests/Test_Engine.h b/src/engine/tests/Test_Engine.h index 48a1b9c74..e111cae6e 100644 --- a/src/engine/tests/Test_Engine.h +++ b/src/engine/tests/Test_Engine.h @@ -133,7 +133,10 @@ class EngineTestSuite : public CxxTest::TestSuite TS_ASSERT( tableau->initializeTableauCalled ); TS_ASSERT( costFunctionManager->initializeWasCalled ); - TS_ASSERT( rowTightener->initializeWasCalled ); + TS_ASSERT( rowTightener->setDimensionsWasCalled ); + + TS_ASSERT_EQUALS( tableau->lastResizeWatchers.size(), 1U ); + TS_ASSERT_EQUALS( *( tableau->lastResizeWatchers.begin() ), rowTightener ); TS_ASSERT_EQUALS( tableau->lastCostFunctionManager, costFunctionManager ); diff --git a/src/engine/tests/Test_RowBoundTightener.h b/src/engine/tests/Test_RowBoundTightener.h index 81d494f83..570f3aea6 100644 --- a/src/engine/tests/Test_RowBoundTightener.h +++ b/src/engine/tests/Test_RowBoundTightener.h @@ -60,7 +60,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, -100 ); tableau->setUpperBound( 4, 100 ); - tightener.initialize(); + tightener.setDimensions(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -122,7 +122,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, 0 ); tableau->setUpperBound( 4, 0 ); - tightener.initialize(); + tightener.setDimensions(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -168,7 +168,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, 0 ); tableau->setUpperBound( 4, 0 ); - tightener.initialize(); + tightener.setDimensions(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -213,7 +213,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite tableau->setLowerBound( 4, -100 ); tableau->setUpperBound( 4, 100 ); - tightener.initialize(); + tightener.setDimensions(); TableauRow row( 3 ); // 1 - x0 - x1 + 2x2 @@ -267,7 +267,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite x1 >= 1.5 */ - tightener.initialize(); + tightener.setDimensions(); double A[] = { 1, -2, 0, 1, 2 }; double b[] = { 1 }; @@ -332,7 +332,7 @@ class RowBoundTightenerTestSuite : public CxxTest::TestSuite Equation 2 gives us that x2 <= 2, >= 1 */ - tightener.initialize(); + tightener.setDimensions(); double A[] = { 1, 0, -2, -2, From b023d476de48d0a58957e0bd65abbe0338826fe9 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 23 May 2018 12:41:44 +0300 Subject: [PATCH 13/14] oops --- src/engine/Tableau.cpp | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index edd688a05..f8cd56ead 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -1251,10 +1251,12 @@ void Tableau::addEquation( const Equation &equation ) double auxCoefficient; double coefficient; unsigned variable; + for ( const auto &addend : equation._addends ) { coefficient = addend._coefficient; variable = addend._variable; + if ( variable == equation._auxVariable ) { auxCoefficient = coefficient; @@ -1272,17 +1274,16 @@ void Tableau::addEquation( const Equation &equation ) ub -= coefficient * _upperBounds[variable]; } } - - if ( FloatUtils::isPositive( auxCoefficient ) ) - { - setLowerBound( equation._auxVariable, lb / auxCoefficient ); - setUpperBound( equation._auxVariable, ub / auxCoefficient ); - } - else - { - setLowerBound( equation._auxVariable, ub / auxCoefficient ); - setUpperBound( equation._auxVariable, lb / auxCoefficient ); - } + } + if ( FloatUtils::isPositive( auxCoefficient ) ) + { + setLowerBound( equation._auxVariable, lb / auxCoefficient ); + setUpperBound( equation._auxVariable, ub / auxCoefficient ); + } + else + { + setLowerBound( equation._auxVariable, ub / auxCoefficient ); + setUpperBound( equation._auxVariable, lb / auxCoefficient ); } // Populate the new row of A and b From a67260121f3ddb7fcd6367ad1b4dafcd7b318161 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Wed, 23 May 2018 15:08:03 +0300 Subject: [PATCH 14/14] oops --- src/engine/Tableau.cpp | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 1acd1c9b4..f8cd56ead 100644 --- a/src/engine/Tableau.cpp +++ b/src/engine/Tableau.cpp @@ -1331,17 +1331,6 @@ void Tableau::addEquation( const Equation &equation ) if ( FloatUtils::isZero( _basicAssignment[_m - 1] ) ) _basicAssignment[_m - 1] = 0.0; - // Refactorize the basis - try - { - _basisFactorization->refactorizeBasis(); - } - catch ( MalformedBasisException & ) - { - log( "addEquation failed - could not refactorize basis" ); - throw ReluplexError( ReluplexError::FAILURE_TO_ADD_NEW_EQUATION ); - } - // Notify about the new variable's assignment and compute its status notifyVariableValue( _basicIndexToVariable[_m - 1], _basicAssignment[_m - 1] ); computeBasicStatus( _m - 1 );