More accurate cost function computation (#75)

* 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

src/configuration/GlobalConfiguration.cpp

@@ -27,6 +27,7 @@ const double GlobalConfiguration::GAUSSIAN_ELIMINATION_PIVOT_SCALE_THRESHOLD = 0
 const unsigned GlobalConfiguration::MAX_SIMPLEX_PIVOT_SEARCH_ITERATIONS = 5;
 const unsigned GlobalConfiguration::CONSTRAINT_VIOLATION_THRESHOLD = 20;
 const unsigned GlobalConfiguration::BOUND_TIGHTING_ON_CONSTRAINT_MATRIX_FREQUENCY = 100;
+const double GlobalConfiguration::COST_FUNCTION_ERROR_THRESHOLD = 0.0000000001;
 
 const bool GlobalConfiguration::PREPROCESS_INPUT_QUERY = true;
 const bool GlobalConfiguration::PREPROCESSOR_ELIMINATE_VARIABLES = true;

src/configuration/GlobalConfiguration.h

@@ -72,6 +72,9 @@ class GlobalConfiguration
     // How often should we perform full bound tightening, on the entire contraints matrix A.
     static const unsigned BOUND_TIGHTING_ON_CONSTRAINT_MATRIX_FREQUENCY;
 
+    // If the cost function error exceeds this threshold, it is recomputed
+    static const double COST_FUNCTION_ERROR_THRESHOLD;
+
     // How often should projected steepest edge reset the reference space?
     static const unsigned PSE_ITERATIONS_BEFORE_RESET;
 

src/engine/CostFunctionManager.cpp

@@ -219,10 +219,12 @@ const double *CostFunctionManager::getCostFunction() const
     return _costFunction;
 }
 
-void CostFunctionManager::updateCostFunctionForPivot( unsigned enteringVariableIndex,
-                                                      unsigned leavingVariableIndex,
-                                                      double pivotElement,
-                                                      const TableauRow *pivotRow )
+double CostFunctionManager::updateCostFunctionForPivot( unsigned enteringVariableIndex,
+                                                        unsigned leavingVariableIndex,
+                                                        double pivotElement,
+                                                        const TableauRow *pivotRow,
+                                                        const double *changeColumn
+                                                        )
 {
     /*
       This method is invoked when the non-basic _enteringVariable and
@@ -247,36 +249,39 @@ void CostFunctionManager::updateCostFunctionForPivot( unsigned enteringVariableI
     ASSERT( _tableau->getM() == _m );
     ASSERT( _tableau->getN() == _n );
 
+    /*
+      The current reduced cost of the entering variable is stored in
+      _costFunction, but since we have the change column we can compute a
+      more accurate version from scratch
+    */
+    double enteringVariableCost = 0;
+    for ( unsigned i = 0; i < _m; ++i )
+        enteringVariableCost -= _basicCosts[i] * changeColumn[i];
+
+    double normalizedError =
+        FloatUtils::abs( enteringVariableCost - _costFunction[enteringVariableIndex] ) /
+        ( FloatUtils::abs( enteringVariableCost ) + 1.0 );
+
     // Update the cost of the new non-basic
-    _costFunction[enteringVariableIndex] /= pivotElement;
+    _costFunction[enteringVariableIndex] = enteringVariableCost / pivotElement;
 
     for ( unsigned i = 0; i < _n - _m; ++i )
     {
         if ( i != enteringVariableIndex )
             _costFunction[i] -= (*pivotRow)[i] * _costFunction[enteringVariableIndex];
     }
 
-    unsigned leavingVariableStatus = _tableau->getBasicStatusByIndex( leavingVariableIndex );
+    /*
+      The leaving variable might have contributed to the cost function, but it will
+      soon be made within bounds. So, we adjust the reduced costs accordingly.
+    */
+    _costFunction[enteringVariableIndex] -= _basicCosts[leavingVariableIndex];
 
-    // Update the basic cost for the leaving variable, which may have changed
-    // since we last computed it
-    switch ( leavingVariableStatus )
-    {
-    case ITableau::ABOVE_UB:
-        _basicCosts[leavingVariableIndex] = 1;
-        break;
-    case ITableau::BELOW_LB:
-        _basicCosts[leavingVariableIndex] = -1;
-        break;
-    default:
-        _basicCosts[leavingVariableIndex] = 0;
-        break;
-    }
+    // The entering varibale is non-basic, so it is within bounds.
+    _basicCosts[leavingVariableIndex] = 0;
 
-    // If the leaving variable was previously out of bounds, this is no longer
-    // the case. Adjust the non-basic cost.
-    _costFunction[enteringVariableIndex] -= _basicCosts[leavingVariableIndex];
     _costFunctionStatus = ICostFunctionManager::COST_FUNCTION_UPDATED;
+    return normalizedError;
 }
 
 bool CostFunctionManager::costFunctionInvalid() const

