--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-explicit-cohomology.cpp Tue Jan 26 11:14:50 2010 -0800
@@ -0,0 +1,201 @@
+#include <topology/cohomology-persistence.h>
+#include <topology/rips.h>
+
+#include <geometry/l2distance.h>
+#include <geometry/distances.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/property-maps.h>
+#include <utilities/timer.h>
+#include <utilities/log.h>
+#include <utilities/counter.h>
+#include <utilities/memory.h>
+
+#include <sys/resource.h>
+#include <malloc.h>
+
+#include <string>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "wrappers.h"
+
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef boost::tuple<Dimension, DistanceType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo, Wrapper<unsigned> > Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+typedef Persistence::CocyclePtr CocyclePtr;
+
+typedef ExplicitDistances< PairDistances > Distances;
+typedef Rips<Distances, Simplex<Vertex, Index> > Generator;
+typedef Generator::Simplex Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef SimplexVector::const_iterator SV_const_iterator;
+
+typedef std::map<Smplx, Index, Smplx::VertexComparison> Complex;
+
+#include "output.h" // for output_*()
+
+void program_options(int argc, char* argv[], std::string& infilename, unsigned& pt_size, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ Dimension skeleton;
+ DistanceType max_distance;
+ ZpField::Element prime;
+ unsigned pt_size;
+ std::string infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
+
+ program_options(argc, argv, infilename, pt_size, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
+ std::ofstream bdry_out(boundary_name.c_str());
+ std::ofstream vertices_out(vertices_name.c_str());
+ std::ofstream diagram_out(diagram_name.c_str());
+ std::cout << "# of points: " << pt_size << std::endl;
+ std::cout << "Boundary matrix: " << boundary_name << std::endl;
+ std::cout << "Cocycles: " << cocycle_prefix << "*.ccl" << std::endl;
+ std::cout << "Vertices: " << vertices_name << std::endl;
+ std::cout << "Diagram: " << diagram_name << std::endl;
+
+ Timer total_timer; total_timer.start();
+
+ // Read the pairwise distances
+ std::ifstream in(infilename.c_str());
+ Distances distances(pt_size);
+ for (unsigned i = 0; i < pt_size; ++i)
+ for (unsigned j = 0; j < pt_size; ++j)
+ in >> distances(i,j);
+
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ Generator::Comparison cmp(distances);
+ SimplexVector v;
+
+ Timer rips_timer; rips_timer.start();
+ rips.generate(skeleton, max_distance, make_push_back_functor(v));
+ std::sort(v.begin(), v.end(), Smplx::VertexComparison());
+
+ std::vector<unsigned> index_in_v(v.size());
+ for (unsigned idx = 0; idx < v.size(); ++idx)
+ index_in_v[idx] = idx;
+ std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
+
+ BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
+
+ rips_timer.stop();
+ std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
+
+ output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
+ output_vertex_indices(vertices_out, v);
+
+ Timer persistence_timer; persistence_timer.start();
+ ZpField zp(prime);
+ Persistence p(zp);
+ boost::progress_display show_progress(v.size());
+
+ for (unsigned j = 0; j < index_in_v.size(); ++j)
+ {
+ SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
+ std::vector<Index> boundary;
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
+ boundary.push_back(map_of_v[*bcur]->data());
+
+ Index idx; Death d; CocyclePtr ccl;
+ bool store = cur->dimension() < skeleton;
+ boost::tie(idx, d, ccl) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
+
+ // c[*cur] = idx;
+ if (store)
+ map_of_v[*cur]->data() = idx;
+
+ if (d && (size(*cur) - d->get<1>()) > 0)
+ {
+ AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+ diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+ }
+ ++show_progress;
+ }
+ // output infinte persistence pairs
+ for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+ diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
+ persistence_timer.stop();
+
+
+ // p.show_cocycles();
+ // Output alive cocycles of dimension 1
+ if (!cocycle_prefix.empty())
+ {
+ unsigned i = 0;
+ for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ if (cur->birth.get<0>() != 1) continue;
+ output_cocycle(cocycle_prefix, i, v, cur->birth, cur->zcolumn, prime);
+ // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+ ++i;
+ }
+ }
+ total_timer.stop();
+ rips_timer.check("Rips timer");
+ persistence_timer.check("Persistence timer");
+ total_timer.check("Total timer");
+}
+
+void program_options(int argc, char* argv[], std::string& infilename, unsigned& pt_size, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose Rips zigzag we want to compute")
+ ("num-pts", po::value<unsigned>(&pt_size), "Number of points in the input");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("prime,p", po::value<ZpField::Element>(&prime)->default_value(11), "Prime p for the field F_p")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction")
+ ("boundary,b", po::value<std::string>(&boundary_name), "Filename where to output the boundary matrix")
+ ("cocycle,c", po::value<std::string>(&cocycle_prefix), "Prefix of the filename where to output the 1-dimensional cocycles")
+ ("vertices,v", po::value<std::string>(&vertices_name), "Filename where to output the simplex-vertex mapping")
+ ("diagram,d", po::value<std::string>(&diagram_name), "Filename where to output the persistence diagram");
+#if LOGGING
+ std::vector<std::string> log_channels;
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
+
+ po::positional_options_description pos;
+ pos.add("input-file", 1);
+ pos.add("num-pts", 1);
+
+ po::options_description all; all.add(visible).add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(pos).run(), vm);
+ po::notify(vm);
+
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
+ if (vm.count("help") || !vm.count("input-file")|| !vm.count("num-pts"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file num-pts" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}