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: +// 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] );