src/engine/CostFunctionManager.h

@@ -56,10 +56,11 @@ class CostFunctionManager : public ICostFunctionManager
       Update the cost fucntion just before a coming pivot step, to avoid having to compute
       it from scratch afterwards.
     */
-    void updateCostFunctionForPivot( unsigned enteringVariableIndex,
-                                     unsigned leavingVariableIndex,
-                                     double pivotElement,
-                                     const TableauRow *pivotRow );
+    double updateCostFunctionForPivot( unsigned enteringVariableIndex,
+                                       unsigned leavingVariableIndex,
+                                       double pivotElement,
+                                       const TableauRow *pivotRow,
+                                       const double *changeColumn );
 
     /*
       For debugging purposes: dump the cost function.

src/engine/ICostFunctionManager.h

@@ -35,10 +35,11 @@ class ICostFunctionManager
     virtual const double *getCostFunction() const = 0;
     virtual void dumpCostFunction() const = 0;
     virtual void setCostFunctionStatus( ICostFunctionManager::CostFunctionStatus status ) = 0;
-    virtual void updateCostFunctionForPivot( unsigned enteringVariableIndex,
-                                             unsigned leavingVariableIndex,
-                                             double pivotElement,
-                                             const TableauRow *pivotRow ) = 0;
+    virtual double updateCostFunctionForPivot( unsigned enteringVariableIndex,
+                                               unsigned leavingVariableIndex,
+                                               double pivotElement,
+                                               const TableauRow *pivotRow,
+                                               const double *changeColumn ) = 0;
 
     virtual bool costFunctionInvalid() const = 0;
     virtual bool costFunctionJustComputed() const = 0;

src/engine/Tableau.cpp

@@ -1977,10 +1977,14 @@ void Tableau::updateCostFunctionForPivot()
         return;
 
     double pivotElement = -_changeColumn[_leavingVariable];
-    _costFunctionManager->updateCostFunctionForPivot( _enteringVariable,
-                                                      _leavingVariable,
-                                                      pivotElement,
-                                                      _pivotRow );
+    double normalizedError = _costFunctionManager->updateCostFunctionForPivot( _enteringVariable,
+                                                                               _leavingVariable,
+                                                                               pivotElement,
+                                                                               _pivotRow,
+                                                                               _changeColumn );
+
+    if ( FloatUtils::gt( normalizedError, GlobalConfiguration::COST_FUNCTION_ERROR_THRESHOLD ) )
+        _costFunctionManager->invalidateCostFunction();
 }
 
 ITableau::BasicAssignmentStatus Tableau::getBasicAssignmentStatus() const

src/engine/tests/MockCostFunctionManager.h

@@ -92,11 +92,14 @@ class MockCostFunctionManager : public ICostFunctionManager
     {
     }
 
-    void updateCostFunctionForPivot( unsigned /* enteringVariableIndex */,
-                                     unsigned /* leavingVariableIndex */,
-                                     double /* pivotElement */,
-                                     const TableauRow */* pivotRow */ )
+    double updateCostFunctionForPivot( unsigned /* enteringVariableIndex */,
+                                       unsigned /* leavingVariableIndex */,
+                                       double /* pivotElement */,
+                                       const TableauRow */* pivotRow */,
+                                       const double */* changeColumn */
+                                       )
     {
+        return 0;
     }
 
     bool costFunctionInvalid() const