examples/cech-complex/cech-complex.cpp
 author Dmitriy Morozov Fri, 01 May 2009 15:36:25 -0700 branch dev changeset 124 af54138b98fd parent 100 884f70adc576 child 179 d15c6d144645 permissions -rw-r--r--
```
#include <utilities/log.h>

#include <topology/simplex.h>
#include <topology/filtration.h>
#include <topology/static-persistence.h>
//#include <topology/dynamic-persistence.h>
#include <topology/persistence-diagram.h>
#include "Miniball_dynamic_d.h"
#include <algorithm>

#include <iostream>
#include <fstream>

typedef         std::vector<Point>                                      PointContainer;
typedef         unsigned int                                            PointIndex;
typedef         Simplex<PointIndex, double>                             Smplx;
typedef         std::vector<Smplx>                                      SimplexVector;
typedef         Filtration<SimplexVector>                               CechFiltration;
typedef         StaticPersistence<>                                     Persistence;
//typedef         DynamicPersistenceTrails<>                              Persistence;
typedef         PersistenceDiagram<>            PDgm;

int choose(int n, int k)
{
if (k > n/2) k = n-k;
int numerator = 1, denominator = 1;
for (int i = 0; i < k; ++i)
{ numerator *= (n-i); denominator *= (i+1); }
return numerator/denominator;
}

void add_simplices(SimplexVector& sv, int d, const PointContainer& points)
{
PointIndex indices[d+1];
for (int i = 0; i < d+1; ++i)
indices[i] = d - i;

while(indices[d] < points.size() - d)
{
Miniball mb(points[indices[0]].dim());
Smplx s;
for (int i = 0; i < d+1; ++i)
{
mb.check_in(points[indices[i]]);
}
mb.build();
sv.push_back(s);

for (int i = 0; i < d+1; ++i)
{
++indices[i];
if (indices[i] < points.size() - i)
{
for (int j = i-1; j >= 0; --j)
indices[j] = indices[j+1] + 1;
break;
}
}
}
}

int main(int argc, char** argv)
{
#ifdef LOGGING

stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
#endif

SetFrequency(GetCounter("persistence/pair"), 10000);
SetTrigger(GetCounter("persistence/pair"), GetCounter(""));

// Read in the point set and compute its Delaunay triangulation
std::istream& in = std::cin;
int ambient_d, homology_d;
in >> ambient_d >> homology_d;

rInfo("Ambient dimension: %d", ambient_d);
rInfo("Will compute PD up to dimension: %d", homology_d);

PointContainer points;
while(in)
{
Point p(ambient_d);
for (int i = 0; i < ambient_d; ++i)
in >> p[i];
points.push_back(p);
}

// Compute simplices with their Cech values
SimplexVector    sv;
int num_simplices = 0;
for (int i = 0; i <= homology_d + 1; ++i)
num_simplices += choose(points.size(), i+1);
sv.reserve(num_simplices);
rInfo("Reserved SimplexVector of size: %d", num_simplices);

for (int i = 0; i <= homology_d + 1; ++i)
rInfo("Size of SimplexVector: %d", sv.size());

std::sort(sv.begin(), sv.end(), Smplx::VertexComparison());

// Set up the filtration
rInfo("Filtration initialized");

// Compute persistence
Persistence p(cf);
rInfo("Persistence initialized");
p.pair_simplices();
rInfo("Simplices paired");

std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
evaluate_through_map(OffsetMap<Persistence::OrderIndex, CechFiltration::Index>(p.begin(), cf.begin()),
evaluate_through_filtration(cf, Smplx::DataEvaluator())),
evaluate_through_map(OffsetMap<Persistence::OrderIndex, CechFiltration::Index>(p.begin(), cf.begin()),
evaluate_through_filtration(cf, Smplx::DimensionExtractor())));

for (int i = 0; i <= homology_d; ++i)
{
std::cout << i << std::endl << dgms[i] << std::endl;
}
}

```