Bug fixes

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() )
     {