Switched to a new architecture (Vineyard is a visitor for Filtration),
added LowerStarFilation, examples/grid (pdbdistance-vineyard in particular)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/README Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,19 @@
+Dependencies
+ CGAL-3.2 - for alpha-shapes and kinetic data structures
+ DSR-PDB - for reading in PDB files
+ scons - for controlling the build process
+ boost - great set of C++ libraries
+ Doxygen - for building documentation
+
+Configuration
+ Local configuration is read from config/HOSTNAME.py
+ The path to and architecture of CGAL are important variables (cgal_path and
+ cgal_architecture respectively).
+
+Building
+ To build examples run "scons" at the top level, or "scons -u" in the specific
+ examples/ subdirectory.
+ "scons debug=1" turns on debugging. "scons counters=1" turns on counters.
+
+Author
+ Dmitriy Morozov <morozov@cs.duke.edu>
--- a/SConstruct Sat Dec 16 15:34:46 2006 -0500
+++ b/SConstruct Thu Dec 21 13:32:34 2006 -0500
@@ -19,21 +19,23 @@
toolpath=["sconstools"])
Help (optns.GenerateHelpText(base_env))
-opt = base_env.Copy(CPPFLAGS = ['-O'])
dbg = base_env.Copy(CPPFLAGS = ['-g'])
-opt.Append (CPPDEFINES = ['NDEBUG'])
+opt = base_env.Copy(CPPDEFINES = ['NDEBUG'])
if ARGUMENTS.get('debug', 0):
- env = dbg
- env.Append(CPPDEFINES = ['CWDEBUG'])
- env.Append(LIBS = ['cwd'])
+ env = dbg
+ env.Append (CPPDEFINES = ['CWDEBUG'])
+ env.Append (LIBS = ['cwd'])
+elif ARGUMENTS.get('optimize', 0):
+ opt.Append (CPPFLAGS = ['-O' + str(ARGUMENTS['optimize'])])
+ env = opt
else:
- env = opt
+ opt.Append (CPPFLAGS = ['-O'])
+ env = opt
if ARGUMENTS.get('counters', 0):
env.Append (CPPDEFINES = ['COUNTERS'])
-
# Don't create .sconsign files all over the place
SConsignFile()
--- a/examples/SConscript Sat Dec 16 15:34:46 2006 -0500
+++ b/examples/SConscript Thu Dec 21 13:32:34 2006 -0500
@@ -1,6 +1,6 @@
Import('*')
-SConscript (['triangle/SConscript'],
+SConscript (['triangle/SConscript',
# 'alphashapes/SConscript',
-# 'grid/SConscript'],
+ 'grid/SConscript'],
exports=['env', 'external_sources'])
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/SConscript Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,13 @@
+Import('*')
+
+# Libraries
+libraries = ['dsrpdb']
+
+# Sources
+sources = ['test-grid2D.cpp', 'pdbdistance-vineyard.cpp']
+
+local_env = env.Copy()
+local_env.Append (LIBS = libraries)
+
+for s in sources:
+ Default(local_env.Program([s] + external_sources))
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/grid2D.h Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,83 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2006
+ */
+
+#ifndef __GRID2D_H__
+#define __GRID2D_H__
+
+#include <memory>
+#include <vector>
+#include <map>
+#include <set>
+#include <limits>
+//#include <cmath>
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/split_member.hpp>
+#include <boost/serialization/nvp.hpp>
+
+#include "types.h"
+
+#include <boost/serialization/export.hpp>
+
+/**
+ * Grid2D stores a grid
+ */
+class Grid2D
+{
+ public:
+ typedef RealType ValueType;
+ typedef unsigned int CoordinateIndex;
+
+ typedef std::vector<ValueType> ValueVector;
+
+ public:
+ Grid2D(CoordinateIndex xx = 1, CoordinateIndex yy = 1);
+
+ /// Sets the grid dimensions to (xx,yy)
+ void change_dimensions(CoordinateIndex xx, CoordinateIndex yy);
+
+ ValueType& operator()(CoordinateIndex i, CoordinateIndex j) { return data[i*x + j]; }
+ const ValueType& operator()(CoordinateIndex i, CoordinateIndex j) const { return data[i*x + j]; }
+ ValueType& operator()(CoordinateIndex i) { return data[i]; }
+ const ValueType& operator()(CoordinateIndex i) const { return data[i]; }
+
+ CoordinateIndex xsize() const { return x; }
+ CoordinateIndex ysize() const { return y; }
+ CoordinateIndex size() const { return x*y; }
+
+ /* Given a sequential index of an element return its coordinates */
+ CoordinateIndex xpos(CoordinateIndex i) const { return i / x; }
+ CoordinateIndex ypos(CoordinateIndex i) const { return i % x; }
+ CoordinateIndex seq(CoordinateIndex i, CoordinateIndex j) const;
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ static const CoordinateIndex INVALID_INDEX = -1;
+
+ private:
+ CoordinateIndex x,y;
+ ValueVector data;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ template<class Archive> void save(Archive& ar, version_type ) const;
+ template<class Archive> void load(Archive& ar, version_type );
+
+ BOOST_SERIALIZATION_SPLIT_MEMBER()
+#endif
+};
+//BOOST_CLASS_EXPORT(Grid2D)
+
+
+std::ostream& operator<<(std::ostream& out, const Grid2D& grid) { return grid.operator<<(out); }
+
+
+#include "grid2D.hpp"
+
+#endif // __GRID2D_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/grid2D.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,65 @@
+#include <iostream>
+#include <limits>
+
+/* Implementations */
+
+Grid2D::
+Grid2D(CoordinateIndex xx, CoordinateIndex yy):
+ x(xx), y(yy), data(x*y)
+{}
+
+void
+Grid2D::
+change_dimensions(CoordinateIndex xx, CoordinateIndex yy)
+{
+ x = xx; y = yy;
+ data.resize(x*y);
+}
+
+Grid2D::CoordinateIndex
+Grid2D::
+seq(CoordinateIndex i, CoordinateIndex j) const
+{
+ // Do not forget to check if less than 0, if Index is made signed --- dangerous
+ if (i >= x || j >= y)
+ return INVALID_INDEX;
+
+ return i*x + j;
+}
+
+std::ostream&
+Grid2D::
+operator<<(std::ostream& out) const
+{
+ for (Grid2D::CoordinateIndex i = 0; i < xsize(); ++i)
+ {
+ for (Grid2D::CoordinateIndex j = 0; j < ysize(); ++j)
+ std::cout << operator()(i, j) << ' ';
+ std::cout << std::endl;
+ }
+ return out;
+}
+
+#if 0
+using boost::serialization::make_nvp;
+
+template<class Archive>
+void
+Grid2D::
+save(Archive& ar, version_type ) const
+{
+ ar << BOOST_SERIALIZATION_NVP(x);
+ ar << BOOST_SERIALIZATION_NVP(y);
+ ar << make_nvp("data", data);
+}
+
+template<class Archive>
+void
+Grid2D::
+load(Archive& ar, version_type )
+{
+ ar >> make_nvp("x", x);
+ ar >> make_nvp("y", y);
+ ar >> make_nvp("data", data);
+}
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/grid2Dvineyard.h Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,198 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2006
+ */
+
+#ifndef __GRID2DVINEYARD_H__
+#define __GRID2DVINEYARD_H__
+
+#include "sys.h"
+#include "debug.h"
+
+#include "grid2D.h"
+#include "lowerstarfiltration.h"
+
+#include <CGAL/Kinetic/Inexact_simulation_traits_1.h>
+#include <CGAL/Kinetic/Sort.h>
+#include <CGAL/Kinetic/Sort_visitor_base.h>
+
+#include <vector>
+
+
+class Grid2DVineyard
+{
+ public:
+ typedef Grid2DVineyard Self;
+
+ class VertexType;
+ typedef std::vector<VertexType> VertexVector;
+ typedef VertexVector::iterator VertexIndex;
+
+ typedef LowerStarFiltration<VertexIndex> LSFiltration;
+
+ class StaticEvaluator;
+ class KineticEvaluator;
+ class VertexComparison;
+
+ typedef Grid2D::CoordinateIndex CoordinateIndex;
+ typedef Grid2D::ValueType ValueType;
+
+ typedef LSFiltration::Index Index;
+ typedef LSFiltration::Simplex Simplex;
+ typedef LSFiltration::VertexOrderIndex VertexOrderIndex;
+ typedef LSFiltration::Vineyard Vineyard;
+ typedef Vineyard::Evaluator Evaluator;
+
+ class SortVisitor;
+ typedef CGAL::Kinetic::Inexact_simulation_traits_1 Traits;
+ typedef CGAL::Kinetic::Sort<Traits, SortVisitor> Sort;
+ typedef Traits::Simulator Simulator;
+ typedef Traits::Active_points_1_table ActivePointsTable;
+ typedef ActivePointsTable::Key Key;
+ typedef std::map<Key, VertexOrderIndex> KeyOrderMap;
+
+ typedef std::vector<Grid2D*> GridStackContainer;
+
+ public:
+ Grid2DVineyard(Grid2D* g);
+ ~Grid2DVineyard();
+
+ void compute_pairing();
+ void compute_vineyard(Grid2D* grid);
+
+ Grid2D* grid() const { return grid_stack_.back(); }
+ Grid2D* grid(int i) const { return grid_stack_[i]; }
+ int num_grids() const { return grid_stack_.size(); }
+ const LSFiltration* filtration() const { return filtration_; }
+ const Vineyard* vineyard() const { return vineyard_; }
+
+ public:
+ // For Kinetic Sort
+ void swap(Key a, Key b);
+
+ protected:
+ // Do something cleverer
+ virtual bool neighbors(VertexIndex v1, VertexIndex v2) const { return true; }
+
+ private:
+ void add_simplices();
+ void change_evaluator(Evaluator* eval);
+
+ private:
+ GridStackContainer grid_stack_;
+ VertexVector vertices_;
+ LSFiltration* filtration_;
+ Vineyard* vineyard_;
+ Evaluator* evaluator_;
+
+ KeyOrderMap kinetic_map_;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ Grid2DVineyard() {}
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type )
+ {
+ ar & BOOST_SERIALIZATION_NVP(grid_stack_);
+ ar & BOOST_SERIALIZATION_NVP(vertices_);
+ ar & BOOST_SERIALIZATION_NVP(filtration_);
+ };
+#endif
+};
+
+//BOOST_CLASS_EXPORT(Grid2DVineyard)
+
+class Grid2DVineyard::VertexType
+{
+ public:
+ VertexType(CoordinateIndex ii = 0): i_(ii) {}
+
+ CoordinateIndex index() const { return i_; }
+ void set_index(CoordinateIndex i) { i_ = i; }
+ VertexOrderIndex get_order() const { return order_; }
+ void set_order(const VertexOrderIndex& o) { order_ = o; }
+
+ Key kinetic_key() const { return key_; }
+ void set_kinetic_key(Key k) { key_ = k; }
+
+ private:
+ CoordinateIndex i_;
+ VertexOrderIndex order_;
+ Key key_;
+};
+
+std::ostream& operator<<(std::ostream& out, const Grid2DVineyard::VertexIndex& vi) { return out << vi->index(); }
+
+class Grid2DVineyard::VertexComparison
+{
+ public:
+ VertexComparison(const Grid2D* g): grid(g) {}
+ bool operator()(VertexIndex i, VertexIndex j) const { return (*grid)(i->index()) <
+ (*grid)(j->index()); }
+
+ private:
+ const Grid2D* grid;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ VertexComparison() {}
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type ) { ar & BOOST_SERIALIZATION_NVP(grid); }
+#endif
+};
+
+class Grid2DVineyard::StaticEvaluator: public Evaluator
+{
+ public:
+ StaticEvaluator(Grid2D* grid, RealType time):
+ time_(time), grid_(grid) {}
+
+ virtual RealType time() { return time_; }
+ virtual RealType value(const Simplex& s) { return (*grid_)(s.get_attachment()->index()); }
+
+ private:
+ RealType time_;
+ Grid2D* grid_;
+};
+
+class Grid2DVineyard::KineticEvaluator: public Evaluator
+{
+ public:
+ KineticEvaluator(Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset):
+ sp_(sp), apt_(apt), time_offset_(time_offset) {}
+
+ virtual RealType time() { return time_offset_ + CGAL::to_double(get_time()); }
+ virtual RealType value(const Simplex& s) { return CGAL::to_double(apt_->at(s.get_attachment()->kinetic_key()).x()(get_time())); }
+
+ private:
+ Simulator::Time get_time() { return sp_->current_time(); }
+
+ Simulator::Handle sp_;
+ ActivePointsTable::Handle apt_;
+ RealType time_offset_;
+};
+
+
+class Grid2DVineyard::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
+{
+ public:
+ SortVisitor(Grid2DVineyard* gv): gv_(gv) {}
+
+ template<class Vertex_handle>
+ void before_swap(Vertex_handle a, Vertex_handle b) const { gv_->swap(*a,*b); }
+
+ private:
+ Grid2DVineyard* gv_;
+};
+
+#include "grid2Dvineyard.hpp"
+
+#endif // __GRID2DVINEYARD_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/grid2Dvineyard.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,141 @@
+/* Implementation */
+
+Grid2DVineyard::
+Grid2DVineyard(Grid2D* g): vertices_(g->size())
+{
+ grid_stack_.push_back(g);
+ for (CoordinateIndex i = 0; i < g->size(); ++i)
+ vertices_[i].set_index(i);
+
+ evaluator_ = new StaticEvaluator(grid(), 0);
+ vineyard_ = new Vineyard(evaluator_);
+
+ filtration_ = new LSFiltration(vertices_.begin(), vertices_.end(), VertexComparison(grid()), vineyard_);
+ add_simplices();
+}
+
+Grid2DVineyard::
+~Grid2DVineyard()
+{
+ delete filtration_;
+ delete vineyard_;
+ delete evaluator_;
+}
+
+void
+Grid2DVineyard::
+compute_pairing()
+{
+ filtration_->fill_simplex_index_map();
+ filtration_->pair_simplices(filtration_->begin(), filtration_->end());
+ vineyard_->start_vines(filtration_->begin(), filtration_->end());
+}
+
+void
+Grid2DVineyard::
+compute_vineyard(Grid2D* g)
+{
+ AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
+
+ typedef Traits::Kinetic_kernel::Point_1 Point;
+ typedef Traits::NT NT;
+ typedef Traits::Kinetic_kernel::Function_kernel::Construct_function CF;
+ typedef Traits::Kinetic_kernel::Motion_function F;
+
+ Traits tr;
+ Simulator::Handle sp = tr.simulator_handle();
+ ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
+ Sort sort(tr, SortVisitor(this));
+
+ // Setup the (linear) trajectories
+ CF cf;
+ kinetic_map_.clear();
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ {
+ ValueType val0 = (*grid())(cur->index());
+ ValueType val1 = (*g)(cur->index());
+ F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t
+ Point p(x);
+ cur->set_kinetic_key(apt->insert(p));
+ kinetic_map_[cur->kinetic_key()] = cur->get_order();
+ }
+
+ // Process all the events (compute the vineyard in the process)
+ change_evaluator(new KineticEvaluator(sp, apt, num_grids() - 1));
+ /*
+ while (sp->next_event_time() < 1)
+ {
+ std::cout << "Next event time: " << sp->next_event_time() << std::endl;
+ sp->set_current_event_number(sp->current_event_number() + 1);
+ std::cout << "Processed event" << std::endl;
+ }
+ */
+ sp->set_current_time(1.0);
+ std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
+
+ // Add the grid to the stack
+ grid_stack_.push_back(g);
+ change_evaluator(new StaticEvaluator(grid(), num_grids() - 1));
+}
+
+void
+Grid2DVineyard::
+swap(Key a, Key b)
+{
+ VertexOrderIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
+ AssertMsg(filtration_->get_vertex_cmp()(ao, bo), "In swap(a,b), a must precede b");
+ filtration_->transpose_vertices(ao);
+ AssertMsg(filtration_->get_vertex_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
+}
+
+void
+Grid2DVineyard::
+add_simplices()
+{
+ // Takes advantage of LowerStarFiltration's smart append (which allows faces
+ // to be inserted after cofaces, since everything is rearranged in the
+ // proper lower star order anyway). Also note that vertices were added by
+ // LowerStarFiltration's constructor
+ for (CoordinateIndex x = 0; x < grid()->xsize() - 1; ++x)
+ for (CoordinateIndex y = 0; y < grid()->ysize() - 1; ++y)
+ {
+ VertexIndex v(&vertices_[grid()->seq(x,y)]);
+ VertexIndex vh(&vertices_[grid()->seq(x+1,y)]);
+ VertexIndex vv(&vertices_[grid()->seq(x,y+1)]);
+ VertexIndex vd(&vertices_[grid()->seq(x+1,y+1)]);
+
+ Simplex sh(v);
+ sh.add(vh); filtration_->append(sh); // Horizontal edge
+ sh.add(vd); filtration_->append(sh); // "Horizontal" triangle
+
+ Simplex sv(v);
+ sv.add(vv); filtration_->append(sv); // Vertical edge
+ sv.add(vd); filtration_->append(sv); // "Vertical" triangle
+
+ Simplex sd(v);
+ sd.add(vd); filtration_->append(sd); // Diagonal edge
+
+ if (y == grid()->ysize() - 2)
+ {
+ Simplex s(vv);
+ s.add(vd); filtration_->append(s); // Top edge
+ }
+ if (x == grid()->xsize() - 2)
+ {
+ Simplex s(vh);
+ s.add(vd); filtration_->append(s); // Right edge
+ }
+ }
+}
+
+void
+Grid2DVineyard::
+change_evaluator(Evaluator* eval)
+{
+ AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
+
+ delete evaluator_;
+ evaluator_ = eval;
+ vineyard_->set_evaluator(evaluator_);
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/pdbdistance-vineyard.cpp Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,78 @@
+//#include <boost/archive/binary_oarchive.hpp>
+#include "sys.h"
+#include "debug.h"
+
+#include "pdbdistance.h"
+#include "grid2Dvineyard.h"
+
+#include <fstream>
+#include <string>
+#include <sstream>
+
+std::string frame_filename(const std::string& prefix, int frame, int subframe)
+{
+ std::ostringstream os;
+ os << prefix << frame << "_" << subframe << ".pdb";
+ return os.str();
+}
+
+int main(int argc, char** argv)
+{
+#ifdef CWDEBUG
+ Debug(dc::filtration.off());
+ Debug(dc::cycle.off());
+ Debug(dc::vineyard.off());
+ Debug(dc::transpositions.off());
+ Debug(dc::lsfiltration.off());
+
+ dionysus::debug::init();
+#endif
+
+ if (argc < 5)
+ {
+ std::cout << "Usage: pdbdistance FILENAME LASTFRAME LASTSUBFRAME OUTFILENAME" << std::endl;
+ std::cout << " FILENAME - prefix of the filenames of the PDB frames" << std::endl;
+ std::cout << " LASTFRAME - the last frame number" << std::endl;
+ std::cout << " LASTSUBFRAME - the last subframe number" << std::endl;
+ std::cout << " OUTFILENAME - filename prefix for the resulting vineyards" << std::endl;
+ std::cout << std::endl;
+ std::cout << "Computes a vineyard of the pairwise distance function for a sequence of PDB frames." << std::endl;
+ std::cout << "Frames are in files FILENAME#1_#2.pdb, where #1 goes from 0 to LASTFRAME, " << std::endl;
+ std::cout << "and #2 goes from 0 to LASTSUBFRAME." << std::endl;
+ exit(0);
+ }
+ std::string infilename = argv[1];
+ int lastframe; std::istringstream(argv[2]) >> lastframe;
+ int lastsubframe; std::istringstream(argv[3]) >> lastsubframe;
+ std::string outfilename = argv[4];
+
+ // Compute initial filtration
+ int f = 0; int sf = 0;
+ std::ifstream in(frame_filename(infilename, f, sf++).c_str());
+ Grid2DVineyard v(new PDBDistanceGrid(in));
+ in.close();
+ std::cout << "Filtration generated, size: " << v.filtration()->size() << std::endl;
+ v.compute_pairing();
+ std::cout << "Pairing computed" << std::endl;
+
+ // Process frames computing the vineyard
+ while (f <= lastframe)
+ {
+ std::string fn = frame_filename(infilename, f, sf++);
+ std::cout << "Processing " << fn << std::endl;
+ in.open(fn.c_str());
+ v.compute_vineyard(new PDBDistanceGrid(in));
+ in.close();
+ if (sf == lastsubframe) { sf = 0; ++f; }
+ }
+ std::cout << "Vineyard computed" << std::endl;
+
+ v.vineyard()->save_edges(outfilename);
+
+#if 0
+ std::ofstream ofs(outfilename.c_str(), std::ios::binary);
+ boost::archive::binary_oarchive oa(ofs);
+ oa << make_nvp("Filtration", pgf);
+ ofs.close();
+#endif
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/pdbdistance.h Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,75 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2006
+ *
+ * Depends on Daniel Russel's "simple C++ PDB" (aka DSR-PDB).
+ */
+
+#ifndef __PDBDISTANCE_H__
+#define __PDBDISTANCE_H__
+
+#include <fstream>
+#include <string>
+#include <dsrpdb/Protein.h>
+#include <dsrpdb/iterator.h>
+#include <cmath>
+
+#include <boost/serialization/access.hpp>
+
+#include "types.h"
+#include "grid2D.h"
+
+#include <boost/serialization/export.hpp>
+
+
+class PDBDistanceGrid: public Grid2D
+{
+ public:
+ PDBDistanceGrid()
+ {}
+
+ PDBDistanceGrid(std::istream& in)
+ {
+ load_stream(in);
+ }
+
+ void load_stream(std::istream& in)
+ {
+ dsrpdb::Protein p(in);
+ std::vector<dsrpdb::Point> CAs(ca_coordinates_begin(p), ca_coordinates_end(p));
+ std::cout << "CAs created, size: " << CAs.size() << std::endl;
+
+ Grid2D::change_dimensions(CAs.size(), CAs.size());
+ for (Grid2D::CoordinateIndex i = 0; i < CAs.size(); ++i)
+ for (Grid2D::CoordinateIndex j = 0; j < CAs.size(); ++j)
+ {
+ if (i < j)
+ Grid2D::operator()(i,j) = distance(CAs[i], CAs[j]);
+ else
+ Grid2D::operator()(i,j) = 0;
+ }
+ }
+
+ private:
+ Grid2D::ValueType distance(dsrpdb::Point p1, dsrpdb::Point p2) const
+ {
+ dsrpdb::Vector v = p1 - p2;
+ return std::sqrt(v*v);
+ }
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type version)
+ {
+ ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Grid2D);
+ }
+#endif
+};
+
+//BOOST_CLASS_EXPORT(PDBDistanceGrid)
+
+#endif // __PDBDISTANCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/test-grid2D.cpp Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,16 @@
+#include "grid2D.h"
+#include <iostream>
+
+int main()
+{
+ Grid2D grid(40,60);
+ int i = 0;
+ for (int x = 0; x < 40; ++x)
+ for (int y = 0; y < 60; ++y)
+ {
+ grid(x,y) = i++;
+ }
+
+ std::cout << grid(20,30) << std::endl;
+ std::cout << grid << std::endl;
+}
--- a/examples/triangle/SConscript Sat Dec 16 15:34:46 2006 -0500
+++ b/examples/triangle/SConscript Thu Dec 21 13:32:34 2006 -0500
@@ -1,20 +1,3 @@
Import('*')
Default(env.Program(['triangle.cpp'] + external_sources))
-
-## Common variables
-##libraries = ['dsrpdb', 'boost_serialization']
-#
-## Sources
-#sources = [['triangle.cpp', []]]
-##, 'grid2D.cpp', 'alphashapes.cpp', 'euclidean.cpp',
-## 'combustion_vineyard.cpp', 'folding.cpp', 'pdbdistance.cpp',
-## 'geodesic_vineyard.cpp']
-#
-#local_env = env.Copy()
-#
-## Build instructions
-#for sl in sources:
-# s,l = sl
-# o = local_env.Object(s)
-# t = Default(local_env.Program(o + external_sources))
--- a/examples/triangle/triangle.cpp Sat Dec 16 15:34:46 2006 -0500
+++ b/examples/triangle/triangle.cpp Thu Dec 21 13:32:34 2006 -0500
@@ -1,11 +1,10 @@
-#include "vineyard.h"
+#include "filtration.h"
#include "simplex.h"
#include <vector>
#include <iostream>
-typedef SimplexWithValue<int> Simplex;
-typedef Vineyard<Simplex> TriangleFiltration;
-//typedef Filtration<Simplex> TriangleFiltration;
+typedef SimplexWithValue<int> Simplex;
+typedef Filtration<Simplex> TriangleFiltration;
void fillTriangleSimplices(TriangleFiltration& f)
{
@@ -23,8 +22,6 @@
f.append(Simplex(bg + 1, bg + 3, 2.9)); // BC
f.append(Simplex(bg + 2, end, 3.5)); // CA
f.append(Simplex(bg, bg + 3, 5)); // ABC
-
- f.fill_simplex_index_map();
}
int main()
@@ -33,14 +30,19 @@
dionysus::debug::init();
Debug(dc::filtration.on());
- Debug(dc::cycle.on());
+ Debug(dc::cycle.off());
Debug(dc::vineyard.on());
+ Debug(dc::transpositions.on());
#endif
- TriangleFiltration tf;
+ Evaluator<Simplex> e;
+ TriangleFiltration::Vineyard v(&e);
+ TriangleFiltration tf(&v);
fillTriangleSimplices(tf);
- tf.pair_simplices(tf.begin());
+ tf.fill_simplex_index_map();
+ tf.pair_simplices(tf.begin(), tf.end());
+ v.start_vines(tf.begin(), tf.end());
std::cout << "Filtration size: " << tf.size() << std::endl;
std::cout << tf << std::endl;
@@ -52,6 +54,8 @@
std::cout << *tf.get_index(BC) << std::endl;
tf.transpose(tf.get_index(BC));
std::cout << tf;
+ std::cout << AB << std::endl;
+ std::cout << *tf.get_index(AB) << std::endl;
tf.transpose(tf.get_index(AB));
std::cout << tf;
#endif
--- a/include/circular_list.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/circular_list.h Thu Dec 21 13:32:34 2006 -0500
@@ -202,8 +202,8 @@
inline T& front() { return *begin(); }
inline const T& front() const { return *begin(); }
- inline T& back() { return m_node->data; }
- inline const T& back() const { return m_node->data; }
+ inline T& back() { return m_node->m_data; }
+ inline const T& back() const { return m_node->m_data; }
inline iterator begin() { return iterator(m_node->m_next); }
inline iterator end() { return iterator(m_node); }
--- a/include/consistencylist.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/consistencylist.h Thu Dec 21 13:32:34 2006 -0500
@@ -51,14 +51,16 @@
typedef T value_type;
typedef T& reference;
typedef const T& const_reference;
- typedef ConsistencyListIterator<T> iterator;
- typedef const_ConsistencyListIterator<T> const_iterator;
+ typedef ConsistencyListIterator<T> iterator;
+ typedef const_ConsistencyListIterator<T> const_iterator;
ConsistencyList(): sz(0) {}
~ConsistencyList() { clear(); }
/// \name Order operations
void swap(iterator i, iterator j); ///< Exchanges the order of simplices pointed to by i and j
+ template<class BinaryPredicate>
+ void sort(BinaryPredicate cmp); ///< Sorts the elements in accordance with cmp
/// \name Container operations
/// @{
--- a/include/consistencylist.hpp Sat Dec 16 15:34:46 2006 -0500
+++ b/include/consistencylist.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -9,6 +9,19 @@
}
template<class T>
+template<class BinaryPredicate>
+void
+ConsistencyList<T>::
+sort(BinaryPredicate cmp)
+{
+ Parent::sort(cmp);
+ OrderType cur_order = 0;
+ for (typename Parent::iterator cur = begin(); cur != end(); ++cur)
+ cur->consistency = cur_order++;
+}
+
+
+template<class T>
typename ConsistencyList<T>::iterator
ConsistencyList<T>::
push_back(const_reference x)
--- a/include/cycle.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/cycle.h Thu Dec 21 13:32:34 2006 -0500
@@ -63,7 +63,7 @@
/// \name Modifiers
/// @{
using CycleRepresentation::erase;
- void append(const_reference x);
+ void append(const_reference x, const ConsistencyComparison& cmp);
/// @}
/// \name Accessors
@@ -96,6 +96,7 @@
typedef std::list<Item> TemporaryCycleRepresenation;
using CycleRepresentation::push_back;
+ using CycleRepresentation::insert;
private:
size_t sz;
--- a/include/cycle.hpp Sat Dec 16 15:34:46 2006 -0500
+++ b/include/cycle.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -22,8 +22,23 @@
template<class I, class OrderCmp, class ConsistencyCmp>
void
Cycle<I,OrderCmp,ConsistencyCmp>::
-append(const_reference x)
-{ push_back(x); }
+append(const_reference x, const ConsistencyCmp& cmp)
+{
+ // First try the special cases that x goes at the end
+ const_reference last = CycleRepresentation::back();
+ if (empty() || cmp(last, x))
+ {
+ push_back(x);
+ return;
+ }
+
+ for (iterator cur = begin(); cur != end(); ++cur)
+ if (cmp(x, *cur))
+ {
+ insert(cur, x);
+ return;
+ }
+}
template<class I, class OrderCmp, class ConsistencyCmp>
typename Cycle<I,OrderCmp,ConsistencyCmp>::const_reference
@@ -140,7 +155,7 @@
{
for (const_iterator cur = begin(); cur != end(); ++cur)
{
- out << **cur << " ";
+ out << **cur << ", ";
}
// out << "(last: " << *last << ")"; // For debugging only
return out;
@@ -158,7 +173,7 @@
Cycle<I, OrderCmp, ConsistencyCmp>::
add(const Self& c, const ConsistencyCmp& cmp)
{
- Dout(dc::cycle, "Adding cycles");
+ Dout(dc::cycle, "Adding cycles: " << *this << " + " << c);
iterator cur1 = begin();
const_iterator cur2 = c.begin();
@@ -168,12 +183,14 @@
if (cur1 == end())
{
while (cur2 != c.end())
- append(*cur2++);
+ push_back(*cur2++);
+ Dout(dc::cycle, "After addition: " << *this);
return *this;
}
// mod 2
int res = cmp.compare(*cur1, *cur2);
+ Dout(dc::cycle, "Comparison result: " << res);
if (res == 0) // *cur1 == *cur2
{
Dout(dc::cycle, "Equality");
--- a/include/debug.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/debug.h Thu Dec 21 13:32:34 2006 -0500
@@ -65,8 +65,6 @@
extern channel_ct vineyard;
extern channel_ct cycle;
extern channel_ct lsfiltration;
- extern channel_ct lsvineyard;
- extern channel_ct vertex_transpositions;
extern channel_ct orderlist;
} // namespace dc
} // namespace DEBUGCHANNELS
@@ -76,8 +74,8 @@
#define AssertMsg(TEST,MSG) \
( (TEST) ? (void)0 \
- : (std::cerr << __FILE__ " (" << __LINE__ \
- << "): Assertion failed " #TEST \
+ : (std::cerr << __FILE__ ":" << __LINE__ \
+ << ": Assertion failed " #TEST \
<< " - " << MSG << std::endl,abort()))
/*
#define AssertMsg(TEST,STRM,MSG) \
--- a/include/eventqueue.h Sat Dec 16 15:34:46 2006 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2006
- */
-
-#ifndef __EVENTQUEUE_H__
-#define __EVENTQUEUE_H__
-
-#include <map>
-
-template<class _Key, class _Event>
-class EventQueue
-{
-
- public:
- typedef _Key Key;
- typedef _Event Event;
-
- // multimap just in case
- typedef std::multimap<Key,Event> QueueRepresentation;
- typedef typename QueueRepresentation::iterator iterator;
- typedef typename QueueRepresentation::const_iterator const_iterator;
- typedef std::pair<Key, Event> KeyValuePair;
-
- EventQueue() {}
-
- const_iterator top() const { assert(!empty()); return queue.begin(); }
- iterator push(Key k, Event e) { return queue.insert(std::make_pair(k,e)); }
- void pop() { assert(!empty()); queue.erase(queue.begin()); }
- void remove(iterator i) { queue.erase(i); }
-
- iterator end() { return queue.end(); }
- const_iterator end() const { return queue.end(); }
- bool empty() const { return queue.empty(); }
- size_t size() const { return queue.size(); }
-
- private:
- QueueRepresentation queue;
-};
-
-#endif // __EVENTQUEUE_H__
--- a/include/filtration.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/filtration.h Thu Dec 21 13:32:34 2006 -0500
@@ -11,6 +11,7 @@
#include "filtrationcontainer.h"
#include "filtrationsimplex.h"
+#include "vineyard.h"
#include <map>
#include <vector>
@@ -23,15 +24,14 @@
* for the simplex order stored in the filtration. Iterators remain valid
* through all the operations.
*/
-template<class Smplx, class FltrSmplx = FiltrationSimplex<Smplx> >
+template<class Smplx, class FltrSmplx = FiltrationSimplex<Smplx>, class Vnrd = Vineyard<FltrSmplx> >
class Filtration: public FltrSmplx::Container
{
public:
typedef Smplx Simplex;
typedef FltrSmplx FiltrationSimplex;
- typedef Filtration<Simplex, FiltrationSimplex> Self;
-
- public:
+ typedef Vnrd Vineyard;
+
/// \name Container Types
/// @{
/** The actual container type (which is the parent of the Filtration) */
@@ -51,28 +51,35 @@
typedef typename Trail::iterator TrailIterator;
/// @}
+ typedef Filtration<Simplex, FiltrationSimplex, Vineyard> Self;
typedef FiltrationContainer Parent;
public:
- Filtration();
- template<class OtherFltrSmplx> Filtration(const Filtration<Simplex,OtherFltrSmplx>& filtration);
- virtual ~Filtration();
+ Filtration(Vineyard* vineyard);
/// \name Core Functionality
/// @{
- /** Computes pairing of the simplices (RU decomposition) starting from bg, assuming that everything before bg has been paired */
- void pair_simplices(Index bg);
+ /// Computes RU decomposition of the simplices in [bg, end) range, assuming that everything before bg has been paired
+ void pair_simplices(Index bg, Index end);
+ bool transpose(Index i);
bool is_paired() const;
Index append(const Simplex& s); ///< Appends s to the filtration
+ Index insert(Index prior, const Simplex& s); ///< Inserts s after prior
const_Index get_index(const Simplex& s) const; /**< Returns the iterator pointing to s
(end() if s not in filtration) */
Index get_index(const Simplex& s); ///< \overload
void fill_simplex_index_map(); ///< Fills the mapping for get_index()
/// @}
+
+ /// \name Accessors
+ /// @{
+ Vineyard* vineyard() { return vineyard_; }
+ const Vineyard* vineyard() const { return vineyard_; }
+ /// @}
protected:
using Parent::swap;
-
+ bool transpose_simplices(Index i);
public:
/// \name Container Operations
@@ -84,10 +91,6 @@
std::ostream& operator<<(std::ostream& out) const;
- /** The usual rebind trick in absence of template typedefs */
- template<class OtherFltrSmplx>
- struct rebind { typedef Filtration<Simplex, OtherFltrSmplx> other; };
-
protected:
/// \name Comparator accessors (protected)
/// @{
@@ -95,20 +98,19 @@
const CyclesComparator& get_cycles_cmp() const { return cycles_cmp; }
const TrailsComparator& get_trails_cmp() const { return trails_cmp; }
/// @}
-
- /*
- template<class Cmp> void sort_simplices(const Cmp& cmp);
- */
private:
typedef std::map<Simplex, Index> SimplexMap;
- /** Initializes the cycle with the indices of the simplices in the boundary, and the trail with the index of this simplex */
+ /// Initializes the cycle with the indices of the simplices in the boundary, and the trail with the index of this simplex
void init_cycle_trail(Index j);
+ void pairing_switch(Index i, Index j);
bool paired;
SimplexMap inverse_simplices;
+ Vineyard* vineyard_;
+
CyclesComparator cycles_cmp;
TrailsComparator trails_cmp;
ConsistencyComparator consistency_cmp;
@@ -116,7 +118,7 @@
private:
/* Serialization */
friend class boost::serialization::access;
-
+
typedef std::map<const_Index, SizeType, ConsistencyComparator> IndexIntMap;
typedef std::vector<Index> IndexVector;
--- a/include/filtration.hpp Sat Dec 16 15:34:46 2006 -0500
+++ b/include/filtration.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -11,36 +11,18 @@
/* Filtration Public */
-template<class S, class FS>
-Filtration<S, FS>::
-Filtration(): paired(false)
+template<class S, class FS, class V>
+Filtration<S, FS, V>::
+Filtration(Vineyard* vnrd = 0): paired(false), vineyard_(vnrd)
{}
-template<class S, class FS>
-template<class OtherFltrSmplx>
-Filtration<S, FS>::
-Filtration(const Filtration<Simplex,OtherFltrSmplx>& filtration)
-{
- for(typename Filtration<Simplex,OtherFltrSmplx>::const_Index cur = filtration.begin(); cur != filtration.end(); ++cur)
- push_back(*cur);
-
- fill_simplex_index_map();
-
- // TODO: finish
-}
-
-template<class S, class FS>
-Filtration<S, FS>::
-~Filtration()
-{}
-
-template<class S, class FS>
+template<class S, class FS, class V>
void
-Filtration<S, FS>::
-pair_simplices(Index bg)
+Filtration<S, FS, V>::
+pair_simplices(Index bg, Index end)
{
Dout(dc::filtration, "Entered: compute_pairing");
- for (Index j = bg; j != end(); ++j)
+ for (Index j = bg; j != end; ++j)
{
Dout(dc::filtration|flush_cf|continued_cf, *j << ": ");
init_cycle_trail(j);
@@ -56,7 +38,6 @@
Dout(dc::filtration, *i << ": " << *(i->pair()));
AssertMsg(!cycles_cmp(i, j), "Simplices in the cycle must precede current simplex: " <<
"(" << *i << " in cycle of " << *j << ")");
- AssertMsg(i->value() <= j->value(), "Values in the cycle must be less than the value of the simplex");
// i is not paired, so we pair j with i
if (i->pair() == i)
@@ -71,34 +52,65 @@
// continue searching --- change the Dout to the continued mode with newlines FIXME
Dout(dc::filtration, " Adding: [" << bdry << "] + ");
Dout(dc::filtration, " [" << i->pair()->cycle() << "]");
- bdry.add(i->pair()->cycle(), consistency_cmp);
- i->pair()->trail().append(j);
+ bdry.add(i->pair()->cycle(), get_consistency_cmp());
+ i->pair()->trail().append(j, get_consistency_cmp());
Dout(dc::filtration, "After addition: " << bdry);
}
+ Dout(dc::filtration, "Finished with " << *j << ": " << *(j->pair()));
}
paired = true;
}
-template<class S, class FS>
+template<class S, class FS, class V>
bool
-Filtration<S, FS>::
+Filtration<S, FS, V>::
is_paired() const
{ return paired; }
-/* Filtration Protected */
-template<class S, class FS>
-typename Filtration<S, FS>::Index
-Filtration<S, FS>::
+/**
+ * Transposes simplices at i and i+1, and records the knee in the vineyard if there is a change in pairing.
+ * Returns true if the pairing changed.
+ */
+template<class S, class FS, class V>
+bool
+Filtration<S,FS,V>::
+transpose(Index i)
+{
+ AssertMsg(vineyard() != 0, "We must have a vineyard for transpositions");
+
+ Index i_orig = i++;
+
+ AssertMsg(i_orig->pair() != i, "Transposing simplices must not be paired");
+ bool result = transpose_simplices(i_orig);
+ AssertMsg(i_orig == boost::next(i), "Wrong indices after transposition");
+
+ if (result) vineyard()->switched(i, i_orig);
+ return result;
+}
+
+template<class S, class FS, class V>
+typename Filtration<S, FS, V>::Index
+Filtration<S, FS, V>::
append(const Simplex& s)
{
Index i = push_back(FiltrationSimplex(s));
- i->set_pair(i);
+ return i;
+}
+
+template<class S, class FS, class V>
+typename Filtration<S, FS, V>::Index
+Filtration<S, FS, V>::
+insert(Index prior, const Simplex& s)
+{
+ Index i = Parent::insert(prior, FiltrationSimplex(s));
+ paired = false;
+
return i;
}
-template<class S, class FS>
-typename Filtration<S, FS>::const_Index
-Filtration<S, FS>::
+template<class S, class FS, class V>
+typename Filtration<S, FS, V>::const_Index
+Filtration<S, FS, V>::
get_index(const Simplex& s) const
{
typename SimplexMap::const_iterator i = inverse_simplices.find(s);
@@ -108,9 +120,9 @@
return i->second;
}
-template<class S, class FS>
-typename Filtration<S, FS>::Index
-Filtration<S, FS>::
+template<class S, class FS, class V>
+typename Filtration<S, FS, V>::Index
+Filtration<S, FS, V>::
get_index(const Simplex& s)
{
typename SimplexMap::const_iterator i = inverse_simplices.find(s);
@@ -120,41 +132,18 @@
return i->second;
}
-template<class S, class FS>
+template<class S, class FS, class V>
void
-Filtration<S, FS>::
+Filtration<S, FS, V>::
fill_simplex_index_map()
{
for (Index i = begin(); i != end(); ++i)
inverse_simplices[*i] = i;
}
-
-/* FIXME
-template<class S>
-template<class Cmp>
-void Filtration<S>::sort_simplices(const Cmp& cmp)
-{
- // Far from efficient, but should work
- typedef std::multiset<Simplex, Cmp> SimplexSet;
- SimplexSet ordered_simplices(cmp);
- for (Index i = begin(); i != end(); ++i)
- ordered_simplices.insert(*i);
-
- simplices.clear();
- inverse_simplices.clear();
- for (typename SimplexSet::const_iterator cur = ordered_simplices.begin();
- cur != ordered_simplices.end();
- ++cur)
- {
- append(*cur);
- }
-}
-*/
-
-template<class S, class FS>
+template<class S, class FS, class V>
std::ostream&
-Filtration<S, FS>::
+Filtration<S, FS, V>::
operator<<(std::ostream& out) const
{
out << "Pairing: " << std::endl;
@@ -168,31 +157,231 @@
return out;
}
+
+/* Filtration Protected */
+/// Transposes simplices at i and i+1. Returns true if the pairing switched.
+template<class S, class FS, class V>
+bool
+Filtration<S,FS,V>::
+transpose_simplices(Index i)
+{
+ AssertMsg(is_paired(), "Pairing must be computed before transpositions");
+ counters.inc("SimplexTransposition");
+
+ Index i_prev = i++;
+
+ if (i_prev->dimension() != i->dimension())
+ {
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Different dimension");
+ counters.inc("Case DiffDim");
+ return false;
+ }
+
+ bool si = i_prev->sign(), sii = i->sign();
+ if (si && sii)
+ {
+ Dout(dc::transpositions, "Trail prev: " << i_prev->trail());
+
+ // Case 1
+ TrailIterator i_prev_second = i_prev->trail().get_second(Filtration::get_trails_cmp());
+ if (*i_prev_second == i)
+ {
+ Dout(dc::transpositions, "Case 1, U[i,i+1] = 1");
+ i_prev->trail().erase(i_prev_second);
+ }
+
+ Index k = i_prev->pair();
+ Index l = i->pair();
+
+ // Explicit treatment of unpaired simplex
+ if (l == i)
+ {
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 1.2 --- unpaired");
+ Dout(dc::transpositions, *i_prev);
+ counters.inc("Case 1.2");
+ return false;
+ } else if (k == i_prev)
+ {
+ if (*(l->cycle().get_second(Filtration::get_cycles_cmp())) != i_prev)
+ {
+ // Case 1.2
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 1.2 --- unpaired");
+ Dout(dc::transpositions, *i_prev);
+ counters.inc("Case 1.2");
+ return false;
+ } else
+ {
+ // Case 1.1.2 --- special version (plain swap, but pairing switches)
+ swap(i_prev, i);
+ pairing_switch(i_prev, i);
+ Dout(dc::transpositions, "Case 1.1.2 --- unpaired");
+ Dout(dc::transpositions, *i_prev);
+ counters.inc("Case 1.1.2");
+ return true;
+ }
+ }
+
+ Dout(dc::transpositions, "l cycle: " << l->cycle());
+ if (*(l->cycle().get_second(Filtration::get_cycles_cmp())) != i_prev)
+ {
+ // Case 1.2
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 1.2");
+ counters.inc("Case 1.2");
+ return false;
+ } else
+ {
+ // Case 1.1
+ if (trails_cmp(k,l))
+ {
+ // Case 1.1.1
+ swap(i_prev, i);
+ l->cycle().add(k->cycle(), Filtration::get_consistency_cmp()); // Add column k to l
+ k->trail().add(l->trail(), Filtration::get_consistency_cmp()); // Add row l to k
+ Dout(dc::transpositions, "Case 1.1.1");
+ counters.inc("Case 1.1.1");
+ return false;
+ } else
+ {
+ // Case 1.1.2
+ swap(i_prev, i);
+ k->cycle().add(l->cycle(), Filtration::get_consistency_cmp()); // Add column l to k
+ l->trail().add(k->trail(), Filtration::get_consistency_cmp()); // Add row k to l
+ pairing_switch(i_prev, i);
+ Dout(dc::transpositions, "Case 1.1.2");
+ counters.inc("Case 1.1.2");
+ return true;
+ }
+ }
+ } else if (!si && !sii)
+ {
+ // Case 2
+ if (*(i_prev->trail().get_second(Filtration::get_trails_cmp())) != i)
+ {
+ // Case 2.2
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 2.2");
+ counters.inc("Case 2.2");
+ return false;
+ } else
+ {
+ // Case 2.1
+ Index low_i = i_prev->pair();
+ Index low_ii = i->pair();
+ i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
+ i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
+ swap(i_prev, i);
+ if (Filtration::get_trails_cmp()(low_ii, low_i))
+ {
+ // Case 2.1.2
+ i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i to i_prev (after transposition)
+ i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
+ pairing_switch(i_prev, i);
+ Dout(dc::transpositions, "Case 2.1.2");
+ counters.inc("Case 2.1.2");
+ return true;
+ }
+
+ // Case 2.1.1
+ Dout(dc::transpositions, "Case 2.1.1");
+ counters.inc("Case 2.1.1");
+ return false;
+ }
+ } else if (!si && sii)
+ {
+ // Case 3
+ if (*(i_prev->trail().get_second(Filtration::get_trails_cmp())) != i)
+ {
+ //AssertMsg(pair(i)->cycle_get_second(cycles_cmp) != i, dc::transpositions, "Problem in Case 3");
+ // Case 3.2
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 3.2");
+ counters.inc("Case 3.2");
+ return false;
+ } else
+ {
+ // Case 3.1
+ i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
+ i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
+ swap(i_prev, i);
+ i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i (after transposition)
+ i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
+ pairing_switch(i_prev, i);
+ Dout(dc::transpositions, "Case 3.1");
+ counters.inc("Case 3.1");
+ return true;
+ }
+ } else if (si && !sii)
+ {
+ // Case 4
+ TrailIterator i_prev_second = i_prev->trail().get_second(Filtration::get_trails_cmp());
+ if (*i_prev_second == i)
+ {
+ Dout(dc::transpositions, "Case 4, U[i,i+1] = 1");
+ i_prev->trail().erase(i_prev_second);
+ }
+ swap(i_prev, i);
+ Dout(dc::transpositions, "Case 4");
+ counters.inc("Case 4");
+ return false;
+ }
+
+ return false; // to avoid compiler complaints, should never reach this point
+}
+
+
/* Filtration Private */
-template<class S, class FS>
+template<class S, class FS, class V>
void
-Filtration<S, FS>::
+Filtration<S, FS, V>::
init_cycle_trail(Index j)
{
typename Simplex::Cycle bdry = j->boundary();
for (typename Simplex::Cycle::const_iterator cur = bdry.begin(); cur != bdry.end(); ++cur)
{
+ Dout(dc::filtration, "Appending in init_cycle_trail(): " << *cur);
AssertMsg(get_index(*cur) != end(), "Non-existent simplex in the cycle");
- Dout(dc::filtration, "Appending in init_cycle_trail(): " << *cur);
- j->cycle().append(get_index(*cur));
+ j->cycle().append(get_index(*cur), get_consistency_cmp());
}
- j->cycle().sort(consistency_cmp);
-
- j->trail().append(j);
+ j->trail().append(j, get_consistency_cmp());
+ j->set_pair(j);
}
+/// Update the pairing, so that whoever was paired with i is now paired with j and vice versa.
+template<class S, class FS, class V>
+void
+Filtration<S,FS,V>::
+pairing_switch(Index i, Index j)
+{
+ Index i_pair = i->pair();
+ Index j_pair = j->pair();
+
+ if (i_pair == i)
+ j->set_pair(j);
+ else
+ {
+ j->set_pair(i_pair);
+ i_pair->set_pair(j);
+ }
+
+ if (j_pair == j)
+ i->set_pair(i);
+ else
+ {
+ i->set_pair(j_pair);
+ j_pair->set_pair(i);
+ }
+}
/* Serialization */
-template<class S, class FS>
+template<class S, class FS, class V>
template<class Archive>
void
-Filtration<S, FS>::
+Filtration<S, FS, V>::
save(Archive& ar, version_type ) const
{
ar << BOOST_SERIALIZATION_NVP(paired);
@@ -221,10 +410,10 @@
Dout(dc::filtration, count << " simplices serialized");
}
-template<class S, class FS>
+template<class S, class FS, class V>
template<class Archive>
void
-Filtration<S, FS>::
+Filtration<S, FS, V>::
load(Archive& ar, version_type )
{
Dout(dc::filtration, "Starting to read filtration");
@@ -257,9 +446,9 @@
Dout(dc::filtration, "In Filtration: simplices read");
}
-template<class S, class FS>
+template<class S, class FS, class V>
std::ostream&
-operator<<(std::ostream& out, const Filtration<S, FS>& f)
+operator<<(std::ostream& out, const Filtration<S, FS, V>& f)
{ return f.operator<<(out); }
--- a/include/filtrationsimplex.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/filtrationsimplex.h Thu Dec 21 13:32:34 2006 -0500
@@ -10,6 +10,8 @@
#include "debug.h"
#include "filtrationcontainer.h"
+#include "vineyard.h"
+#include "types.h"
#include <list>
@@ -20,32 +22,50 @@
#endif
/**
+ * Evaluator is a base class for the structures that are able to return a value
+ * given a simplex.
+ */
+template<class Smplx>
+class Evaluator
+{
+ public:
+ typedef Smplx Simplex;
+
+ virtual RealType time() { return 0; }
+ virtual RealType value(const Simplex& s) { return 0; }
+};
+
+/**
* FiltrationSimplex stores information needed for the RU-decomposition:
* cycle (column of R), trail (row of U), and pair.
*/
-template<class Smplx, class ContainerSimplex = void>
+template<class Smplx>
class FiltrationSimplex: public Smplx
{
public:
typedef Smplx Simplex;
typedef FiltrationSimplex<Simplex> Self;
- typedef FiltrationContainer<ContainerSimplex> Container;
+ typedef FiltrationContainer<Self> Container;
typedef Simplex Parent;
+ typedef Vine<Simplex> Vine;
typedef typename Container::Cycle Cycle;
typedef typename Container::Trail Trail;
typedef typename Container::Index Index;
+ typedef Evaluator<Simplex> Evaluator;
FiltrationSimplex(const Simplex& s):
- Simplex(s) {}
+ Simplex(s), vine_(0) {}
/// \name Core functionality
/// @{
void set_pair(Index pair) { pair_ = pair; }
bool sign() const { return cycles_column.empty(); }
- using Parent::dimension;
+ bool is_paired() const { return pair() != pair()->pair(); }
+ void set_vine(Vine* v) { vine_ = v; }
+ using Parent::dimension;
/// @}
@@ -56,122 +76,14 @@
const Cycle& cycle() const { return cycles_column; }
const Trail& trail() const { return trails_row; }
Index pair() const { return pair_; }
- /// @}
-
- private:
- Cycle cycles_column;
- Trail trails_row;
- Index pair_;
-};
-
-/** Specialization for ContainerSimplex = void */
-template<class Smplx>
-class FiltrationSimplex<Smplx, void> : public Smplx
-{
- public:
- typedef Smplx Simplex;
- typedef FiltrationSimplex<Simplex> Self;
- typedef FiltrationContainer<Self> Container;
- typedef Simplex Parent;
-
- typedef typename Container::Cycle Cycle;
- typedef typename Container::Trail Trail;
- typedef typename Container::Index Index;
-
-
- FiltrationSimplex(const Simplex& s):
- Simplex(s) {}
-
-
- /// \name Core functionality
- /// @{
- void set_pair(Index pair) { pair_ = pair; }
- bool sign() const { return cycles_column.empty(); }
- using Parent::dimension;
- /// @}
-
-
- /// \name Accessors
- /// @{
- Cycle& cycle() { return cycles_column; }
- Trail& trail() { return trails_row; }
- const Cycle& cycle() const { return cycles_column; }
- const Trail& trail() const { return trails_row; }
- Index pair() const { return pair_; }
+ Vine* vine() const { return vine_; }
/// @}
private:
Cycle cycles_column;
Trail trails_row;
Index pair_;
+ Vine* vine_;
};
-#if 0 // FIXME
-template<class S>
-class Filtration<S>::FiltrationSimplexSerialization: public Simplex
-{
- public:
- typedef std::list<IntegerIndex> IntegerIndexList;
-
- // Default constructor for serialization
- FiltrationSimplexSerialization() {}
-
- FiltrationSimplexSerialization(const FiltrationSimplex& fs, const IndexIntMap& im):
- Simplex(fs)
- {
- pairing = im.find(fs.pair())->second;
- for (typename FiltrationCycle::const_iterator cur = fs.cycle().begin();
- cur != fs.cycle().end();
- ++cur)
- { cycle.push_back(im.find(*cur)->second); }
-
- for (typename FiltrationTrail::const_iterator cur = fs.trail().begin();
- cur != fs.trail().end();
- ++cur)
- { trail.push_back(im.find(*cur)->second); }
-
- vine = fs.get_vine();
- }
-
- void set_filtration_simplex(FiltrationSimplex& s, const IndexVector& index_vector) const
- {
- s = *this;
-
- s.pair_with(index_vector[pairing]);
-
- // Just in case
- s.cycles_column.clear();
- s.trails_row.clear();
-
- for (IntegerIndexList::const_iterator cur = cycle.begin(); cur != cycle.end(); ++cur)
- { s.cycles_column.append(index_vector[*cur]); }
- for (IntegerIndexList::const_iterator cur = trail.begin(); cur != trail.end(); ++cur)
- { s.trails_row.append(index_vector[*cur]); }
-
- s.set_vine(vine);
- }
-
- private:
- IntegerIndexList cycle;
- IntegerIndexList trail;
- IntegerIndex pairing;
- Vine* vine;
-
- private:
- // Serialization
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type )
- {
- ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Simplex);
-
- ar & BOOST_SERIALIZATION_NVP(cycle);
- ar & BOOST_SERIALIZATION_NVP(trail);
- ar & BOOST_SERIALIZATION_NVP(pairing);
- ar & BOOST_SERIALIZATION_NVP(vine);
- }
-};
-#endif
-
#endif // __FILTRATIONSIMPLEX_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/lowerstarfiltration.h Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,114 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2006
+ */
+
+#ifndef __LOWERSTARFILTRATION_H__
+#define __LOWERSTARFILTRATION_H__
+
+#include "sys.h"
+#include "debug.h"
+
+#include "filtration.h"
+#include "simplex.h"
+#include "consistencylist.h"
+#include <boost/utility.hpp>
+#include <list>
+#include "types.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/map.hpp>
+#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/nvp.hpp>
+
+
+template<class VI,
+ class Smplx = SimplexWithAttachment<VI>,
+ class FltrSmplx = FiltrationSimplex<Smplx>,
+ class Vnrd = Vineyard<FltrSmplx> >
+class LowerStarFiltration: public Filtration<Smplx, FltrSmplx, Vnrd>
+{
+ public:
+ // Treat VertexIndex as an iterator
+ typedef VI VertexIndex;
+ typedef Smplx Simplex;
+ typedef Filtration<Simplex> Parent;
+ typedef typename Parent::Vineyard Vineyard;
+
+ typedef typename Parent::Index Index;
+ typedef typename Parent::const_Index const_Index;
+ typedef typename Parent::Cycle Cycle;
+ typedef typename Parent::Trail Trail;
+ typedef typename Simplex::Cycle SimplexBoundaryCycle;
+
+ struct VertexDescriptor;
+ typedef ConsistencyList<VertexDescriptor> VertexOrder;
+ typedef typename VertexOrder::iterator VertexOrderIndex;
+ typedef typename VertexOrder::const_iterator const_VertexOrderIndex;
+ typedef typename VertexOrder::LessThanComparison VertexOrderComparison;
+ struct SimplexAttachmentComparison;
+
+ public:
+ template<class VertexCmp>
+ LowerStarFiltration(VertexIndex begin, VertexIndex end, const VertexCmp& cmp, Vineyard* vineyard);
+
+ using Parent::size;
+ using Parent::begin;
+ using Parent::end;
+ VertexIndex num_vertices() const { return vertex_order.size(); }
+ const VertexOrderComparison&
+ get_vertex_cmp() const { return vertex_cmp; }
+
+ Index append(const Simplex& s);
+ bool transpose_vertices(const VertexOrderIndex& voi);
+
+ protected:
+ /// Hint function: if unsure, should return true
+ virtual bool neighbors(VertexIndex v1, VertexIndex v2) const { return true; }
+
+ private:
+ bool transpose(Index i);
+ void assert_pairing(Index i);
+
+ private:
+ VertexOrder vertex_order;
+ VertexOrderComparison vertex_cmp;
+
+ /* Serialization */
+ protected:
+ LowerStarFiltration() {}
+
+ private:
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void save(Archive& ar, version_type ) const { ar << BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent); }
+
+ template<class Archive>
+ void load(Archive& ar, version_type );
+
+ BOOST_SERIALIZATION_SPLIT_MEMBER()
+};
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+struct LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::VertexDescriptor
+{
+ VertexDescriptor(VertexIndex vi, Index si):
+ vertex_index(vi), simplex_index(si)
+ {}
+
+ VertexIndex vertex_index;
+ Index simplex_index;
+};
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+struct LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::SimplexAttachmentComparison
+{
+ bool operator()(const Simplex& first, const Simplex& second) const;
+ VertexOrderComparison vertex_cmp;
+};
+
+#include "lowerstarfiltration.hpp"
+
+#endif // __LOWERSTARFILTRATION_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/lowerstarfiltration.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -0,0 +1,215 @@
+/* Implementations */
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+template<class VertexCmp>
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+LowerStarFiltration(VertexIndex begin, VertexIndex end, const VertexCmp& cmp, Vineyard* vineyard):
+ Parent(vineyard)
+{
+ // Record VertexIndexes in a temporary list
+ typedef std::list<VertexIndex> VertexIndexList;
+ VertexIndexList tmp_list;
+ while (begin != end)
+ tmp_list.push_back(begin++);
+
+ // Sort the temporary list
+ tmp_list.sort(cmp);
+
+ // Record vertex order
+ for(typename VertexIndexList::const_iterator cur = tmp_list.begin(); cur != tmp_list.end(); ++cur)
+ (*cur)->set_order(vertex_order.push_back(VertexDescriptor(*cur, Parent::append(Simplex(*cur)))));
+}
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+typename LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::Index
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+append(const Simplex& s)
+{
+ AssertMsg(s.dimension() != 0, "All vertices must have been inserted in the constructor");
+
+ // Find the highest vertex
+ typename Simplex::VertexContainer::const_iterator cur = s.vertices().begin(), max = cur++;
+ for (; cur != s.vertices().end(); ++cur)
+ if (!vertex_cmp((*cur)->get_order(), (*max)->get_order()))
+ max = cur;
+
+ Index ms = (*max)->get_order()->simplex_index; Index prior;
+ do { prior = ms++; } while (ms->dimension() <= s.dimension() && ms != Parent::end() && ms->get_attachment() == *max);
+
+ Index i = Parent::insert(prior, s);
+ i->set_attachment(*max);
+
+ return i;
+}
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+bool
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::SimplexAttachmentComparison::
+operator()(const Simplex& first, const Simplex& second) const
+{
+ int cmp = vertex_cmp.compare(first.get_attachment()->get_order(), second.get_attachment()->get_order());
+ if (cmp == 0)
+ return first.dimension() < second.dimension();
+ else
+ return cmp == -1;
+}
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+bool
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+transpose_vertices(const VertexOrderIndex& order)
+{
+ counters.inc("VertexTransposition");
+
+ if ((counters.lookup("VertexTransposition") % 1000000) == 0)
+ {
+ Dout(dc::lsfiltration, "Vertex transpositions: " << counters.lookup("VertexTransposition"));
+ Dout(dc::lsfiltration, "Simplex transpositions: " << counters.lookup("SimplexTransposition"));
+ Dout(dc::lsfiltration, "Attachment changed: " << counters.lookup("ChangedAttachment"));
+ Dout(dc::lsfiltration, "Regular disconnected: " << counters.lookup("RegularDisconnected"));
+ Dout(dc::lsfiltration, "Pairing Changed: " << counters.lookup("ChangedPairing"));
+ Dout(dc::lsfiltration, "------------------------");
+ }
+
+ Dout(dc::lsfiltration, "Transposing vertices (" << order->vertex_index << ", "
+ << boost::next(order)->vertex_index << ")");
+
+ Index i = order->simplex_index;
+ Index i_prev = boost::prior(i);
+ Index i_next = boost::next(order)->simplex_index;
+ Index i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
+ Index j = boost::next(i_next);
+
+ const VertexIndex& v_i = order->vertex_index;
+ const VertexIndex& v_i_next = boost::next(order)->vertex_index;
+ bool nbghrs = neighbors(v_i, v_i_next);
+
+ bool result = false;
+
+ // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
+ while (i_next_prev != i_prev)
+ {
+ result |= transpose(i_next_prev);
+ i_next_prev = boost::prior(i_next);
+ }
+ Dout(dc::lsfiltration, "Done moving the vertex");
+
+ // Second, move the simplices attached to it
+ Dout(dc::lsfiltration, "Moving attached simplices");
+ while (j != Parent::end() && j->get_attachment() == v_i_next)
+ {
+ Dout(dc::lsfiltration, " Considering " << *j);
+ if (nbghrs && j->contains(v_i)) // short circuit
+ {
+ counters.inc("ChangedAttachment");
+ Dout(dc::lsfiltration, " Attachment changed for " << *j);
+ j->set_attachment(v_i);
+ ++j;
+ continue;
+ }
+
+ Index j_prev = j; ++j;
+ while ((--j_prev)->get_attachment() != v_i_next) // i.e., until we have reached v_i_next
+ // (and the simplices that follow it) again
+ {
+ Dout(dc::lsfiltration, " Moving: " << *j_prev << ", " << *boost::next(j_prev));
+ AssertMsg(j_prev->get_attachment() == v_i, "Simplex preceding the one being moved must be attached to v_i");
+ result |= transpose(j_prev);
+ --j_prev;
+ }
+ }
+ Dout(dc::lsfiltration, "Done moving attached simplices");
+ vertex_order.swap(order, boost::next(order));
+
+ return result;
+}
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+bool
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+transpose(Index i)
+{
+ Index j = boost::next(i);
+
+ Dout(dc::lsfiltration, " Transposing (" << *i << ", " << *(i->pair()) << ") and ("
+ << *j << ", " << *(j->pair()) << ")");
+
+ assert_pairing(i);
+ assert_pairing(j);
+
+ bool res = Parent::transpose(i);
+ Dout(dc::lsfiltration, " " << *j << ": " << *(j->pair()) << ", " << *i << ": " << *(i->pair()));
+
+ assert_pairing(j);
+ assert_pairing(i);
+
+ return res;
+}
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+void
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+assert_pairing(Index i)
+{
+#ifndef NDEBUG
+ AssertMsg(i != Parent::end(), "Cannot assert pairing of end()");
+ if (!i->sign())
+ {
+ if (i->pair() != i->cycle().top(Parent::get_cycles_cmp()))
+ {
+ Dout(dc::lsfiltration, "i (negative): " << *i);
+ Dout(dc::lsfiltration, "pair(i): " << *(i->pair()));
+ Dout(dc::lsfiltration, "i->cycle().top(): " << *(i->cycle().top(Parent::get_cycles_cmp())));
+ DoutFatal(dc::fatal, "Pairing not matching the matrix at " << *i);
+ }
+ } else
+ {
+ if (i->pair() != i)
+ {
+ if (i->pair()->cycle().top(Parent::get_cycles_cmp()) != i)
+ {
+ Dout(dc::lsfiltration, "i (positive): " << *i);
+ Dout(dc::lsfiltration, "pair(i): " << *(i->pair()));
+ Dout(dc::lsfiltration, "pair(i)->cycle(): " << i->pair()->cycle());
+ Dout(dc::lsfiltration, "pair->cycle().top(): " << *(i->pair()->cycle().top(Parent::get_cycles_cmp())));
+ DoutFatal(dc::fatal, "Pairing not matching the matrix at " << *(i->pair()));
+ }
+ }
+ }
+#endif
+}
+
+
+template<class VI, class Smplx, class FltrSmplx, class Vnrd>
+template<class Archive>
+void
+LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
+load(Archive& ar, version_type )
+{
+/*
+ ar >> BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
+
+ // Count the number of vertices
+ VertexIndex num_vertices = 0;
+ for (Index cur = begin(); cur != end(); ++cur)
+ if (dimension(cur) == 0)
+ num_vertices++;
+
+ // Second pass to record vertex_order
+ vertex_order.resize(num_vertices);
+ inverse_vertex_order.resize(num_vertices);
+ num_vertices = 0;
+ for (Index cur = begin(); cur != end(); ++cur)
+ {
+ if (dimension(cur) == 0)
+ {
+ vertex_order[num_vertices].index = cur;
+ vertex_order[num_vertices].vertex_index = *(cur->get_vertices().begin());
+ inverse_vertex_order[vertex_order[num_vertices].vertex_index] = num_vertices;
+ ++num_vertices;
+ }
+ }
+*/
+}
+
+
--- a/include/orderlist.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/orderlist.h Thu Dec 21 13:32:34 2006 -0500
@@ -64,6 +64,8 @@
/// \name Order operations
void swap(iterator i, iterator j); ///< Exchanges the order of simplices pointed to by i and j
+ template<class BinaryPredicate>
+ void sort(BinaryPredicate cmp); ///< Sorts the elements in accordance with cmp
/// \name Container operations
/// @{
@@ -116,7 +118,19 @@
T data;
OrderType tag;
- std::ostream& operator<<(std::ostream& out) const { return out << data << ": " << tag; }
+ std::ostream& operator<<(std::ostream& out) const { return out << data << ": " << tag; }
+};
+
+template<class T, class BinaryPredicate>
+class OrderListNodeComparison
+{
+ public:
+ typedef OrderListNode<T> Node;
+ OrderListNodeComparison(BinaryPredicate cmp): cmp_(cmp) {}
+ bool operator()(const Node& a, const Node& b) const { return cmp_(a.data, b.data); }
+
+ private:
+ BinaryPredicate cmp_;
};
template<class T>
--- a/include/orderlist.hpp Sat Dec 16 15:34:46 2006 -0500
+++ b/include/orderlist.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -16,6 +16,18 @@
}
template<class T>
+template<class BinaryPredicate>
+void
+OrderList<T>::
+sort(BinaryPredicate cmp)
+{
+ Parent::sort(OrderListNodeComparison<T, BinaryPredicate>(cmp));
+ OrderType cur_order = 0;
+ for (typename Parent::iterator cur = begin(); cur != end(); ++cur)
+ cur->tag = cur_order++;
+}
+
+template<class T>
typename OrderList<T>::iterator
OrderList<T>::
push_back(const_reference x)
--- a/include/simplex.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/simplex.h Thu Dec 21 13:32:34 2006 -0500
@@ -50,7 +50,6 @@
/// \name Core
/// @{
Cycle boundary() const;
- virtual RealType value() const { return 0; }
Dimension dimension() const { return vertices_.size()-1; }
/// @}
@@ -113,11 +112,10 @@
/// \name Core
/// @{
void set_value(Value value) { val = value; }
- virtual Value value() const { return val; }
+ Value get_value() const { return val; }
/// @}
const Self& operator=(const Self& s);
-
std::ostream& operator<<(std::ostream& out) const;
private:
@@ -150,7 +148,7 @@
SimplexWithAttachment(const Parent& s):
Parent(s) {}
SimplexWithAttachment(VertexIndex vi):
- Parent(vi) {}
+ Parent(vi), attachment(vi) {}
/// @}
void set_attachment(VertexIndex v) { attachment = v; }
--- a/include/types.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/types.h Thu Dec 21 13:32:34 2006 -0500
@@ -1,7 +1,7 @@
#ifndef __TYPES_H__
#define __TYPES_H__
-#include <iostream>
+#include <limits>
/* Types */
typedef bool Sign;
@@ -11,41 +11,8 @@
typedef double RealType;
typedef unsigned int SizeType;
+static RealType Infinity = std::numeric_limits<RealType>::infinity();
+
typedef const unsigned int& version_type;
-
-/**
- * Rational number type
- */
-template<typename T>
-class Rational
-{
- public:
- typedef Rational<T> Self;
-
- Rational(T v);
- Rational(const Self& other);
-
- Self& operator/=(const Self& rhs);
-
- /// \name Comparisons
- /// Assume denominator is positive.
- /// @{
- bool operator<(const Self& rhs) const;
- bool operator<=(const Self& rhs) const;
- bool operator>(const Self& rhs) const;
- bool operator>=(const Self& rhs) const;
- /// @}
-
- Self& operator=(const Self& rhs);
- RealType to_real() const;
-
- std::ostream& operator<<(std::ostream& out) const;
-
- private:
- T numerator, denominator;
-};
-
-#include "types.hpp"
-
#endif // __TYPES_H__
--- a/include/types.hpp Sat Dec 16 15:34:46 2006 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-/* Implementations */
-
-template<typename T>
-Rational<T>::
-Rational(T v):
- numerator(v), denominator(1)
-{}
-
-template<typename T>
-Rational<T>::
-Rational(const Self& other):
- numerator(other.numerator), denominator(other.denominator)
-{}
-
-template<typename T>
-typename Rational<T>::Self&
-Rational<T>::
-operator/=(const Self& rhs)
-{
- numerator *= rhs.denominator;
- denominator *= rhs.numerator;
-
- if (denominator < 0)
- {
- numerator = -numerator;
- denominator = -denominator;
- }
-
- return *this;
-}
-
-template<typename T>
-bool
-Rational<T>::
-operator<(const Self& rhs) const
-{ return (numerator * rhs.denominator < denominator * rhs.numerator); }
-
-template<typename T>
-bool
-Rational<T>::
-operator<=(const Self& rhs) const
-{ return (numerator * rhs.denominator <= denominator * rhs.numerator); }
-
-template<typename T>
-bool
-Rational<T>::
-operator>(const Self& rhs) const
-{ return rhs < (*this); }
-
-template<typename T>
-bool
-Rational<T>::
-operator>=(const Self& rhs) const
-{ return rhs <= (*this); }
-
-template<typename T>
-typename Rational<T>::Self&
-Rational<T>::
-operator=(const Self& rhs)
-{ numerator = rhs.numerator; denominator = rhs.denominator; return *this; }
-
-template<typename T>
-std::ostream&
-Rational<T>::
-operator<<(std::ostream& out) const
-{ out << numerator << " / " << denominator << " = " << (numerator/denominator); return out; }
-
-template<typename T>
-RealType
-Rational<T>::
-to_real() const
-{ return numerator/denominator; }
-
-
-template<typename T>
-std::ostream& operator<<(std::ostream& out, const Rational<T>& r)
-{
- return r.operator<<(out);
-}
--- a/include/vineyard.h Sat Dec 16 15:34:46 2006 -0500
+++ b/include/vineyard.h Thu Dec 21 13:32:34 2006 -0500
@@ -6,10 +6,9 @@
#ifndef __VINEYARD_H__
#define __VINEYARD_H__
-#include "filtration.h"
-#include "vineyardsimplex.h"
#include "types.h"
#include <list>
+#include <string>
#include <boost/serialization/access.hpp>
#include <boost/serialization/vector.hpp>
@@ -18,70 +17,129 @@
#include <boost/serialization/is_abstract.hpp>
-template<class S> class VineyardSimplex;
+template<class Smplx> class Knee;
+template<class Smplx> class Vine;
/**
- * Vineyard class. Provides transpose() function, and keeps track of vines and knees in the process.
- * Derives from Filtration, which represents the current state of the filtration (after the performed
- * transpositions).
+ * Vineyard class. Keeps track of vines and knees. switched() is the key function called
+ * by filtration when pairing switches after a Filtration::transpose().
*/
-template<class Smplx,
- class Fltr = Filtration<Smplx>,
- class VnrdSmplx = VineyardSimplex<Smplx> >
-class Vineyard: public Fltr::template rebind<VnrdSmplx>::other
+template<class FltrSmplx>
+class Vineyard
+{
+ public:
+ typedef FltrSmplx FiltrationSimplex;
+ typedef typename FiltrationSimplex::Simplex Simplex;
+ typedef Vine<Simplex> Vine;
+ typedef Knee<Simplex> Knee;
+ typedef std::list<Vine> VineList;
+ typedef std::vector<VineList> VineListVector;
+ typedef typename FiltrationSimplex::Cycle Cycle;
+
+ typedef typename FiltrationSimplex::Index Index;
+ typedef typename FiltrationSimplex::Evaluator Evaluator;
+
+ public:
+ Vineyard(Evaluator* eval = 0):
+ evaluator(eval) {}
+
+ void start_vines(Index bg, Index end);
+ void switched(Index i, Index j);
+ void record_diagram(Index bg, Index end);
+
+ void set_evaluator(Evaluator* eval) { evaluator = eval; }
+
+ void save_edges(const std::string& filename) const;
+
+ protected:
+ typename Knee::SimplexList resolve_cycle(Index i) const;
+
+ private:
+ void record_knee(Index i);
+ void start_vine(Index i);
+
+ private:
+ VineListVector vines;
+ Evaluator* evaluator;
+};
+
+/**
+ * Knee class stores the knee in R^3 as well as the cycle that is associated with the Simplex starting from the Knee.
+ */
+template<class S>
+class Knee
{
public:
- typedef Smplx Simplex;
- typedef VnrdSmplx VineyardSimplex;
- typedef typename Fltr::template rebind<VineyardSimplex>::other Filtration;
-
- /// \name Vineyard types
- /// @{
- typedef typename VineyardSimplex::Knee Knee;
- typedef typename VineyardSimplex::Vine Vine;
- /// @}
+ typedef S Simplex;
+ typedef std::list<Simplex> SimplexList;
+
+ RealType birth;
+ RealType death;
+ RealType time;
+ SimplexList cycle;
+
+ // Default parameters for serialization
+ Knee(RealType b = 0, RealType d = 0, RealType t = 0):
+ birth(b), death(d), time(t)
+ {}
+ Knee(const Knee& other):
+ birth(other.birth), death(other.death), time(other.time)
+ {}
+
+ bool is_diagonal() const { return birth == death; }
+ bool is_infinite() const { return (death == Infinity) || (birth == Infinity); }
+ void set_cycle(const SimplexList& lst) { cycle = lst; }
- /// \name Container Types
- /// @{
- typedef typename Filtration::Index Index;
- typedef typename Filtration::const_Index const_Index;
- /// @}
-
- /// \name Cycles and Trails
- /// @{
- typedef typename Filtration::CyclesComparator CyclesComparator;
- typedef typename Filtration::Cycle Cycle;
- typedef typename Filtration::TrailIterator TrailIterator;
- /// @}
+ std::ostream& operator<<(std::ostream& out) const { return out << "(" << birth << ", "
+ << death << ", "
+ << time << ")"; }
+
+ private:
+ friend class boost::serialization::access;
-
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
+};
+
+template<class S>
+std::ostream& operator<<(std::ostream& out, const Knee<S>& k) { return k.operator<<(out); }
+
+/**
+ * Vine is a list of Knees
+ */
+template<class S>
+class Vine: public std::list<Knee<S> >
+{
public:
- Vineyard();
- Vineyard(const Vineyard& vineyard); ///< Copy-constructor
- template<class OtherFltrSmplx> Vineyard(const typename Filtration::template rebind<OtherFltrSmplx>::other& filtration);
- ///< Pseudo from-base constructor
+ typedef S Simplex;
+ typedef Knee<Simplex> Knee;
+ typedef std::list<Knee> VineRepresentation;
+ typedef typename VineRepresentation::const_iterator const_knee_iterator;
- /// \name Core functionality
- /// @{
- void pair_simplices(Index bg);
- void record_frame(Index bg);
- bool transpose(Index i);
- /// @}
+ Vine() {}
+ Vine(const Knee& k) { add(k); }
+
+ void add(RealType b, RealType d, RealType t) { push_back(Knee(b,d,t)); }
+ void add(const Knee& k) { push_back(k); }
- using Filtration::is_paired;
- using Filtration::begin;
- using Filtration::end;
- using Filtration::size;
-
+ using VineRepresentation::begin;
+ using VineRepresentation::end;
+ using VineRepresentation::front;
+ using VineRepresentation::back;
+ using VineRepresentation::size;
+ using VineRepresentation::empty;
+
protected:
- bool transpose_simplices(Index i);
- typename Knee::SimplexList resolve_cycle(Index i) const;
-
+ using VineRepresentation::push_back;
+
private:
- void start_vines(Index bg);
- void pairing_switch(Index i, Index j);
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
+
#include "vineyard.hpp"
#endif // __VINEYARD_H__
--- a/include/vineyard.hpp Sat Dec 16 15:34:46 2006 -0500
+++ b/include/vineyard.hpp Thu Dec 21 13:32:34 2006 -0500
@@ -1,279 +1,152 @@
/* Implementations */
-
-template<class S, class F, class VS>
-Vineyard<S,F,VS>::
-Vineyard()
-{
-}
-
-template<class S, class F, class VS>
-Vineyard<S,F,VS>::
-Vineyard(const Vineyard& vineyard):
- Filtration(vineyard)
-{
- // TODO: copy vines
-}
-
-template<class S, class F, class VS>
-template<class OtherFltrSmplx>
-Vineyard<S,F,VS>::
-Vineyard(const typename Filtration::template rebind<OtherFltrSmplx>::other& filtration):
- Filtration(filtration)
-{
- // TODO (if anything?)
-}
-
+#include <fstream>
+#include <sstream>
-template<class S, class F, class VS>
-void
-Vineyard<S,F,VS>::
-pair_simplices(Index bg)
+template<class FS>
+void
+Vineyard<FS>::
+start_vines(Index bg, Index end)
{
- Filtration::pair_simplices(bg);
- start_vines(bg);
- record_frame(bg);
-}
+ AssertMsg(evaluator != 0, "Cannot start vines with a null evaluator");
+ for (Index cur = bg; cur != end; ++cur)
+ {
+ if (!cur->sign()) continue;
+ Dimension dim = cur->dimension();
+
+ if (dim >= vines.size())
+ {
+ AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
+ vines.push_back(std::list<Vine>());
+ }
-/** Starts a new frame by recording the current diagram in the vineyard */
-template<class S, class F, class VS>
-void
-Vineyard<S,F,VS>::
-record_frame(Index bg)
-{
- Dout(dc::vineyard, "Entered: record_frame()");
-
- AssertMsg(is_paired(), "Pairing must be computed before a vine frame can be recorded");
- for (Index i = begin(); i != end(); ++i)
- {
- if (i->sign()) continue;
- Knee* k = i->new_frame(Knee(i->pair()->value(), i->value(), 0));
- k->set_cycle(resolve_cycle(i));
+ start_vine(cur);
+ record_knee(cur);
}
}
-/**
- * Transposes simplices at i and i+1, and records the knee in the vineyard if there is a change in pairing.
- * Returns true if the pairing changed.
- */
-template<class S, class F, class VS>
-bool
-Vineyard<S,F,VS>::
-transpose(Index i)
+template<class FS>
+void
+Vineyard<FS>::
+switched(Index i, Index j)
+{
+ Vine* i_vine = i->vine();
+ Vine* j_vine = j->vine();
+ i->set_vine(j_vine);
+ j->set_vine(i_vine);
+
+ // Since the pairing has already been updated, the following assertions should be true
+ AssertMsg(i->vine() == i->pair()->vine(), "i's vine must match the vine of its pair");
+ AssertMsg(j->vine() == j->pair()->vine(), "j's vine must match the vine of its pair");
+
+ if (!i->sign()) i = i->pair();
+ if (!j->sign()) j = j->pair();
+ record_knee(i);
+ record_knee(j);
+}
+
+template<class FS>
+void
+Vineyard<FS>::
+start_vine(Index i)
{
- bool result = transpose_simplices(i);
- if (result)
+ Dout(dc::vineyard, "Starting new vine");
+ AssertMsg(i->sign(), "Can only start vines for positive simplices");
+
+ Dimension dim = i->dimension();
+ vines[dim].push_back(Vine());
+ i->set_vine(&vines[dim].back());
+ i->pair()->set_vine(i->vine());
+}
+
+/// Records the current diagram in the vineyard
+template<class FS>
+void
+Vineyard<FS>::
+record_diagram(Index bg, Index end)
+{
+ Dout(dc::vineyard, "Entered: record_diagram()");
+ AssertMsg(evaluator != 0, "Cannot record diagram with a null evaluator");
+
+ for (Index i = bg; i != end; ++i)
{
- // Record pairing: TODO
+ AssertMsg(i->vine() != 0, "Cannot process a null vine in record_diagram");
+ if (!i->sign()) continue;
+ record_knee(i);
}
- return result;
}
-/// Transposes simplices at i and i+1. Returns true if the pairing switched.
-template<class S, class F, class VS>
-bool
-Vineyard<S,F,VS>::
-transpose_simplices(Index i)
+
+template<class FS>
+void
+Vineyard<FS>::
+save_edges(const std::string& filename) const
{
- AssertMsg(is_paired(), "Pairing must be computed before transpositions");
- counters.inc("SimplexTransposition");
+ for (int i = 0; i < vines.size(); ++i)
+ {
+ std::ostringstream os; os << i;
+ std::string fn = filename + os.str() + ".edg";
+ std::ofstream out(fn.c_str());
+ for (typename VineList::const_iterator vi = vines[i].begin(); vi != vines[i].end(); ++vi)
+ for (typename Vine::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
+ {
+ out << kiprev->birth << ' ' << kiprev->death << ' ' << kiprev->time << std::endl;
+ out << ki->birth << ' ' << ki->death << ' ' << ki->time << std::endl;
+ }
+ out.close();
+ }
+}
+
+/// Records a knee for the given simplex
+template<class FS>
+void
+Vineyard<FS>::
+record_knee(Index i)
+{
+ Dout(dc::vineyard, "Entered record_knee()");
+ AssertMsg(evaluator != 0, "Cannot record knee with a null evaluator");
+ AssertMsg(i->vine() != 0, "Cannot add a knee to a null vine");
+ AssertMsg(i->sign(), "record_knee() must be called on a positive simplex");
- Index i_prev = i++;
-
- if (i_prev->dimension() != i->dimension())
+ if (!i->is_paired())
+ i->vine()->add(evaluator->value(*i), Infinity, evaluator->time());
+ else
{
- swap(i_prev, i);
- Dout(dc::transpositions, "Different dimension");
- counters.inc("Case DiffDim");
- return false;
+ Dout(dc::vineyard, "Creating knee");
+ Knee k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
+ Dout(dc::vineyard, "Knee created: " << k);
+
+ if (!k.is_diagonal() || i->vine()->empty()) // non-diagonal k, or empty vine
+ {
+ Dout(dc::vineyard, "Extending a vine");
+ i->vine()->add(k);
+ }
+ else if (i->vine()->back().is_diagonal()) // last knee is diagonal
+ {
+ AssertMsg(i->vine()->size() == 1, "Only first knee may be diagonal for a live vine");
+ Dout(dc::vineyard, "Overwriting first diagonal knee");
+ i->vine()->back() = k;
+ } else // finish this vine
+ {
+ Dout(dc::vineyard, "Finishing a vine");
+ i->vine()->add(k);
+ start_vine(i);
+ i->vine()->add(k);
+ }
}
- bool si = i_prev->sign(), sii = i->sign();
- if (si && sii)
- {
- Dout(dc::transpositions, "Trail prev: " << i_prev->trail());
-
- // Case 1
- TrailIterator i_prev_second = i_prev->trail().get_second(Filtration::get_trails_cmp());
- if (*i_prev_second == i)
- {
- Dout(dc::transpositions, "Case 1, U[i,i+1] = 1");
- i_prev->trail().erase(i_prev_second);
- }
-
- Index k = i_prev->pair();
- Index l = i->pair();
-
- // Explicit treatment of unpaired simplex
- if (l == i)
- {
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 1.2 --- unpaired");
- Dout(dc::transpositions, *i_prev);
- counters.inc("Case 1.2");
- return false;
- } else if (k == i_prev)
- {
- if (*(l->cycle().get_second(Filtration::get_cycles_cmp())) != i_prev)
- {
- // Case 1.2
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 1.2 --- unpaired");
- Dout(dc::transpositions, *i_prev);
- counters.inc("Case 1.2");
- return false;
- } else
- {
- // Case 1.1.2 --- special version (plain swap, but pairing switches)
- swap(i_prev, i);
- pairing_switch(i_prev, i);
- Dout(dc::transpositions, "Case 1.1.2 --- unpaired");
- Dout(dc::transpositions, *i_prev);
- counters.inc("Case 1.1.2");
- return true;
- }
- }
-
- Dout(dc::transpositions, "l cycle: " << l->cycle());
- if (*(l->cycle().get_second(Filtration::get_cycles_cmp())) != i_prev)
- {
- // Case 1.2
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 1.2");
- counters.inc("Case 1.2");
- return false;
- } else
- {
- // Case 1.1
- if (trails_cmp(k,l))
- {
- // Case 1.1.1
- swap(i_prev, i);
- l->cycle().add(k->cycle(), Filtration::get_consistency_cmp()); // Add column k to l
- k->trail().add(l->trail(), Filtration::get_consistency_cmp()); // Add row l to k
- Dout(dc::transpositions, "Case 1.1.1");
- counters.inc("Case 1.1.1");
- return false;
- } else
- {
- // Case 1.1.2
- swap(i_prev, i);
- k->cycle().add(l->cycle(), Filtration::get_consistency_cmp()); // Add column l to k
- l->trail().add(k->trail(), Filtration::get_consistency_cmp()); // Add row k to l
- pairing_switch(i_prev, i);
- Dout(dc::transpositions, "Case 1.1.2");
- counters.inc("Case 1.1.2");
- return true;
- }
- }
- } else if (!si && !sii)
- {
- // Case 2
- if (*(i_prev->trail().get_second(Filtration::get_trails_cmp())) != i)
- {
- // Case 2.2
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 2.2");
- counters.inc("Case 2.2");
- return false;
- } else
- {
- // Case 2.1
- Index low_i = i_prev->pair();
- Index low_ii = i->pair();
- i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
- swap(i_prev, i);
- if (Filtration::get_trails_cmp()(low_ii, low_i))
- {
- // Case 2.1.2
- i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i to i_prev (after transposition)
- i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- pairing_switch(i_prev, i);
- Dout(dc::transpositions, "Case 2.1.2");
- counters.inc("Case 2.1.2");
- return true;
- }
-
- // Case 2.1.1
- Dout(dc::transpositions, "Case 2.1.1");
- counters.inc("Case 2.1.1");
- return false;
- }
- } else if (!si && sii)
- {
- // Case 3
- if (*(i_prev->trail().get_second(Filtration::get_trails_cmp())) != i)
- {
- //AssertMsg(pair(i)->cycle_get_second(cycles_cmp) != i, dc::transpositions, "Problem in Case 3");
- // Case 3.2
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 3.2");
- counters.inc("Case 3.2");
- return false;
- } else
- {
- // Case 3.1
- i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
- swap(i_prev, i);
- i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i (after transposition)
- i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- pairing_switch(i_prev, i);
- Dout(dc::transpositions, "Case 3.1");
- counters.inc("Case 3.1");
- return true;
- }
- } else if (si && !sii)
- {
- // Case 4
- TrailIterator i_prev_second = i_prev->trail().get_second(Filtration::get_trails_cmp());
- if (*i_prev_second == i)
- {
- Dout(dc::transpositions, "Case 4, U[i,i+1] = 1");
- i_prev->trail().erase(i_prev_second);
- }
- swap(i_prev, i);
- Dout(dc::transpositions, "Case 4");
- counters.inc("Case 4");
- return false;
- }
-
- return false; // to avoid compiler complaints, should never reach this point
+ i->vine()->back().set_cycle(resolve_cycle(i));
+ Dout(dc::vineyard, "Leaving record_knee()");
}
-template<class S, class F, class VS>
-typename Vineyard<S,F,VS>::Knee::SimplexList
-Vineyard<S,F,VS>::
+template<class FS>
+typename Vineyard<FS>::Knee::SimplexList
+Vineyard<FS>::
resolve_cycle(Index i) const
{
- Dout(dc::filtration, "Entered: resolve_cycle");
+ Dout(dc::vineyard, "Entered resolve_cycle");
const Cycle& cycle = i->cycle();
-
- const CyclesComparator& cmp = Filtration::get_cycles_cmp();
-#if 0
- // Make into canonical cycle
- bool done = false;
- while (!done)
- {
- done = true;
- for (typename Cycle::const_iterator cur = cycle.begin(); cur != cycle.end(); ++cur)
- {
- if (!((*cur)->sign()))
- continue;
-
- if (cmp(i,(*cur)->pair()))
- {
- done = false;
- cycle.add((*cur)->pair()->cycle(), get_consistency_cmp());
- break;
- }
- }
- }
- std::cout << "Canonical cycle computed" << std::endl;
-#endif
-
// Resolve simplices
typename Knee::SimplexList lst;
for (typename Cycle::const_iterator cur = cycle.begin(); cur != cycle.end(); ++cur)
@@ -281,46 +154,3 @@
return lst;
}
-
-
-/* Vineyard Private */
-/** Initializes vines in VineyardSimplices; should be called after pair_simplices() */
-template<class S, class F, class VS>
-void
-Vineyard<S,F,VS>::
-start_vines(Index bg)
-{
- for (Index i = bg; i != end(); ++i)
- {
- if (!i->sign()) continue;
- Vine* v = new Vine;
- i->set_vine(v);
- i->pair()->set_vine(v); // if i is unpaired i->pair() == i
- }
-}
-
-/// Update the pairing, so that whoever was paired with i is now paired with j and vice versa.
-template<class S, class F, class VS>
-void
-Vineyard<S,F,VS>::
-pairing_switch(Index i, Index j)
-{
- Index i_pair = i->pair();
- Index j_pair = j->pair();
-
- if (i_pair == i)
- j->set_pair(j);
- else
- {
- j->set_pair(i_pair);
- i_pair->set_pair(j);
- }
-
- if (j_pair == j)
- i->set_pair(i);
- else
- {
- i->set_pair(j_pair);
- j_pair->set_pair(i);
- }
-}
--- a/include/vineyardsimplex.h Sat Dec 16 15:34:46 2006 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,125 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2006
- */
-
-#ifndef __VINEYARDSIMPLEX_H__
-#define __VINEYARDSIMPLEX_H__
-
-#include "types.h"
-
-#include <boost/serialization/access.hpp>
-
-template<class S> class Knee;
-template<class S> class Vine;
-
-/**
- * VineyardSimplex class stores a Vine as a list of VineFrames, each of which is a list of Knees.
- */
-template<class Smplx>
-class VineyardSimplex: public FiltrationSimplex<Smplx, VineyardSimplex<Smplx> >
-{
- public:
- typedef Smplx Simplex;
- typedef VineyardSimplex Self;
- typedef FiltrationSimplex<Simplex, Self> FiltrationSimplex;
-
- typedef typename FiltrationSimplex::Container Container;
- typedef typename FiltrationSimplex::Index Index;
-
- typedef Knee<Simplex> Knee;
- typedef Vine<Simplex> Vine;
-
- VineyardSimplex(const FiltrationSimplex& fs):
- FiltrationSimplex(fs) {}
-
- Vine* vine() const { return vine_; }
- void set_vine(Vine* v) { vine_ = v; }
- Knee* add_knee(const Knee& k) { assert(vine_); return vine_->add(k); }
- Knee* new_frame(const Knee& k) { assert(vine_); return vine_->new_frame(k); }
- void swap_vine(Index i) { std::swap(vine_, i->vine_); }
-
- using FiltrationSimplex::set_pair;
-
- private:
- Vine* vine_;
-};
-
-/**
- * Knee class stores the knee in R^3 as well as the cycle that is associated with the Simplex starting from the Knee.
- */
-template<class S>
-class Knee
-{
- public:
- typedef S Simplex;
- typedef std::list<Simplex> SimplexList;
-
- RealType birth;
- RealType death;
- RealType time;
- SimplexList cycle;
-
- // Default parameters for serialization
- Knee(RealType b = 0, RealType d = 0, RealType t = 0):
- birth(b), death(d), time(t)
- {}
-
- bool is_diagonal() const { return birth == death; }
- void set_cycle(const SimplexList& lst) { cycle = lst; }
-
- private:
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
-};
-
-/**
- * VineFrame is a list of Knees.
- */
-template<class S>
-class VineFrame: public std::list<Knee<S> >
-{};
-
-
-/**
- * Vine is a list of VineFrames.
- */
-template<class S>
-class Vine: public std::list<VineFrame<S> >
-{
- public:
- typedef S Simplex;
- typedef Knee<Simplex> Knee;
- typedef VineFrame<Simplex> VineFrame;
- typedef std::list<VineFrame> VineRepresentation;
- typedef typename VineRepresentation::const_iterator const_frame_iterator;
-
- Vine() {}
-
- Knee* add(RealType b, RealType d, RealType t) { back().push_back(Knee(b,d,t)); return &(back().back()); }
- Knee* add(const Knee& k) { back().push_back(k); return &(back().back()); }
- Knee* new_frame(const Knee& k) { push_back(VineFrame()); return add(k); }
- unsigned int num_frames() const { return size(); }
-
- // TODO: Do a proper frame_iterator here, and knee_iterator
- using VineRepresentation::begin;
- using VineRepresentation::end;
- using VineRepresentation::front;
- using VineRepresentation::back;
- using VineRepresentation::size;
-
- protected:
- using VineRepresentation::push_back;
-
- private:
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
-};
-
-#include "vineyardsimplex.hpp"
-
-#endif // __VINEYARDSIMPLEX_H__
--- a/include/vineyardsimplex.hpp Sat Dec 16 15:34:46 2006 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,27 +0,0 @@
-#include <boost/serialization/vector.hpp>
-#include <boost/serialization/nvp.hpp>
-#include <boost/serialization/list.hpp>
-
-using boost::serialization::make_nvp;
-
-/* Implementations */
-template<class S>
-template<class Archive>
-void
-Knee<S>::
-serialize(Archive& ar, version_type )
-{
- ar & BOOST_SERIALIZATION_NVP(birth);
- ar & BOOST_SERIALIZATION_NVP(death);
- ar & BOOST_SERIALIZATION_NVP(time);
- ar & BOOST_SERIALIZATION_NVP(cycle);
-}
-
-template<class S>
-template<class Archive>
-void
-Vine<S>::
-serialize(Archive& ar, version_type )
-{ ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(VineRepresentation); }
-
-
--- a/src/debug.cpp Sat Dec 16 15:34:46 2006 -0500
+++ b/src/debug.cpp Thu Dec 21 13:32:34 2006 -0500
@@ -13,12 +13,10 @@
{
// Add new debug channels here.
channel_ct filtration("FILTRATION");
- channel_ct transpositions("TRANSPOSITIONS");
+ channel_ct transpositions("TRANSPOS");
channel_ct vineyard("VINEYARD");
channel_ct cycle("CYCLE");
channel_ct lsfiltration("LSFILTRATION");
- channel_ct lsvineyard("LSVINEYARD");
- channel_ct vertex_transpositions("VERTEX_TRANSPOS");
channel_ct orderlist("ORDERLIST");
}
} // namespace DEBUGCHANNELS
--- a/tests/SConscript Sat Dec 16 15:34:46 2006 -0500
+++ b/tests/SConscript Thu Dec 21 13:32:34 2006 -0500
@@ -1,6 +1,6 @@
Import('*')
-sources = ['test-orderlist.cpp', 'test-consistencylist.cpp']
+sources = ['test-orderlist.cpp', 'test-consistencylist.cpp', 'test-realtype.cpp']
for s in sources:
Default(env.Program([s] + external_sources))