cleaned up unused code

1 files changed, 62 insertions, 312 deletions

diff --git a/ripser.cpp b/ripser.cpp
index 8a5b3ec..c4310c5 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -7,7 +7,6 @@
 #include <queue>
 #include <unordered_map>
 #include "prettyprint.hpp"
-//#include <boost/numeric/ublas/symmetric.hpp>
 
 typedef float value_t;
 typedef long index_t;
@@ -53,15 +52,12 @@ public:
     }
     
     inline index_t operator()(index_t n, index_t k) const {
-//        std::cout << "B(" << n << "," << k << ")\n";
         assert(n <= n_max);
         assert(k <= k_max);
         return B[n][k];
     }
 };
 
-
-
 template<typename OutputIterator>
 OutputIterator get_simplex_vertices( index_t idx, const index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out  )
 {
@@ -126,17 +122,10 @@ private:
     
 public:
     rips_filtration_comparator(
-    const DistanceMatrix& _dist,
-    const index_t _dim,
-    const binomial_coeff_table& _binomial_coeff
-//    ,
-//    bool _reverse = true
-    ):
-    dist(_dist), dim(_dim), vertices(_dim + 1),
-    binomial_coeff(_binomial_coeff)
-//    ,
-//    reverse(_reverse)
-    {};
+        const DistanceMatrix& _dist,
+        const index_t _dim,
+        const binomial_coeff_table& _binomial_coeff
+    ): dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff) {};
     
     inline dist_t diameter(const index_t index) const {
         dist_t diam = 0;
@@ -155,7 +144,6 @@ public:
         assert(b < binomial_coeff(dist.size(), dim + 1));
         
         dist_t a_diam = 0, b_diam = 0;
-
     
         b_diam = diameter(b);
         
@@ -169,99 +157,8 @@ public:
        
         return (a_diam == b_diam) && (a > b);
     }
-    
 };
 
