1 --- a/examples/rips/rips-pairwise.cpp Sat Nov 28 16:45:42 2009 -0800
2 +++ b/examples/rips/rips-pairwise.cpp Mon Jul 25 23:21:29 2011 -0700
3 @@ -21,10 +21,9 @@
4
5 typedef Rips<PairDistances> 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 void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name);
17 @@ -45,44 +44,45 @@
18 PairDistances distances(points);
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(skeleton, max_distance, make_push_back_functor(complex));
26 - std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
27 - std::cout << "Generated complex of size: " << complex.size() << std::endl;
28 + rips.generate(skeleton, max_distance, make_push_back_functor(f));
29 + std::cout << "# Generated complex of size: " << f.size() << std::endl;
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
35 Timer persistence_timer; persistence_timer.start();
36 Persistence p(f);
37 p.pair_simplices();
38 persistence_timer.stop();
39
40 +#if 1
41 // Output cycles
42 - for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
43 + Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
44 + for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
45 {
46 - Persistence::Cycle& cycle = cur->cycle;
47 + const Persistence::Cycle& cycle = cur->cycle;
48
49 if (!cur->sign()) // only negative simplices have non-empty cycles
50 {
51 Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
52
53 - const Smplx& b = f.simplex(f.begin() + (birth - p.begin())); // eventually this will be encapsulated behind an interface
54 - const Smplx& d = f.simplex(f.begin() + (cur - p.begin()));
55 + const Smplx& b = m[birth];
56 + const Smplx& d = m[cur];
57
58 // if (b.dimension() != 1) continue;
59 // std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
60 if (b.dimension() >= skeleton) continue;
61 diagram_out << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
62 - } else if (cur->pair == cur) // positive could be unpaired
63 + } else if (cur->unpaired()) // positive could be unpaired
64 {
65 - const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
66 + const Smplx& b = m[cur];
67 // if (b.dimension() != 1) continue;
68
69 // std::cout << "Unpaired birth: " << size(b) << std::endl;
70 - // cycle = cur->chain;
71 + // cycle = cur->chain; // TODO
72 if (b.dimension() >= skeleton) continue;
73 diagram_out << b.dimension() << " " << size(b) << " inf" << std::endl;
74 }
75 @@ -91,13 +91,14 @@
76 // for (Persistence::Cycle::const_iterator si = cycle.begin();
77 // si != cycle.end(); ++si)
78 // {
79 - // const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
80 + // const Smplx& s = m[*si];
81 // //std::cout << s.dimension() << std::endl;
82 // const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
83 // AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
84 // std::cout << vertices[0] << " " << vertices[1] << std::endl;
85 // }
86 }
87 +#endif
88
89 persistence_timer.check("# Persistence timer");
90 }