1 --- a/examples/rips/rips.cpp Fri May 01 15:09:38 2009 -0700
2 +++ b/examples/rips/rips.cpp Mon Jul 25 23:21:29 2011 -0700
3 @@ -25,10 +25,9 @@
4 //typedef Rips<ExplicitDistances<Distances> > Generator;
5 typedef Rips<Distances> Generator;
6 typedef Generator::Simplex Smplx;
7 -typedef std::vector<Smplx> SimplexVector;
8 -typedef Filtration<SimplexVector, unsigned> Fltr;
9 -//typedef StaticPersistence<> Persistence;
10 -typedef DynamicPersistenceChains<> Persistence;
11 +typedef Filtration<Smplx> Fltr;
12 +typedef StaticPersistence<> Persistence;
13 +// typedef DynamicPersistenceChains<> Persistence;
14 typedef PersistenceDiagram<> PDgm;
15
16
17 @@ -50,26 +49,25 @@
18
19 Generator rips(distances);
20 Generator::Evaluator size(distances);
21 - SimplexVector complex;
22 + Fltr f;
23
24 // Generate 2-skeleton of the Rips complex for epsilon = 50
25 - rips.generate(2, 10, make_push_back_functor(complex));
26 - std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
27 - rInfo("Generated complex of size: %d", complex.size());
28 + rips.generate(2, 10, make_push_back_functor(f));
29 + rInfo("Generated complex of size: %d", f.size());
30
31 // Generate filtration with respect to distance and compute its persistence
32 - Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
33 + f.sort(Generator::Comparison(distances));
34 Persistence p(f);
35 p.pair_simplices();
36 rInfo("Simplices paired");
37
38 + Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
39 +
40 // Record the persistence intervals in the persistence diagrams
41 std::map<Dimension, PDgm> dgms;
42 init_diagrams(dgms, p.begin(), p.end(),
43 - evaluate_through_map(make_offset_map(p.begin(), f.begin()),
44 - evaluate_through_filtration(f, size)),
45 - evaluate_through_map(make_offset_map(p.begin(), f.begin()),
46 - evaluate_through_filtration(f, Smplx::DimensionExtractor())));
47 + evaluate_through_map(m, size),
48 + evaluate_through_map(m, Smplx::DimensionExtractor()));
49
50 // Serialize the diagrams to a file
51 std::ofstream ofs("rips-diagrams");
52 @@ -77,33 +75,32 @@
53 oa << dgms;
54
55 // Output cycles
56 - for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
57 + for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
58 {
59 - Persistence::Cycle& cycle = cur->cycle;
60 + const Persistence::Cycle& cycle = cur->cycle;
61
62 if (!cur->sign()) // only negative simplices have non-empty cycles
63 {
64 Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
65
66 - const Smplx& b = f.simplex(f.begin() + (birth - p.begin())); // eventually this will be encapsulated behind an interface
67 - const Smplx& d = f.simplex(f.begin() + (cur - p.begin()));
68 + const Smplx& b = m[birth];
69 + const Smplx& d = m[cur];
70
71 if (b.dimension() != 1) continue;
72 std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
73 - } else if (cur->pair == cur) // positive could be unpaired
74 + } else if (cur->unpaired()) // positive could be unpaired
75 {
76 - const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
77 + const Smplx& b = m[cur];
78 if (b.dimension() != 1) continue;
79
80 std::cout << "Unpaired birth: " << size(b) << std::endl;
81 - cycle = cur->chain;
82 + // cycle = cur->chain;
83 }
84
85 // Iterate over the cycle
86 - for (Persistence::Cycle::const_iterator si = cycle.begin();
87 - si != cycle.end(); ++si)
88 + for (Persistence::Cycle::const_iterator si = cycle.begin(); si != cycle.end(); ++si)
89 {
90 - const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
91 + const Smplx& s = m[*si];
92 //std::cout << s.dimension() << std::endl;
93 const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
94 AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");