diff --git a/src/basis_factorization/GaussianEliminator.cpp b/src/basis_factorization/GaussianEliminator.cpp index f569181b5..57ddf017a 100644 --- a/src/basis_factorization/GaussianEliminator.cpp +++ b/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]. diff --git a/src/engine/Engine.cpp b/src/engine/Engine.cpp index 9367bd62d..e1a044748 100644 --- a/src/engine/Engine.cpp +++ b/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 { diff --git a/src/engine/PrecisionRestorer.cpp b/src/engine/PrecisionRestorer.cpp index 422f94fb4..0cb52b714 100644 --- a/src/engine/PrecisionRestorer.cpp +++ b/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 currentBasics = tableau.getBasicVariables(); - if ( restoreBasics == RESTORE_BASICS ) + if ( dimensionsRestored && restoreBasics == RESTORE_BASICS ) { List 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 ) diff --git a/src/engine/ProjectedSteepestEdge.cpp b/src/engine/ProjectedSteepestEdge.cpp index f4a8302e7..9e1d7f6de 100644 --- a/src/engine/ProjectedSteepestEdge.cpp +++ b/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 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 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 ); diff --git a/src/engine/ProjectedSteepestEdge.h b/src/engine/ProjectedSteepestEdge.h index 46ba68885..c52d2a6d5 100644 --- a/src/engine/ProjectedSteepestEdge.h +++ b/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. diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp index 23580d691..369fd8218 100644 --- a/src/engine/Tableau.cpp +++ b/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() ) {