Modified cohomology code (unweighted and weighted) so that it doesn't use std::map for boundary computations, which should optimize memory usage. dev
authorChristos Mantoulidis <cmad@stanford.edu>
Tue Aug 04 16:58:44 2009 -0700 (2009-08-04)
branchdev
changeset 157700cbac5b23c
parent 156 f75fb57d2831
child 158 ce2aa05994f0
Modified cohomology code (unweighted and weighted) so that it doesn't use std::map for boundary computations, which should optimize memory usage.
examples/cohomology/CMakeLists.txt
examples/cohomology/output.h
examples/cohomology/rips-pairwise-cohomology.cpp
examples/cohomology/rips-weighted-cohomology.cpp
examples/cohomology/wrappers.h
     1.1 --- a/examples/cohomology/CMakeLists.txt	Tue Aug 04 13:23:16 2009 -0700
     1.2 +++ b/examples/cohomology/CMakeLists.txt	Tue Aug 04 16:58:44 2009 -0700
     1.3 @@ -1,7 +1,7 @@
     1.4  set                         (targets                        
     1.5                               rips-cohomology
     1.6                               rips-pairwise-cohomology
     1.7 -			     rips-weighted-cohomology
     1.8 +                             rips-weighted-cohomology
     1.9                               triangle-cohomology
    1.10                              )
    1.11                               
     2.1 --- a/examples/cohomology/output.h	Tue Aug 04 13:23:16 2009 -0700
     2.2 +++ b/examples/cohomology/output.h	Tue Aug 04 16:58:44 2009 -0700
     2.3 @@ -13,14 +13,16 @@
     2.4      return cmp(s1, s2) || cmp(s2, s1);
     2.5  }
     2.6  
     2.7 -unsigned index(const SimplexVector& v, const Smplx& s, const Generator::Comparison& cmp)
     2.8 +template<class Comparison>
     2.9 +unsigned index(const SimplexVector& v, const Smplx& s, const Comparison& cmp)
    2.10  {
    2.11      SimplexVector::const_iterator it = std::lower_bound(v.begin(), v.end(), s, cmp);
    2.12      while (neq(*it, s)) ++it;
    2.13      return it - v.begin();
    2.14  }
    2.15  
    2.16 -void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Generator::Comparison& cmp)
    2.17 +template<class Comparison>
    2.18 +void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Comparison& cmp)
    2.19  {
    2.20      unsigned i = 0;
    2.21      for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
    2.22 @@ -49,7 +51,7 @@
    2.23      }
    2.24  }
    2.25  
    2.26 -void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime, const Generator::Comparison& cmp)
    2.27 +void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime)
    2.28  {
    2.29      std::ostringstream istr; istr << '-' << i;
    2.30      std::string filename = cocycle_prefix + istr.str() + ".ccl";
    2.31 @@ -57,8 +59,8 @@
    2.32      out << "# Cocycle born at " << c.birth.get<1>() << std::endl;
    2.33      for (Persistence::ZColumn::const_iterator zcur = c.zcolumn.begin(); zcur != c.zcolumn.end(); ++zcur)
    2.34      {
    2.35 -        const Smplx& s = **(zcur->si);
    2.36 +        //const Smplx& s = **(zcur->si);
    2.37          out << (zcur->coefficient <= prime/2 ? zcur->coefficient : zcur->coefficient - prime) << " ";
    2.38 -        out << index(v, s, cmp) << "\n";
    2.39 +        out << zcur->si->getValue() << "\n";
    2.40      }
    2.41  }
     3.1 --- a/examples/cohomology/rips-pairwise-cohomology.cpp	Tue Aug 04 13:23:16 2009 -0700
     3.2 +++ b/examples/cohomology/rips-pairwise-cohomology.cpp	Tue Aug 04 16:58:44 2009 -0700
     3.3 @@ -5,6 +5,7 @@
     3.4  #include <geometry/distances.h>
     3.5  
     3.6  #include <utilities/containers.h>           // for BackInsertFunctor
     3.7 +#include <utilities/property-maps.h>
     3.8  #include <utilities/timer.h>
     3.9  #include <utilities/log.h>
    3.10  
    3.11 @@ -13,19 +14,22 @@
    3.12  #include <boost/tuple/tuple.hpp>
    3.13  #include <boost/program_options.hpp>
    3.14  
    3.15 +#include "wrappers.h"
    3.16 +
    3.17  typedef     PairwiseDistances<PointContainer, L2Distance>           PairDistances;
    3.18  typedef     PairDistances::DistanceType                             DistanceType;
    3.19  typedef     PairDistances::IndexType                                Vertex;
    3.20   
    3.21 -typedef     Rips<PairDistances>                                     Generator;
    3.22 +typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
    3.23 +typedef     CohomologyPersistence<BirthInfo, Wrapper<unsigned> >    Persistence;
    3.24 +typedef     Persistence::SimplexIndex                               Index;
    3.25 +typedef     Persistence::Death                                      Death;
    3.26 +
    3.27 +typedef     Rips<PairDistances, Simplex<Vertex, Index> >            Generator;
    3.28  typedef     Generator::Simplex                                      Smplx;
    3.29  typedef     std::vector<Smplx>                                      SimplexVector;
    3.30  typedef     SimplexVector::const_iterator                           SV_const_iterator;
    3.31  
    3.32 -typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
    3.33 -typedef     CohomologyPersistence<BirthInfo, SV_const_iterator>     Persistence;
    3.34 -typedef     Persistence::SimplexIndex                               Index;
    3.35 -typedef     Persistence::Death                                      Death;
    3.36  typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
    3.37  
    3.38  #include "output.h"         // for output_*()
    3.39 @@ -63,33 +67,41 @@
    3.40      Generator::Evaluator    size(distances);
    3.41      Generator::Comparison   cmp(distances);
    3.42      SimplexVector           v;
    3.43 -    Complex                 c;
    3.44      
    3.45      Timer rips_timer; rips_timer.start();
    3.46      rips.generate(skeleton, max_distance, make_push_back_functor(v));
    3.47 -    std::sort(v.begin(), v.end(), cmp);
    3.48 +    std::sort(v.begin(), v.end(), Smplx::VertexComparison());
    3.49 +
    3.50 +    std::vector<unsigned> index_in_v(v.size());
    3.51 +    for (unsigned idx = 0; idx < v.size(); ++idx)
    3.52 +        index_in_v[idx] = idx;
    3.53 +    std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
    3.54 +
    3.55 +    BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
    3.56 +
    3.57      rips_timer.stop();
    3.58      std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
    3.59  
    3.60 -    output_boundary_matrix(bdry_out, v, cmp);
    3.61 +    output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
    3.62      output_vertex_indices(vertices_out, v);
    3.63  
    3.64      Timer persistence_timer; persistence_timer.start();
    3.65      ZpField                 zp(prime);
    3.66      Persistence             p(zp);
    3.67 -    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
    3.68 +    for (unsigned j = 0; j < index_in_v.size(); ++j)
    3.69      {
    3.70 +        SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
    3.71          std::vector<Index>      boundary;
    3.72          for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
    3.73 -            boundary.push_back(c[*bcur]);
    3.74 +            boundary.push_back(map_of_v[*bcur]->data());
    3.75          
    3.76          Index idx; Death d;
    3.77          bool store = cur->dimension() < skeleton;
    3.78 -        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
    3.79 +        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
    3.80          
    3.81          // c[*cur] = idx;
    3.82          if (store)
    3.83 -            c.insert(std::make_pair(*cur, idx));
    3.84 +            map_of_v[*cur]->data() = idx;
    3.85  
    3.86          if (d && (size(*cur) - d->get<1>()) > 0)
    3.87          {
    3.88 @@ -109,7 +121,7 @@
    3.89      for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
    3.90      {
    3.91          if (cur->birth.get<0>() != 1) continue;
    3.92 -        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
    3.93 +        output_cocycle(cocycle_prefix, i, v, *cur, prime);
    3.94          // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
    3.95          ++i;
    3.96      }
     4.1 --- a/examples/cohomology/rips-weighted-cohomology.cpp	Tue Aug 04 13:23:16 2009 -0700
     4.2 +++ b/examples/cohomology/rips-weighted-cohomology.cpp	Tue Aug 04 16:58:44 2009 -0700
     4.3 @@ -5,6 +5,7 @@
     4.4  #include <geometry/distances.h>
     4.5  
     4.6  #include <utilities/containers.h>           // for BackInsertFunctor
     4.7 +#include <utilities/property-maps.h>
     4.8  #include <utilities/timer.h>
     4.9  #include <utilities/log.h>
    4.10  
    4.11 @@ -14,21 +15,22 @@
    4.12  #include <boost/program_options.hpp>
    4.13  #include <boost/progress.hpp>
    4.14  
    4.15 +#include "wrappers.h"
    4.16 +
    4.17  typedef     PairwiseDistances<PointContainer, WeightedL2Distance>   PairDistances;
    4.18  typedef     PairDistances::DistanceType                             DistanceType;
    4.19  typedef     PairDistances::IndexType                                Vertex;
    4.20   
    4.21 -typedef     WeightedRips<PairDistances>                             Generator;
    4.22 +typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
    4.23 +typedef     CohomologyPersistence<BirthInfo, Wrapper<unsigned> >    Persistence;
    4.24 +typedef     Persistence::SimplexIndex                               Index;
    4.25 +typedef     Persistence::Death                                      Death;
    4.26 +
    4.27 +typedef     WeightedRips<PairDistances, Simplex<Vertex, Index> >    Generator;
    4.28  typedef     Generator::Simplex                                      Smplx;
    4.29  typedef     std::vector<Smplx>                                      SimplexVector;
    4.30  typedef     SimplexVector::const_iterator                           SV_const_iterator;
    4.31  
    4.32 -typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
    4.33 -typedef     CohomologyPersistence<BirthInfo, SV_const_iterator>     Persistence;
    4.34 -typedef     Persistence::SimplexIndex                               Index;
    4.35 -typedef     Persistence::Death                                      Death;
    4.36 -typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
    4.37 -
    4.38  #include "output.h"         // for output_*()
    4.39  
    4.40  void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name);
    4.41 @@ -64,41 +66,53 @@
    4.42      Generator::Evaluator    size(distances);
    4.43      Generator::Comparison   cmp(distances);
    4.44      SimplexVector           v;
    4.45 -    Complex                 c;
    4.46 -    
    4.47 +
    4.48      Timer rips_timer; rips_timer.start();
    4.49      rips.generate(skeleton, max_distance, make_push_back_functor(v));
    4.50 -    std::sort(v.begin(), v.end(), cmp);
    4.51 +
    4.52 +    /* Keep simplices sorted lexicographically (so that we can binary search through them) */
    4.53 +    std::sort(v.begin(), v.end(), Smplx::VertexComparison());
    4.54 +
    4.55 +    /* We also need the simplices sorted by value though for the filtration:
    4.56 +       index_in_v[j] refers to the simplex v[index_in_v[j]]                      */
    4.57 +    std::vector<unsigned> index_in_v(v.size());
    4.58 +    for (unsigned idx = 0; idx < v.size(); ++idx)
    4.59 +        index_in_v[idx] = idx;
    4.60 +    std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
    4.61 +
    4.62 +    /* Set up map access to the lexicographically sorted simplices */
    4.63 +    BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
    4.64 +
    4.65      rips_timer.stop();
    4.66      std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
    4.67  
    4.68 -    /*output_boundary_matrix(bdry_out, v, cmp);
    4.69 +    /*output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
    4.70      output_vertex_indices(vertices_out, v);*/
    4.71  
    4.72      Timer persistence_timer; persistence_timer.start();
    4.73      ZpField                 zp(prime);
    4.74      Persistence             p(zp);
    4.75      boost::progress_display show_progress(v.size());
    4.76 -    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
    4.77 +    for (unsigned j = 0; j < index_in_v.size(); ++j)
    4.78      {
    4.79 +        SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
    4.80          std::vector<Index>      boundary;
    4.81          for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
    4.82 -            boundary.push_back(c[*bcur]);
    4.83 +            boundary.push_back(map_of_v[*bcur]->data());
    4.84          
    4.85          Index idx; Death d;
    4.86          bool store = cur->dimension() < skeleton;
    4.87 -        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
    4.88 +        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
    4.89          
    4.90 -        // c[*cur] = idx;
    4.91          if (store)
    4.92 -            c.insert(std::make_pair(*cur, idx));
    4.93 +            map_of_v[*cur]->data() = idx;
    4.94  
    4.95          if (d && (size(*cur) - d->get<1>()) > 0)
    4.96          {
    4.97              AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
    4.98              diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
    4.99          }
   4.100 -	++show_progress;
   4.101 +        ++show_progress;
   4.102      }
   4.103      // output infinte persistence pairs 
   4.104      for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
   4.105 @@ -112,7 +126,7 @@
   4.106      for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
   4.107      {
   4.108          if (cur->birth.get<0>() != 1) continue;
   4.109 -        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
   4.110 +        output_cocycle(cocycle_prefix, i, v, *cur, prime);
   4.111          // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
   4.112          ++i;
   4.113      }
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/examples/cohomology/wrappers.h	Tue Aug 04 16:58:44 2009 -0700
     5.3 @@ -0,0 +1,59 @@
     5.4 +/**
     5.5 + * Wrapper class
     5.6 + * 
     5.7 + * At points we need to wrap primitive data types in classes, so that we can
     5.8 + * pass them as template arguments to classes that like to inherit from their
     5.9 + * template arguments--this is done in CohomologyPersistence, for example.
    5.10 + */
    5.11 +
    5.12 +template<typename Primitive>
    5.13 +class    Wrapper
    5.14 +{
    5.15 +
    5.16 +    public:
    5.17 +
    5.18 +        Wrapper () {}
    5.19 +        Wrapper (Primitive v)                       { value = v;    }
    5.20 +
    5.21 +        void       setValue  (const Primitive &v)   { value = v;    }
    5.22 +        Primitive &getValue  ()                     { return value; }
    5.23 +
    5.24 +        /* provide seemless integration */
    5.25 +        Wrapper   &operator =(const Primitive &v)   { setValue(v);     return *this; }
    5.26 +        operator  Primitive()                       { return getValue;               }
    5.27 +
    5.28 +    protected:
    5.29 +
    5.30 +        Primitive value;
    5.31 +
    5.32 +};
    5.33 +
    5.34 +/**
    5.35 + * IndirectIndexComparison class
    5.36 + * 
    5.37 + * This class serves as a comparison function for arrays that are being sorted
    5.38 + * even though they only contain *indices* to the real data. Therefore, a reference
    5.39 + * to the original data as well as the data comparison function needs to be passed
    5.40 + * to this class for it to be functional.
    5.41 + */
    5.42 +
    5.43 +template<class DataContainer, class DataComparison>
    5.44 +class IndirectIndexComparison
    5.45 +{
    5.46 +
    5.47 +    public:
    5.48 +
    5.49 +        IndirectIndexComparison(const DataContainer &dstor, const DataComparison &dcmp) :
    5.50 +            container(dstor), comparison(dcmp) { }
    5.51 +
    5.52 +        bool operator()(const unsigned &idx_1, const unsigned &idx_2) const
    5.53 +        {
    5.54 +            return comparison(container[idx_1], container[idx_2]);
    5.55 +        }
    5.56 +
    5.57 +    private:
    5.58 +
    5.59 +        const DataContainer  &container;
    5.60 +        const DataComparison &comparison;
    5.61 +
    5.62 +};