From f318b4f06416358b36d9c3460252ee1550794da7 Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 8 May 2018 11:11:10 +0300 Subject: [PATCH 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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 6/6] 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()