Skip to content

Commit

Permalink
better stats for pivoting (#33)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@guykatzz
guykatzz authored May 8, 2018
1 parent 8e15b7a commit 8b6f1ce
Show file tree
Hide file tree
Showing 7 changed files with 207 additions and 50 deletions.
162 changes: 113 additions & 49 deletions src/basis_factorization/ForrestTomlinFactorization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
, _Q( m )
, _invQ( m )
, _U( NULL )
, _invLP( NULL )
, _explicitBasisAvailable( false )
, _workMatrix( NULL )
, _workVector( NULL )
Expand All @@ -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 );
Expand Down Expand Up @@ -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()
Expand All @@ -94,7 +112,13 @@ ForrestTomlinFactorization::~ForrestTomlinFactorization()
_U = NULL;
}

if ( _workMatrix )
if ( _invLP )
{
delete[] _invLP;
_invLP = NULL;
}

if ( _workMatrix )
{
delete[] _workMatrix;
_workMatrix = NULL;
Expand Down Expand Up @@ -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<AlmostIdentityMatrix *>::iterator aIt;
for ( aIt = otherFTFactorization->_A.begin(); aIt != otherFTFactorization->_A.end(); ++aIt )
delete *aIt;
Expand Down Expand Up @@ -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<AlmostIdentityMatrix *>::iterator aIt;
for ( aIt = _A.begin(); aIt != _A.end(); ++aIt )
delete *aIt;
Expand Down Expand Up @@ -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
Expand All @@ -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;
}

Expand Down Expand Up @@ -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 ../.. "
Expand Down
8 changes: 8 additions & 0 deletions src/basis_factorization/ForrestTomlinFactorization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ class ForrestTomlinFactorization : public IBasisFactorization
const EtaMatrix **getU() const;
const List<LPElement *> *getLP() const;
const List<AlmostIdentityMatrix *> *getA() const;
const double *getInvLP() const;

void pushA( const AlmostIdentityMatrix &matrix );
void setQ( const PermutationMatrix &Q );
Expand All @@ -123,6 +124,13 @@ class ForrestTomlinFactorization : public IBasisFactorization
// A1 is the first element of the list
List<AlmostIdentityMatrix *>_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.
*/
Expand Down
11 changes: 11 additions & 0 deletions src/basis_factorization/tests/Test_CompareFactorizations.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
42 changes: 42 additions & 0 deletions src/basis_factorization/tests/Test_ForrestTomlinFactorization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
}

Expand Down
11 changes: 10 additions & 1 deletion src/engine/Statistics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ Statistics::Statistics()
, _numTableauPivots( 0 )
, _numTableauDegeneratePivots( 0 )
, _numTableauDegeneratePivotsByRequest( 0 )
, _timePivotsMicro( 0 )
, _numSimplexPivotSelectionsIgnoredForStability( 0 )
, _numSimplexUnstablePivots( 0 )
, _numTableauBoundHopping( 0 )
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -368,6 +372,11 @@ void Statistics::incNumTableauDegeneratePivotsByRequest()
++_numTableauDegeneratePivotsByRequest;
}

void Statistics::addTimePivots( unsigned long long time )
{
_timePivotsMicro += time;
}

void Statistics::incNumSimplexPivotSelectionsIgnoredForStability()
{
++_numSimplexPivotSelectionsIgnoredForStability;
Expand Down
4 changes: 4 additions & 0 deletions src/engine/Statistics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit 8b6f1ce

Please sign in to comment.