From 8b6f1ce1f0dd87e892580920bb875b23c3a6216f Mon Sep 17 00:00:00 2001 From: Guy Katz Date: Tue, 8 May 2018 14:00:26 +0300 Subject: [PATCH] better stats for pivoting (#33) * 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 * pivoting stats * bug fix in store/restore --- .../ForrestTomlinFactorization.cpp | 162 ++++++++++++------ .../ForrestTomlinFactorization.h | 8 + .../tests/Test_CompareFactorizations.h | 11 ++ .../tests/Test_ForrestTomlinFactorization.h | 42 +++++ src/engine/Statistics.cpp | 11 +- src/engine/Statistics.h | 4 + src/engine/Tableau.cpp | 19 ++ 7 files changed, 207 insertions(+), 50 deletions(-) diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp index ddf760131..3c0d98635 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; @@ -509,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; @@ -546,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; @@ -663,6 +691,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 +764,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 +991,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 ); } 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 )