Bug fixes (#77)

* initial work on sparse matrix representation

* store/restore functionality

* addLastRow functionality

* getRow and getColumn

* column-merging functionality

* added an interface class

* introducing also sparse vectors

* added addLastColumn functionality

* another unittest

* get sparse columns/matrices in dense form

* WIP on storing the constraint matrix inside the tableau in sparse form

* more WIP, fixed a few bugs, still have a couple of failing tests

* fixed some test issues

* initialization

* resize _a along with the rest

* sparse lu factors

* store also the transposed versions of F and V

* starting work on the sparse GE

* some work on changing the values within an existing sparse
representation, needed for sparse factorization.
still WIP

* refactoring and new functionality for CSRMatrix: any kind of
insertions and deletions

* support for empty initialization and counting elements

* sparse GE is now working. minor bug fixes elsewhere

* compute Ft and Vt as part of the G-elimination process

* tests

* basis oracles can return also sparse columns, not just dense

* sparse LU factorization class

* switch to using the sparse factorization in the engine/tableau

* bug fix in mergeColumns, and nicer printing

* bug fix

* bug fix

* bug fix: merging columns does not delete the actual column, just
leaves it empty

* configuration changes

* optimization: since the sparse columns of A are needed all the time,
just compute them once-and-for-all

* a more efficient implementation of sparse vectors

* comments and unit tests

* cleanup

* keep _A in dense column-major format, too, instead of repeatedly
invoking toDense() to gets its columns

* bad deletes

* bug fix in test

* bug fixes: relu constraint propagation, and the handling of merged
variables in the preprocessor

* new test

* compute Ft incrementally, use it whenever sparse columns are required

* did the todo

* valgrind fixes

* debugging

* un-initialized memory

* cleanup

* fixing an issue with cost function degradation

* reinstating the anti-looping trick

* debug prints

* verify invariants

* debug info for crash

* fixing a numerical stability issue
new assertions

* cleanup

* cleanup

* removing the code with the failing assertion

* weaken a couple of too-storng assertions

* bug fix and some fine-tuning of PSE

* bug fix in LU factorization
optimizations to PSE

src/basis_factorization/GaussianEliminator.cpp

@@ -318,6 +318,18 @@ void GaussianEliminator::eliminate()
 {
     unsigned fColumnIndex = _luFactors->_P._columnOrdering[_eliminationStep];
 
+    /*
+      The pivot row is not eliminated per se, but it is excluded
+      from the active submatrix, so we adjust the element counters
+    */
+    _numURowElements[_eliminationStep] = 0;
+    for ( unsigned uColumn = _eliminationStep; uColumn < _m; ++uColumn )
+    {
+        unsigned vColumn = _luFactors->_Q._rowOrdering[uColumn];
+        if ( !FloatUtils::isZero( _luFactors->_V[_vPivotRow*_m + vColumn] ) )
+            --_numUColumnElements[uColumn];
+    }
+
     /*
       Eliminate all entries below the pivot element U[k,k]
       We know that V[_pivotRow, _pivotColumn] = U[k,k].

src/engine/Engine.cpp

@@ -502,7 +502,7 @@ void Engine::fixViolatedPlConstraintIfPossible()
 
     _activeEntryStrategy->prePivotHook( _tableau, false );
     _tableau->performDegeneratePivot();
-    _activeEntryStrategy->prePivotHook( _tableau, false );
+    _activeEntryStrategy->postPivotHook( _tableau, false );
 
     ASSERT( !_tableau->isBasic( fix._variable ) );
     _tableau->setNonBasicAssignment( fix._variable, fix._value, true );
@@ -817,7 +817,7 @@ void Engine::applySplit( const PiecewiseLinearCaseSplit &split )
 
                 _activeEntryStrategy->prePivotHook( _tableau, false );
                 _tableau->performDegeneratePivot();
-                _activeEntryStrategy->prePivotHook( _tableau, false );
+                _activeEntryStrategy->postPivotHook( _tableau, false );
             }
 
             if ( _tableau->isBasic( x2 ) )
@@ -856,13 +856,16 @@ void Engine::applySplit( const PiecewiseLinearCaseSplit &split )
 
                 _activeEntryStrategy->prePivotHook( _tableau, false );
                 _tableau->performDegeneratePivot();
-                _activeEntryStrategy->prePivotHook( _tableau, false );
+                _activeEntryStrategy->postPivotHook( _tableau, false );
 
             }
 
             // Both variables are now non-basic, so we can merge their columns
             _tableau->mergeColumns( x1, x2 );
             DEBUG( _tableau->verifyInvariants() );
+
+            // Reset the entry strategy
+            _activeEntryStrategy->initialize( _tableau );
         }
         else
         {

src/engine/PrecisionRestorer.cpp

@@ -56,13 +56,17 @@ void PrecisionRestorer::restorePrecision( IEngine &engine,
 
         // At this point, the tableau has the appropriate dimensions. Restore the variable bounds
         // and basic variables.
+        // Note that if column merging is enabled, the dimensions may not be precisely those before
+        // the resotration, because merging sometimes fails - in which case an equation is added. If
+        // we fail to restore the dimensions, we cannot restore the basics.
 
-        ASSERT( tableau.getN() == targetN );
-        ASSERT( tableau.getM() == targetM );
+        bool dimensionsRestored = ( tableau.getN() == targetN ) && ( tableau.getM() == targetM );
+
+        ASSERT( dimensionsRestored || GlobalConfiguration::USE_COLUMN_MERGING_EQUATIONS );
 
         Set<unsigned> currentBasics = tableau.getBasicVariables();
 
-        if ( restoreBasics == RESTORE_BASICS )
+        if ( dimensionsRestored && restoreBasics == RESTORE_BASICS )
         {
             List<unsigned> shouldBeBasicList;
             for ( const auto &basic : shouldBeBasic )
@@ -115,8 +119,8 @@ void PrecisionRestorer::restorePrecision( IEngine &engine,
 
         DEBUG({
                 // Same dimensions
-                ASSERT( tableau.getN() == targetN );
-                ASSERT( tableau.getM() == targetM );
+                ASSERT( GlobalConfiguration::USE_COLUMN_MERGING_EQUATIONS || tableau.getN() == targetN );
+                ASSERT( GlobalConfiguration::USE_COLUMN_MERGING_EQUATIONS || tableau.getM() == targetM );
 
                 // Constraints should be in the same state before and after restoration
                 for ( const auto &pair : targetEngineState._plConstraintToState )

src/engine/ProjectedSteepestEdge.cpp

@@ -101,6 +101,7 @@ void ProjectedSteepestEdgeRule::resetReferenceSpace( const ITableau &tableau )
     }
 
     _iterationsUntilReset = GlobalConfiguration::PSE_ITERATIONS_BEFORE_RESET;
+    _errorInGamma = 0.0;
 
     if ( _statistics )
         _statistics->pseIncNumResetReferenceSpace();
@@ -145,7 +146,7 @@ bool ProjectedSteepestEdgeRule::select( ITableau &tableau, const Set<unsigned> &
     unsigned bestCandidate = *it;
     double gammaValue = _gamma[*it];
     double bestValue =
-        !FloatUtils::isPositive( gammaValue ) ? 0 : ( costFunction[*it] * costFunction[*it] ) / _gamma[*it];
+        !FloatUtils::isPositive( gammaValue ) ? 0 : ( costFunction[*it] * costFunction[*it] ) / gammaValue;
 
     ++it;
 
@@ -154,7 +155,7 @@ bool ProjectedSteepestEdgeRule::select( ITableau &tableau, const Set<unsigned> &
         unsigned contender = *it;
         gammaValue = _gamma[*it];
         double contenderValue =
-            !FloatUtils::isPositive( gammaValue ) ? 0 : ( costFunction[*it] * costFunction[*it] ) / _gamma[*it];
+            !FloatUtils::isPositive( gammaValue ) ? 0 : ( costFunction[*it] * costFunction[*it] ) / gammaValue;
 
         if ( FloatUtils::gt( contenderValue, bestValue ) )
         {
@@ -184,12 +185,14 @@ void ProjectedSteepestEdgeRule::prePivotHook( const ITableau &tableau, bool fake
     }
 
     // When this hook is called, the entering and leaving variables have
-    // already been determined. These are the actual varaibles, not the indices.
+    // already been determined.
     unsigned entering = tableau.getEnteringVariable();
     unsigned enteringIndex = tableau.variableToIndex( entering );
     unsigned leaving = tableau.getLeavingVariable();
     unsigned leavingIndex = tableau.variableToIndex( leaving );
 
+    ASSERT( entering != leaving );
+
     const double *changeColumn = tableau.getChangeColumn();
     const TableauRow &pivotRow = *tableau.getPivotRow();
 
@@ -222,7 +225,7 @@ void ProjectedSteepestEdgeRule::prePivotHook( const ITableau &tableau, bool fake
         if ( i == enteringIndex )
             continue;
 
-        if ( FloatUtils::isZero( pivotRow[i] ) )
+        if ( FloatUtils::isZero( pivotRow[i], 1e-9 ) )
             continue;
 
         r = pivotRow[i] / -changeColumn[leavingIndex];
@@ -279,7 +282,7 @@ void ProjectedSteepestEdgeRule::postPivotHook( const ITableau &tableau, bool fak
 
     // If the iteration limit has been exhausted, reset the reference space
     --_iterationsUntilReset;
-    if ( _iterationsUntilReset == 0 )
+    if ( _iterationsUntilReset <= 0 )
     {
         log( "PostPivotHook reseting ref space (iterations)" );
         resetReferenceSpace( tableau );

src/engine/ProjectedSteepestEdge.h

@@ -80,7 +80,7 @@ class ProjectedSteepestEdgeRule : public IProjectedSteepestEdgeRule
     /*
       Remaining iterations before resetting the reference space.
     */
-    unsigned _iterationsUntilReset;
+    int _iterationsUntilReset;
 
     /*
       The error in gamma compuated in the previous iteration.

src/engine/Tableau.cpp

@@ -1861,37 +1861,37 @@ void Tableau::updateAssignmentForPivot()
 
     _basicAssignmentStatus = ITableau::BASIC_ASSIGNMENT_UPDATED;
 
-    // If the change ratio is 0, just maintain the current assignment
-    if ( FloatUtils::isZero( _changeRatio ) )
-    {
-        ASSERT( !performingFakePivot() );
-
-        DEBUG({
-                // This should only happen when the basic variable is pressed against
-                // one of its bounds
-                if ( !( _basicStatus[_leavingVariable] == Tableau::AT_UB ||
-                        _basicStatus[_leavingVariable] == Tableau::AT_LB ||
-                        _basicStatus[_leavingVariable] == Tableau::BETWEEN
-                        ) )
-                {
-                    printf( "Assertion violation!\n" );
-                    printf( "Basic (leaving) variable is: %u\n", _basicIndexToVariable[_leavingVariable] );
-                    printf( "Basic assignment: %.10lf. Bounds: [%.10lf, %.10lf]\n",
-                            _basicAssignment[_leavingVariable],
-                            _lowerBounds[_basicIndexToVariable[_leavingVariable]],
-                            _upperBounds[_basicIndexToVariable[_leavingVariable]] );
-                    printf( "Basic status: %u\n", _basicStatus[_leavingVariable] );
-                    printf( "leavingVariableIncreases = %s", _leavingVariableIncreases ? "yes" : "no" );
-                    exit( 1 );
-                }
-            });
-
-        double basicAssignment = _basicAssignment[_leavingVariable];
-        double nonBasicAssignment = _nonBasicAssignment[_enteringVariable];
-        _basicAssignment[_leavingVariable] = nonBasicAssignment;
-        _nonBasicAssignment[_enteringVariable] = basicAssignment;
-        return;
-    }
+    // // If the change ratio is 0, just maintain the current assignment
+    // if ( FloatUtils::isZero( _changeRatio ) )
+    // {
+    //     ASSERT( !performingFakePivot() );
+
+    //     DEBUG({
+    //             // This should only happen when the basic variable is pressed against
+    //             // one of its bounds
+    //             if ( !( _basicStatus[_leavingVariable] == Tableau::AT_UB ||
+    //                     _basicStatus[_leavingVariable] == Tableau::AT_LB ||
+    //                     _basicStatus[_leavingVariable] == Tableau::BETWEEN
+    //                     ) )
+    //             {
+    //                 printf( "Assertion violation!\n" );
+    //                 printf( "Basic (leaving) variable is: %u\n", _basicIndexToVariable[_leavingVariable] );
+    //                 printf( "Basic assignment: %.10lf. Bounds: [%.10lf, %.10lf]\n",
+    //                         _basicAssignment[_leavingVariable],
+    //                         _lowerBounds[_basicIndexToVariable[_leavingVariable]],
+    //                         _upperBounds[_basicIndexToVariable[_leavingVariable]] );
+    //                 printf( "Basic status: %u\n", _basicStatus[_leavingVariable] );
+    //                 printf( "leavingVariableIncreases = %s", _leavingVariableIncreases ? "yes" : "no" );
+    //                 exit( 1 );
+    //             }
+    //         });
+
+    //     double basicAssignment = _basicAssignment[_leavingVariable];
+    //     double nonBasicAssignment = _nonBasicAssignment[_enteringVariable];
+    //     _basicAssignment[_leavingVariable] = nonBasicAssignment;
+    //     _nonBasicAssignment[_enteringVariable] = basicAssignment;
+    //     return;
+    // }
 
     if ( performingFakePivot() )
     {