-template <typename DistanceMatrix, index_t dim>
-class rips_filtration_comparator_dim;
-
-template <typename DistanceMatrix>
-class rips_filtration_comparator_dim <DistanceMatrix, 2> {
-public:
-    const DistanceMatrix& dist;
-    
-private:
-    mutable std::vector<index_t> vertices;
-    
-    typedef decltype(dist(0,0)) dist_t;
-    
-    bool reverse;
-    
-    const binomial_coeff_table& binomial_coeff;
-    
-public:
-    rips_filtration_comparator_dim(
-    const DistanceMatrix& _dist,
-    const binomial_coeff_table& _binomial_coeff
-    ):
-    dist(_dist), vertices(3),
-    binomial_coeff(_binomial_coeff)
-    {};
-    
-    inline dist_t diameter(const index_t index) const {
-        dist_t diam = 0;
-        get_simplex_vertices(index, 2, dist.size(), binomial_coeff, vertices.begin() );
-        
-        diam = std::max(diam, dist(vertices[0], vertices[1]));
-        diam = std::max(diam, dist(vertices[0], vertices[2]));
-        diam = std::max(diam, dist(vertices[0], vertices[3]));
-        diam = std::max(diam, dist(vertices[1], vertices[2]));
-        diam = std::max(diam, dist(vertices[1], vertices[3]));
-        diam = std::max(diam, dist(vertices[2], vertices[3]));
-
-        return diam;
-    }
-
-    inline bool operator()(const index_t a, const index_t b) const
-    {
-        assert(a < binomial_coeff(dist.size(), 3));
-        assert(b < binomial_coeff(dist.size(), 3));
-        
-        dist_t a_diam = 0, b_diam = 0;
-
-    
-        b_diam = diameter(b);
-        
-        get_simplex_vertices(a, 2, dist.size(), binomial_coeff, vertices.begin() );
-        if ((a_diam = dist(vertices[0], vertices[1])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        if ((a_diam = dist(vertices[0], vertices[2])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        if ((a_diam = dist(vertices[0], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        if ((a_diam = dist(vertices[1], vertices[2])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        if ((a_diam = dist(vertices[1], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        if ((a_diam = dist(vertices[2], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b);
-        
-        return false;
-    }
-    
-};
-
-
-template<typename OutputIterator>
-void inline get_simplex_coboundary( index_t idx, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator coboundary )
-{
-
-    --n;
-    
-    index_t modified_idx = idx;
-    
-    for( index_t k = dim + 1; k != -1 && n != -1; --k ) {
-        while( binomial_coeff( n , k ) > idx ) {
-            *coboundary++ = modified_idx + binomial_coeff( n , k + 1 );
-            if (n==0) break;
-            --n;
-        }
-        
-        idx -= binomial_coeff( n , k );
-            
-        modified_idx -= binomial_coeff( n , k );
-        modified_idx += binomial_coeff( n , k + 1 );
-        
-        --n;
-    }
- 
-    return;
-}
-
 class simplex_coboundary_enumerator {
 private:
     index_t idx, modified_idx, dim, n, k;
@@ -270,14 +167,9 @@ private:
     
 public:
     inline simplex_coboundary_enumerator (
-        index_t _idx,
-        index_t _dim,
-        index_t _n,
+        index_t _idx, index_t _dim, index_t _n,
         const binomial_coeff_table& _binomial_coeff):
-        idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff)
-    {
-        
-    }
+        idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff) {}
     
     inline bool has_next()
     {
@@ -309,8 +201,7 @@ public:
 
 template <typename DistanceMatrix>
 std::vector<value_t> get_diameters (
-    const DistanceMatrix& dist,
-    const index_t dim,
+    const DistanceMatrix& dist, const index_t dim,
     const std::vector<value_t>& previous_diameters,
     const binomial_coeff_table& binomial_coeff
     )
@@ -339,11 +230,6 @@ std::vector<value_t> get_diameters (
     std::cout << "\033[K";
     #endif
     
-//    rips_filtration_comparator<decltype(dist)> comp(dist, dim, binomial_coeff);
-//    for (index_t simplex = 0; simplex < diameters.size(); ++simplex) {
-//        assert(diameters[simplex] == comp.diameter(simplex));
-//    }
-
     return diameters;
 }
 
@@ -359,7 +245,6 @@ public:
     
     typedef value_t dist_t;
     
-    
     const binomial_coeff_table& binomial_coeff;
     
 public:
@@ -386,7 +271,6 @@ public:
 
         return ((a_diam > b_diam) || ((a_diam == b_diam) && (a > b)));
     }
-    
 };
 
 
@@ -401,7 +285,6 @@ public:
     inline size_t size() const {
         return distances.size();
     }
-
 };
 
 
@@ -439,7 +322,6 @@ public:
     inline size_t size() const {
         return n;
     }
-
 };
 
 
@@ -474,8 +356,6 @@ public:
         init();
     }
     
-//    template <typename DistanceMatrix>
-//    compressed_lower_distance_matrix(const DistanceMatrix& mat) : compressed_lower_distance_matrix(mat.size()) {
     compressed_lower_distance_matrix(const compressed_upper_distance_matrix& mat) {
         n = mat.size();
         
@@ -488,7 +368,6 @@ public:
         }
     }
     
-    
     inline value_t operator()(const index_t i, const index_t j) const {
         if (i > j)
             return rows[i][j];
@@ -501,44 +380,6 @@ public:
     inline size_t size() const {
         return n;
     }
-
-};
-
-template <typename ValueType>
-class compressed_sparse_matrix  {
-public:
-    std::vector<size_t> bounds;
-    std::vector<ValueType> entries;
-    
-    
-    inline size_t size() const {
-        return bounds.size();
-    }
-    
-    inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const {
-        assert(index < size());
-        return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1];
-    }
-
-    inline typename std::vector<ValueType>::const_iterator cend(size_t index) const {
-        assert(index < size());
-        return entries.cbegin() + bounds[index];
-    }
-    
-    template <typename Iterator>
-    inline void append(Iterator begin, Iterator end) {
-        for (Iterator it = begin; it != end; ++it) {
-            entries.push_back(*it);
-        }
-        bounds.push_back(entries.size());
-    }
-
-    
-    template <typename Collection>
-    inline void append(const Collection collection) {
-        append(collection.cbegin(), collection.cend());
-    }
-
 };
 
 template <typename Heap>
@@ -571,20 +412,6 @@ inline index_t get_pivot(Heap& column)
     return max_element;
 }
 
-template <typename Heap>
-inline std::vector<index_t> move_to_column_vector(Heap& column)
-{
-    std::vector<index_t> temp_col;
-    index_t pivot = pop_pivot( column );
-    while( pivot != -1 ) {
-        temp_col.push_back( pivot );
-        pivot = pop_pivot( column );
-    }
-    return temp_col;
-}
-
-
-
 void print_help_and_exit()
 {
     std::cerr << "Usage: " << "ripser " << "[options] input_filename output_filename" << std::endl;
@@ -608,85 +435,73 @@ void compute_pairs(
     const binomial_coeff_table& binomial_coeff
 ) {
     
-        std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
+    std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
+    
+    for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
+        index_t index = columns_to_reduce[i];
+        
+        #ifdef ASSEMBLE_REDUCTION_COLUMN
+        std::priority_queue<index_t, std::vector<index_t>, decltype(comp) > reduction_column(comp);
+        #endif
+        
+        std::priority_queue<index_t, std::vector<index_t>, decltype(comp) > working_coboundary(comp);
+        
+        #ifdef INDICATE_PROGRESS
+        std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size()
+        << " (diameter " << comp_prev.diameter(index) << ")"
+        << std::flush << "\r";
+        #endif
+        
+        index_t pivot, column = index;
+        
+        std::vector<index_t> coboundary;
+
+        do {
         
-        for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
-            index_t index = columns_to_reduce[i];
-            
             #ifdef ASSEMBLE_REDUCTION_COLUMN
-            std::priority_queue<index_t, std::vector<index_t>, decltype(comp) > reduction_column(comp);
-            #endif
-            
-            std::priority_queue<index_t, std::vector<index_t>, decltype(comp) > working_coboundary(comp);
-            
-            #ifdef INDICATE_PROGRESS
-            std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size()
-            << " (diameter " << comp_prev.diameter(index) << ")"
-            << std::flush << "\r";
+            reduction_column.push( column );
             #endif
-            
-            index_t pivot, column = index;
-            
-            std::vector<index_t> coboundary;
 
-            do {
+            simplex_coboundary_enumerator coboundaries(column, dim, n, binomial_coeff);
+            while (coboundaries.has_next()) {
+                index_t coface = coboundaries.next();
+                if (comp.diameter(coface) <= threshold)
+                    working_coboundary.push(coface);
+            }
+                           
+            pivot = get_pivot(working_coboundary);
             
-                #ifdef ASSEMBLE_REDUCTION_COLUMN
-                reduction_column.push( column );
-                #endif
-
-                simplex_coboundary_enumerator coboundaries(column, dim, n, binomial_coeff);
-                while (coboundaries.has_next()) {
-                    index_t coface = coboundaries.next();
-                    if (comp.diameter(coface) <= threshold)
-                        working_coboundary.push(coface);
-                }
-                               
-                pivot = get_pivot(working_coboundary);
+            if (pivot != -1) {
+                auto pair = pivot_column_index.find(pivot);
                 
-
-                //add boundary column at index birth to working_coboundary
-                
-                //since the boundary column is not stored explicitly,
-                //add the boundary of the column at index birth in the reduction matrix instead
-                
-                //space-efficient variant: add just the boundary column at index birth instead
-                //this avoids having to store the reduction matrix
-                
-                if (pivot != -1) {
-                    auto pair = pivot_column_index.find(pivot);
-                    
-                    if (pair == pivot_column_index.end()) {
-                        pivot_column_index.insert(std::make_pair(pivot, index));
-                        
-                        #ifdef PRINT_PERSISTENCE_PAIRS
-                        value_t birth = comp_prev.diameter(index), death = comp.diameter(pivot);
-                        if (birth != death)
-                            std::cout << "\033[K" << " [" << birth << "," << death << ")" << std::endl << std::flush;
-                        #endif
-
-                        break;
-                    }
+                if (pair == pivot_column_index.end()) {
+                    pivot_column_index.insert(std::make_pair(pivot, index));
                     
-                    column = pair->second;
-                }
+                    #ifdef PRINT_PERSISTENCE_PAIRS
+                    value_t birth = comp_prev.diameter(index), death = comp.diameter(pivot);
+                    if (birth != death)
+                        std::cout << "\033[K" << " [" << birth << "," << death << ")" << std::endl << std::flush;
+                    #endif
 
-            } while ( pivot != -1 );
-            
-            #ifdef PRINT_PERSISTENCE_PAIRS
-            if ( pivot == -1 ) {
-                value_t birth = comp_prev.diameter(index);
-                std::cout << "\033[K" << " [" << birth << ", )" << std::endl << std::flush;
+                    break;
+                }
+                
+                column = pair->second;
             }
-            #endif
 
+        } while ( pivot != -1 );
+        
+        #ifdef PRINT_PERSISTENCE_PAIRS
+        if ( pivot == -1 ) {
+            value_t birth = comp_prev.diameter(index);
+            std::cout << "\033[K" << " [" << birth << ", )" << std::endl << std::flush;
         }
-    
-        std::cout << "\033[K";
-
-}
+        #endif
 
+    }
 
+    std::cout << "\033[K";
+}
 
 
 int main( int argc, char** argv ) {
@@ -740,7 +555,6 @@ int main( int argc, char** argv ) {
         exit(-1);
     }
 
-   
     int64_t file_type;
     input_stream.read( (char*)&file_type, sizeof( int64_t ) );
     if( file_type != 7 ) {
@@ -751,18 +565,12 @@ int main( int argc, char** argv ) {
     int64_t n;
     input_stream.read( (char*)&n, sizeof( int64_t ) );
     
-
-    //std::vector<value_t> distances =
-
     distance_matrix dist;
     dist.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n));
     
     for( int i = 0; i < n; ++i ) {
         input_stream.read( (char*)&dist.distances[i][0], n * sizeof(int64_t) );
     }
-    
-    //std::cout << dist.distances << std::endl;
-    
     #endif
 
 
@@ -773,8 +581,6 @@ int main( int argc, char** argv ) {
     while(std::getline(input_stream, value_string, ','))
         distances.push_back(std::stod(value_string));
     
-//    std::cout << "distances: " << distances << std::endl;
-    
     compressed_upper_distance_matrix dist_upper(std::move(distances));
     
     index_t n = dist_upper.size();
@@ -782,44 +588,15 @@ int main( int argc, char** argv ) {
     #endif
     
     std::cout << "distance matrix with " << n << " points" << std::endl;
-    
-//    std::vector<std::vector<value_t>> distance_matrix_full(n, std::vector<value_t>(n));
-//    for (index_t i = 0; i < n; ++i)
-//        for (index_t j = 0; j < n; ++j)
-//            distance_matrix_full[i][j] = dist_upper(i, j);
-//    std::cout << distance_matrix_full << std::endl;
-
-
-
+ 
     compressed_lower_distance_matrix dist(dist_upper);
 
     std::cout << "distance matrix transformed to lower triangular form" << std::endl;
 
-//    for (index_t i = 0; i < n; ++i)
-//        for (index_t j = 0; j < n; ++j)
-//            distance_matrix_full[i][j] = dist(i, j);
-//    std::cout << distance_matrix_full << std::endl;
-    
-    
-//    std::cout << "distances (lower): " << dist.distances << std::endl;
-
-//    return 0;
-    
     assert(dim_max < n - 2);
     
     binomial_coeff_table binomial_coeff(n, dim_max + 2);
     
-//    std::vector<index_t> coboundary;
-//    get_simplex_coboundary( 1, 2, n, binomial_coeff, std::back_inserter(coboundary) );
-//    std::cout << coboundary <<std::endl;
-//    
-//    coboundary.clear();
-//    simplex_coboundary_enumerator coboundaries(1, 2, n, binomial_coeff);
-//    while (coboundaries.has_next())
-//        coboundary.push_back(coboundaries.next());
-//    std::cout << coboundary <<std::endl;
-   
-    
     std::vector<index_t> columns_to_reduce;
     
     
@@ -833,30 +610,7 @@ int main( int argc, char** argv ) {
 
 
     std::vector<value_t>& diameters = dist.distances;
-    
-//    std::cout << "precomputing 1-simplex diameters" << std::endl;
-//    
-//    std::vector<value_t> diameters( binomial_coeff( n, 2 ) );
-//    
-//    std::vector<index_t> edge_vertices(2);
-//    for (index_t edge = 0; edge < diameters.size(); ++edge) {
-//        #ifdef INDICATE_PROGRESS
-//        std::cout << "\033[Kstoring diameter of simplex " << edge + 1 << "/" << diameters.size() << std::flush << "\r";
-//        #endif
-//        
-//        edge_vertices.clear();
-//        get_simplex_vertices( edge, 1, n, binomial_coeff, std::back_inserter(edge_vertices) );
-//        diameters[edge] = dist(edge_vertices[0], edge_vertices[1]);
-//    }
-//    #ifdef INDICATE_PROGRESS
-//    std::cout << "\033[K";
-//    #endif
-
     #endif
-
-    
-//    std::cout << "diameters: " << diameters << std::endl;
-    
     
     for (index_t dim = 0; dim < dim_max; ++dim) {
     
@@ -879,7 +633,6 @@ int main( int argc, char** argv ) {
             binomial_coeff
         );
         
-        
         index_t num_simplices = binomial_coeff(n, dim + 2);
         
         columns_to_reduce.clear();
@@ -893,8 +646,6 @@ int main( int argc, char** argv ) {
         
         std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp);
 
-        
-        
         #ifdef PRECOMPUTE_DIAMETERS
         previous_diameters = std::move(diameters);
         
@@ -906,7 +657,6 @@ int main( int argc, char** argv ) {
         std::cout << "precomputing " << dim + 2 << "-simplex diameters" << std::endl;
         diameters = get_diameters( dist, dim + 2, previous_diameters, binomial_coeff );
         #endif
-        
     }
     
     index_t dim = dim_max;