blob: e562443f36ba2f46c55f49d312285c44121d675c (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#include <topology/rips.h>
#include <topology/filtration.h>
#include <geometry/l2distance.h>
#include <geometry/distances.h>
#include <utilities/containers.h> // for BackInsertFunctor
#include <vector>
#include <phat/boundary_matrix.h>
typedef PairwiseDistances<PointContainer, L2Distance> Pair_distances;
typedef Pair_distances::DistanceType Distance_type;
typedef Pair_distances::IndexType Vertex;
typedef Rips<Pair_distances> Generator;
typedef Generator::Simplex Smplx;
typedef Filtration<Smplx> Fltr;
int main(int argc, char* argv[])
{
Dimension skeleton;
Distance_type max_distance;
std::string infilename;
if(argc<2) {
std::cerr << "Requires inputfile" << std::endl;
std::exit(1);
}
infilename = argv[1];
if(argc>=3) {
skeleton=atoi(argv[2]);
if(skeleton==0) {
std::cerr << "# Command line argument 0 ignored" << std::endl;
skeleton=std::numeric_limits<Dimension>::max();
}
} else {
skeleton=std::numeric_limits<Dimension>::max();
}
if(argc>=4) {
max_distance=atof(argv[3]);
} else {
max_distance=std::numeric_limits<Distance_type>::max();
}
PointContainer points;
read_points(infilename, points);
Pair_distances distances(points);
Generator rips(distances);
Generator::Evaluator size(distances);
Fltr f;
// Generate 2-skeleton of the Rips complex for epsilon = 50
rips.generate(skeleton, max_distance, make_push_back_functor(f));
std::cerr << "# Generated complex of size: " << f.size() << std::endl;
// Generate filtration with respect to distance and compute its persistence
f.sort(Generator::Comparison(distances));
std::map<Smplx,phat::index,Smplx::VertexComparison> simplex_map;
phat::index size_of_simplex_map=0;
phat::boundary_matrix< phat::vector_vector > boundary_matrix;
boundary_matrix.set_num_cols( f.size() );
for(Fltr::Index it=f.begin();it!=f.end();it++) {
phat::column boundary_indices;
const Smplx& c = f.simplex(it);
for(Smplx::BoundaryIterator bit = c.boundary_begin(); bit != c.boundary_end(); bit++)
boundary_indices.push_back( simplex_map[*bit] );
std::sort(boundary_indices.begin(),boundary_indices.end());
boundary_matrix.set_col( size_of_simplex_map, boundary_indices );
phat::dimension dim_of_column = boundary_indices.size()==0 ? 0 : boundary_indices.size()-1;
boundary_matrix.set_dim( size_of_simplex_map, dim_of_column );
simplex_map[c] = size_of_simplex_map++;
}
boundary_matrix.save_binary( "rips.bin" );
return 1;
}
|