Moved files around: lsvineyard.h -> include/topology, examples/grid -> examples/pl-functions
--- a/examples/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
+++ b/examples/CMakeLists.txt Sun Dec 20 10:43:00 2009 -0800
@@ -4,7 +4,7 @@
add_subdirectory (consistency)
add_subdirectory (cohomology)
add_subdirectory (fitness)
-add_subdirectory (grid)
+add_subdirectory (pl-functions)
add_subdirectory (triangle)
add_subdirectory (poincare)
add_subdirectory (rips)
--- a/examples/grid/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-set (targets
- test-grid2D
- test-grid2D-vineyard
- combustion-vineyard)
-
-if (use_dsrpdb)
- set (targets ${targets}
- pdbdistance-vineyard)
-endif (use_dsrpdb)
-
-if (CGAL_FOUND)
- # add_definitions (${CGAL_CXX_FLAGS_INIT})
- add_definitions (-frounding-math)
- include_directories (${CGAL_INCLUDE_DIRS})
-
- foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${CGAL_LIBRARY} ${dsrpdb_LIBRARY})
- endforeach (t ${targets})
-
- # Extra special for program_options
- add_executable (pl-vineyard pl-vineyard.cpp)
- target_link_libraries (pl-vineyard ${libraries} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-endif (CGAL_FOUND)
--- a/examples/grid/combustion-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,88 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005
- */
-
-#include <utilities/log.h>
-
-#include <iostream>
-#include <fstream>
-#include <algorithm>
-
-#include "grid2D.h"
-#include "lsvineyard.h"
-
-const int xsize = 600;
-const int ysize = 600;
-const int var = 0; // which variable to use out of nc of them in each record in the file
-
-template<typename T>
-void read(std::ifstream& ifs, T& var)
-{
- ifs.read(reinterpret_cast<char*>(&var), sizeof(T));
-}
-
-int main(int argc, char** argv)
-{
- if (argc < 3)
- {
- std::cout << "Usage: combustion-vineyard FRAME1 FRAME2" << std::endl;
- exit(0);
- }
-
- int size0, nc0;
- int size1, nc1;
-
- std::cout << "Reading: " << argv[1] << std::endl;
- std::ifstream ifs0(argv[1], std::ios::binary);
- std::cout << "Reading: " << argv[2] << std::endl;
- std::ifstream ifs1(argv[2], std::ios::binary);
-
- if (!ifs0 || !ifs1)
- {
- std::cout << "Could not open the frames" << std::endl;
- exit(0);
- }
-
- read(ifs0, size0); read(ifs0, nc0);
- read(ifs1, size1); read(ifs1, nc1);
-
- assert(size0 == size1); assert(nc0 == nc1);
- assert(size0 == xsize*ysize);
-
- Grid2D g0(xsize, ysize), g1(xsize, ysize);
-
- for (int y = 0; y < ysize; ++y)
- for (int x = 0; x < xsize; ++x)
- for (int d = 0; d < nc0; ++d)
- {
- float val0, val1;
- read(ifs0, val0);
- read(ifs1, val1);
- if (d == var)
- {
- g0(x,y) = val0;
- g1(x,y) = val1;
- }
- }
- std::cout << "Grids read" << std::endl;
-
- // Generate the complex, initialize the vineyard (which also computes the pairing)
- typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
-
- Grid2DVineyard::LSFiltration simplices;
- g0.complex_generator(make_push_back_functor(simplices));
- Grid2DVineyard::VertexComparison vcmp(g0);
- Grid2DVineyard::SimplexComparison scmp(vcmp);
- simplices.sort(scmp);
-
- Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
- std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
- std::cout << "Pairing computed" << std::endl;
-
- // Compute vineyard
- v.compute_vineyard(g1, true);
- std::cout << "Vineyard computed" << std::endl;
-
- v.vineyard().save_edges("combustion");
-}
--- a/examples/grid/grid2D.h Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,96 +0,0 @@
-/*
- * 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 <iostream>
-//#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 <boost/serialization/export.hpp>
-
-#include <boost/iterator/counting_iterator.hpp>
-
-#include "utilities/types.h"
-#include "topology/simplex.h"
-
-/**
- * Grid2D stores a grid
- */
-class Grid2D
-{
- public:
- typedef RealType ValueType;
- typedef ValueType result_type;
-
- typedef unsigned int CoordinateIndex;
- typedef boost::counting_iterator<CoordinateIndex> iterator;
-
- typedef Simplex<CoordinateIndex> Smplx;
- 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; }
-
- iterator begin() const { return iterator(0); }
- iterator end() const { return iterator(size()); }
-
- /* 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;
-
- template<class Functor>
- void complex_generator(const Functor& f);
-
- 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__
--- a/examples/grid/grid2D.hpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,121 +0,0 @@
-#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
-
-template<class Functor>
-void
-Grid2D::
-complex_generator(const Functor& f)
-{
- for (CoordinateIndex x = 0; x < xsize() - 1; ++x)
- for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
- {
- CoordinateIndex v(seq(x,y));
- CoordinateIndex vh(seq(x+1,y));
- CoordinateIndex vv(seq(x,y+1));
- CoordinateIndex vd(seq(x+1,y+1));
-
- Smplx sh; sh.add(v);
- f(sh);
- sh.add(vh); f(sh); // Horizontal edge
- sh.add(vd); f(sh); // "Horizontal" triangle
-
- Smplx sv; sv.add(v);
- sv.add(vv); f(sv); // Vertical edge
- sv.add(vd); f(sv); // "Vertical" triangle
-
- Smplx sd; sd.add(v);
- sd.add(vd); f(sd); // Diagonal edge
-
- if (y == ysize() - 2)
- {
- Smplx s; s.add(vv);
- s.add(vd); f(s); // Top edge
- }
- if (x == xsize() - 2)
- {
- Smplx s; s.add(vh);
- s.add(vd); f(s); // Right edge
- }
- }
-
- // Last row
- for (CoordinateIndex x = 0; x < xsize(); ++x)
- {
- std::cout << x << " " << ysize() - 1 << " " << seq(x, ysize() - 1) << std::endl;
- CoordinateIndex v(seq(x, ysize() - 1));
- Smplx s; s.add(v);
- f(s);
- }
-
- // Last column
- for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
- {
- std::cout << xsize() - 1 << " " << y << " " << seq(xsize() - 1, y) << std::endl;
- CoordinateIndex v(seq(xsize() - 1, y));
- Smplx s; s.add(v);
- f(s);
- }
-}
--- a/examples/grid/lsvineyard.h Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,245 +0,0 @@
-/**
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2009
- */
-
-#ifndef __LSVINEYARD_H__
-#define __LSVINEYARD_H__
-
-#include <iostream>
-
-#include "topology/simplex.h"
-#include "topology/dynamic-persistence.h"
-#include "topology/lowerstarfiltration.h"
-#include "topology/vineyard.h"
-
-#include <utilities/indirect.h>
-
-#include <CGAL/Kinetic/Inexact_simulation_traits.h>
-#include <CGAL/Kinetic/Sort.h>
-#include <CGAL/Kinetic/Sort_visitor_base.h>
-
-#include <list>
-
-
-template<class Vertex_, class VertexEvaluator_, class Simplex_ = Simplex<Vertex_>, class Filtration_ = Filtration<Simplex_> >
-class LSVineyard
-{
- public:
- typedef LSVineyard Self;
-
- typedef Vertex_ Vertex;
- typedef VertexEvaluator_ VertexEvaluator;
- typedef typename VertexEvaluator::result_type VertexValue;
-
- typedef Simplex_ Simplex;
- typedef Filtration_ LSFiltration;
- typedef typename LSFiltration::Index LSFIndex;
-
- class KineticVertexType;
- class KineticVertexComparison;
- typedef std::list<KineticVertexType> VertexContainer;
- typedef typename VertexContainer::iterator VertexIndex;
-
- struct AttachmentData: public VineData
- {
- void set_attachment(VertexIndex v) { attachment = v; }
- VertexIndex attachment;
- };
- typedef DynamicPersistenceTrails<AttachmentData> Persistence;
- typedef typename Persistence::OrderIndex Index;
- typedef typename Persistence::iterator iterator;
-
- typedef typename Persistence::template SimplexMap<LSFiltration>
- PFMap;
-
- class Evaluator;
- class StaticEvaluator;
- class KineticEvaluator;
- class DimensionFromIterator;
-
- typedef ThroughEvaluatorComparison<VertexEvaluator> VertexComparison;
- typedef MaxVertexComparison<Simplex, VertexComparison> SimplexComparison;
-
- class TranspositionVisitor;
- friend class TranspositionVisitor;
-
- typedef Vineyard<Index, iterator, Evaluator> Vnrd;
-
- class SortVisitor;
- typedef CGAL::Kinetic::Inexact_simulation_traits 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, VertexIndex> KeyVertexMap;
-
- public:
- template<class VertexIterator>
- LSVineyard(VertexIterator begin, VertexIterator end,
- LSFiltration& filtration,
- const VertexEvaluator& veval = VertexEvaluator());
- ~LSVineyard();
-
- void compute_vineyard(const VertexEvaluator& veval, bool explicit_events = false);
- bool transpose_vertices(VertexIndex vi);
-
- const LSFiltration& filtration() const { return filtration_; }
- const Vnrd& vineyard() const { return vineyard_; }
- const Persistence& persistence() const { return persistence_; }
- const VertexComparison& vertex_comparison() const { return vcmp_; }
- const VertexEvaluator& vertex_evaluator() const { return veval_; }
- const SimplexComparison& simplex_comparison() const { return scmp_; }
-
- VertexValue vertex_value(const Vertex& v) const { return veval_(v); }
- VertexValue simplex_value(const Simplex& s) const { return vertex_value(*std::max_element(s.vertices().begin(), s.vertices().end(), vcmp_)); }
- const Simplex& pfmap(iterator i) const { return pfmap_[i]; }
- const Simplex& pfmap(Index i) const { return pfmap_[i]; }
-
- Index index(iterator i) const { return persistence_.index(i); }
-
- public:
- // For Kinetic Sort
- void swap(Key a, Key b);
-
- private:
- void change_evaluator(Evaluator* eval);
- void set_attachment(iterator i, VertexIndex vi) { persistence_.modifier()(i, boost::bind(&AttachmentData::set_attachment, bl::_1, vi)); }
- void transpose_filtration(iterator i) { filtration_.transpose(filtration_.begin() + (i - persistence_.begin())); }
-
- private:
- VertexContainer vertices_;
-
- VertexEvaluator veval_;
- VertexComparison vcmp_;
- SimplexComparison scmp_;
-
- LSFiltration& filtration_;
- Persistence persistence_;
- PFMap pfmap_;
-
- Vnrd vineyard_;
- Evaluator* evaluator_;
- unsigned time_count_;
-
- KeyVertexMap kinetic_map_;
-
-#if 0
- private:
- // Serialization
- friend class boost::serialization::access;
-
- LSVineyard() {}
-
- 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(LSVineyard)
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::KineticVertexType
-{
- public:
- KineticVertexType(const Vertex& v):
- vertex_(v) {}
-
- Key kinetic_key() const { return key_; }
- void set_kinetic_key(Key k) { key_ = k; }
-
- Vertex vertex() const { return vertex_; }
- void set_vertex(Vertex v) { vertex_ = v; }
-
- iterator simplex_index() const { return simplex_index_; }
- void set_simplex_index(iterator i) { simplex_index_ = i; }
-
- private:
- Key key_;
- Vertex vertex_;
- iterator simplex_index_;
-};
-
-template<class V, class VE, class S, class C>
-std::ostream&
-operator<<(std::ostream& out, const typename LSVineyard<V,VE,S,C>::VertexIndex& vi)
-{ return out << vi->vertex(); }
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::KineticVertexComparison: public std::binary_function<const KineticVertexType&, const KineticVertexType&, bool>
-{
- public:
- KineticVertexComparison(const VertexComparison& vcmp):
- vcmp_(vcmp) {}
-
- bool operator()(const KineticVertexType& v1, const KineticVertexType& v2) const
- { return vcmp_(v1.vertex(), v2.vertex()); }
-
- private:
- VertexComparison vcmp_;
-};
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::TranspositionVisitor: public Persistence::TranspositionVisitor
-{
- public:
- typedef typename Persistence::TranspositionVisitor Parent;
- typedef typename LSVineyard<V,VE,S,C>::iterator iterator;
- typedef typename LSVineyard<V,VE,S,C>::Index Index;
-
- TranspositionVisitor(LSVineyard& v): lsvineyard_(v) {}
-
- void transpose(iterator i) { lsvineyard_.transpose_filtration(i); }
- void switched(iterator i, SwitchType type) { lsvineyard_.vineyard_.switched(index(i), index(boost::next(i))); }
-
- private:
- Index index(iterator i) { return lsvineyard_.index(i); }
-
- LSVineyard& lsvineyard_;
-};
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::Evaluator: public std::unary_function<Index, RealType>
-{
- public:
- virtual RealType time() const =0;
- virtual RealType operator()(Index i) const =0;
- virtual Dimension dimension(Index i) const =0;
- virtual RealType operator()(iterator i) const { return operator()(&*i); }
- virtual Dimension dimension(iterator i) const { return dimension(&*i); }
-};
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
-{
- public:
- SortVisitor(LSVineyard& v):
- vineyard_(v) {}
-
- template<class Vertex_handle>
- void pre_swap(Vertex_handle a, Vertex_handle b) const { vineyard_.swap(*a,*b); }
-
- private:
- LSVineyard& vineyard_;
-};
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::DimensionFromIterator: std::unary_function<iterator, Dimension>
-{
- public:
- DimensionFromIterator(const PFMap& pfmap): pfmap_(pfmap) {}
-
- Dimension operator()(iterator i) const { return pfmap_[i].dimension(); }
-
- private:
- const PFMap& pfmap_;
-};
-
-#include "lsvineyard.hpp"
-
-#endif // __LSVINEYARD_H__
--- a/examples/grid/lsvineyard.hpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,254 +0,0 @@
-#include <utilities/log.h>
-
-#ifdef LOGGING
-static rlog::RLogChannel* rlLSVineyard = DEF_CHANNEL("lsvineyard/info", rlog::Log_Debug);
-static rlog::RLogChannel* rlLSVineyardDebug = DEF_CHANNEL("lsvineyard/debug", rlog::Log_Debug);
-#endif // LOGGING
-
-#ifdef COUNTERS
-static Counter* cVertexTransposition = GetCounter("lsfiltration/transposition"); // counts number of vertex transpositions
-static Counter* cAttachment = GetCounter("lsfiltration/attachment"); // counts the number of attachment changes
-#endif
-
-
-template<class V, class VE, class S, class F>
-template<class VertexIterator>
-LSVineyard<V,VE,S,F>::
-LSVineyard(VertexIterator begin, VertexIterator end,
- LSFiltration& filtration,
- const VertexEvaluator& veval):
- vertices_(begin, end), filtration_(filtration),
- persistence_(filtration_),
- veval_(veval), vcmp_(veval_), scmp_(vcmp_),
- pfmap_(persistence_.make_simplex_map(filtration_)),
- time_count_(0)
-{
- // TODO: LSVineyard really should sort the filtration_ itself
- // filtration_.sort(scmp_);
- // persistence_.init(filtration_);
-
- rLog(rlLSVineyardDebug, "Initializing LSVineyard");
- persistence_.pair_simplices();
- rLog(rlLSVineyardDebug, "Simplices paired");
-
- vertices_.sort(KineticVertexComparison(vcmp_)); // sort vertices w.r.t. vcmp_
-#if LOGGING
- rLog(rlLSVineyardDebug, "Vertex order:");
- for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
- rLog(rlLSVineyardDebug, " %d", cur->vertex());
-#endif
-
- // Initialize simplex_index() and attachment
- VertexIndex vi = boost::prior(vertices_.begin());
- for (iterator i = persistence_.begin(); i != persistence_.end(); ++i)
- {
- const Simplex& s = pfmap_[i];
- if (s.dimension() == 0)
- {
- ++vi;
- AssertMsg(s.vertices().front() == vi->vertex(), "In constructor, simplices and vertices must match.");
- vi->set_simplex_index(i);
- }
- set_attachment(i, vi);
- rLog(rlLSVineyardDebug, "%s attached to %d", tostring(pfmap(i)).c_str(), vi->vertex());
- }
-
- evaluator_ = new StaticEvaluator(*this, time_count_);
- vineyard_.set_evaluator(evaluator_);
- vineyard_.start_vines(persistence_.begin(), persistence_.end());
-}
-
-template<class V, class VE, class S, class F>
-LSVineyard<V,VE,S,F>::
-~LSVineyard()
-{
- delete evaluator_;
-}
-
-template<class V, class VE, class S, class F_>
-void
-LSVineyard<V,VE,S,F_>::
-compute_vineyard(const VertexEvaluator& veval, bool explicit_events)
-{
- typedef Traits::Kinetic_kernel::Point_1 Point;
- typedef Traits::Kinetic_kernel::Function_kernel::Construct_function CF;
- typedef Traits::Kinetic_kernel::Motion_function F;
-
- Traits tr(0,1);
- Simulator::Handle sp = tr.simulator_handle();
- ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
- Sort sort(tr, SortVisitor(*this));
-
- // Setup the (linear) trajectories
- rLog(rlLSVineyard, "Setting up trajectories");
- CF cf;
- kinetic_map_.clear();
- for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
- {
- VertexValue val0 = veval_(cur->vertex());
- VertexValue val1 = veval(cur->vertex());
- rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
- F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
- Point p(x);
- cur->set_kinetic_key(apt->insert(p));
- kinetic_map_[cur->kinetic_key()] = cur;
- rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
- }
-
- // Process all the events (compute the vineyard in the process)
- change_evaluator(new KineticEvaluator(*this, sp, apt, time_count_));
- if (explicit_events)
- {
- while (sp->next_event_time() < 1)
- {
- rLog(rlLSVineyardDebug, "Next event time: %f", sp->next_event_time());
- sp->set_current_event_number(sp->current_event_number() + 1);
- rLog(rlLSVineyardDebug, "Processed event");
- }
- } else
- sp->set_current_time(1.0);
- rLog(rlLSVineyard, "Processed %d events", sp->current_event_number());
-
- veval_ = veval;
- change_evaluator(new StaticEvaluator(*this, ++time_count_));
- vineyard_.record_diagram(persistence_.begin(), persistence_.end());
-}
-
-template<class V, class VE, class S, class F>
-void
-LSVineyard<V,VE,S,F>::
-swap(Key a, Key b)
-{
- rLog(rlLSVineyardDebug, "Entered swap");
- VertexIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
- AssertMsg(vcmp_(ao->vertex(), bo->vertex()), "In swap(a,b), a must precede b");
- transpose_vertices(ao);
- // AssertMsg(vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), b must precede a after the transposition");
-}
-
-template<class V, class VE, class S, class F>
-void
-LSVineyard<V,VE,S,F>::
-change_evaluator(Evaluator* eval)
-{
- AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
-
- delete evaluator_;
- evaluator_ = eval;
- vineyard_.set_evaluator(evaluator_);
-}
-
-#include <boost/function.hpp>
-#include <boost/bind.hpp>
-#include <boost/lambda/lambda.hpp>
-namespace bl = boost::lambda;
-
-template<class V, class VE, class S, class F>
-bool
-LSVineyard<V,VE,S,F>::
-transpose_vertices(VertexIndex vi)
-{
- Count(cVertexTransposition);
- rLog(rlLSVineyard, "Transposing vertices (%d, %d)", vi->vertex(), boost::next(vi)->vertex());
-
- DimensionFromIterator dim(pfmap_);
- TranspositionVisitor visitor(*this);
-
- iterator i = vi->simplex_index();
- iterator i_prev = boost::prior(i);
- iterator i_next = boost::next(vi)->simplex_index();
- iterator i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
- iterator j = boost::next(i_next);
-
- VertexIndex vi_next = boost::next(vi);
- const Vertex& v = vi->vertex();
-
- bool result = false; // has a switch in pairing occurred
-
- // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
- rLog(rlLSVineyardDebug, "Starting to move the vertex");
- while (i_next_prev != i_prev)
- {
- rLog(rlLSVineyardDebug, " Transposing %s %s", tostring(pfmap(i_next_prev)).c_str(),
- tostring(pfmap(boost::next(i_next_prev))).c_str());
- result |= persistence_.transpose(i_next_prev, dim, visitor);
- i_next_prev = boost::prior(i_next);
- }
- rLog(rlLSVineyardDebug, "Done moving the vertex");
-
- // Second, move the simplices attached to it
- rLog(rlLSVineyardDebug, "Moving attached simplices");
- while (j != persistence_.end() && j->attachment == vi_next)
- {
- rLog(rlLSVineyardDebug, " Considering %s", tostring(pfmap(j)).c_str());
- if (pfmap(j).contains(v)) // j becomes attached to v and does not move
- {
- Count(cAttachment);
- rLog(rlLSVineyardDebug, " Attachment changed for ", tostring(pfmap(j)).c_str());
- set_attachment(j, vi);
- ++j;
- continue;
- }
-
- iterator j_prev = j; ++j;
- while ((--j_prev)->attachment != vi_next) // i.e., until we have reached vi_next (and the simplices that follow it) again
- {
- rLog(rlLSVineyardDebug, " Moving: %s, %s",
- tostring(pfmap(j_prev)).c_str(),
- tostring(pfmap(boost::next(j_prev))).c_str());
- AssertMsg(j_prev->attachment == vi, "Simplex preceding the one being moved must be attached to v");
- result |= persistence_.transpose(j_prev, dim, visitor);
- --j_prev;
- }
- }
- rLog(rlLSVineyard, "Done moving attached simplices");
- vertices_.splice(vi, vertices_, vi_next); // swap vi and vi_next
-
- return result;
-}
-
-
-/* Evaluators */
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::StaticEvaluator: public Evaluator
-{
- public:
- StaticEvaluator(const LSVineyard& v, RealType time):
- time_(time), vineyard_(v) {}
-
- virtual RealType time() const { return time_; }
- virtual RealType operator()(Index i) const { return vineyard_.simplex_value(vineyard_.pfmap(i)); }
- virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
-
- private:
- RealType time_;
- const LSVineyard& vineyard_;
-};
-
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::KineticEvaluator: public Evaluator
-{
- public:
- KineticEvaluator(const LSVineyard& v, Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset):
- vineyard_(v), sp_(sp), apt_(apt), time_offset_(time_offset) {}
-
- virtual RealType time() const { return time_offset_ + CGAL::to_double(get_time()); }
- virtual RealType operator()(Index i) const
- {
- rLog(rlLSVineyard, "%s (attached to %d): %s(%f) = %f", tostring(vineyard_.pfmap(i)).c_str(),
- i->attachment->vertex(),
- tostring(apt_->at(i->attachment->kinetic_key()).x()).c_str(),
- get_time(),
- CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())));
- return CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time()));
- }
- virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
-
- private:
- Simulator::Time get_time() const { return sp_->current_time(); }
-
- const LSVineyard& vineyard_;
- Simulator::Handle sp_;
- ActivePointsTable::Handle apt_;
- RealType time_offset_;
-};
--- a/examples/grid/pdbdistance-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,90 +0,0 @@
-//#include <boost/archive/binary_oarchive.hpp>
-#include "utilities/log.h"
-
-#include "pdbdistance.h"
-#include "lsvineyard.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 LOGGING
- rlog::RLogInit(argc, argv);
-
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
-#endif
-
- if (argc < 5)
- {
- std::cout << "Usage: pdbdistance FILENAME LASTFRAME LASTSUBFRAME OUTFILENAME [CAs_ONLY]" << 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 << " CAs_ONLY - only use alpha carbons [1 = true, 0 = false, default: 1]" << 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];
- bool cas_only = true;
- if (argc > 5)
- std::istringstream(argv[5]) >> cas_only;
-
- // Compute initial filtration
- int f = 0; int sf = 0;
- std::ifstream in(frame_filename(infilename, f, sf++).c_str());
- PDBDistanceGrid ginit(in, cas_only);
- in.close();
-
- typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
-
- Grid2DVineyard::LSFiltration simplices;
- ginit.complex_generator(make_push_back_functor(simplices));
- Grid2DVineyard::VertexComparison vcmp(ginit);
- Grid2DVineyard::SimplexComparison scmp(vcmp);
- simplices.sort(scmp);
- std::cout << "Complex generated, size: " << simplices.size() << std::endl;
-
- Grid2DVineyard v(ginit.begin(), ginit.end(), simplices, ginit);
- std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
- 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(PDBDistanceGrid(in, cas_only));
- 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
-}
--- a/examples/grid/pdbdistance.h Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,87 +0,0 @@
-/*
- * 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 "utilities/types.h"
-#include "grid2D.h"
-
-#include <boost/serialization/export.hpp>
-
-
-class PDBDistanceGrid: public Grid2D
-{
- public:
- PDBDistanceGrid()
- {}
-
- PDBDistanceGrid(std::istream& in, bool ca_only = true)
- {
- load_stream(in, ca_only);
- }
-
- void load_stream(std::istream& in, bool ca_only = true)
- {
- dsrpdb::Protein p(in);
- typedef std::vector<dsrpdb::Point> PointVector;
- PointVector coordinates;
- if (ca_only)
- {
- PointVector v(ca_coordinates_begin(p), ca_coordinates_end(p));
- coordinates.swap(v);
- }
- else
- {
- PointVector v(backbone_coordinates_begin(p), backbone_coordinates_end(p));
- coordinates.swap(v);
- }
-
- std::cout << "Coordinatess created, size: " << coordinates.size() << std::endl;
-
- Grid2D::change_dimensions(coordinates.size(), coordinates.size());
- for (Grid2D::CoordinateIndex i = 0; i < coordinates.size(); ++i)
- for (Grid2D::CoordinateIndex j = 0; j < coordinates.size(); ++j)
- {
- if (i < j)
- Grid2D::operator()(i,j) = distance(coordinates[i], coordinates[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__
--- a/examples/grid/pl-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,169 +0,0 @@
-#include <utilities/log.h>
-
-#include <iostream>
-#include <fstream>
-#include <string>
-#include <vector>
-
-#include "lsvineyard.h"
-
-#include <boost/program_options.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/function.hpp>
-#include <boost/bind.hpp>
-#include <boost/lambda/lambda.hpp>
-namespace bl = boost::lambda;
-
-
-typedef float VertexValue;
-typedef unsigned Vertex;
-typedef std::vector<VertexValue> VertexVector;
-struct SubscriptFunctor: public std::unary_function<Vertex, VertexValue>
-{
- SubscriptFunctor(const VertexVector& v): vec(&v) {}
- float operator()(Vertex i) const { return (*vec)[i]; }
- SubscriptFunctor& operator=(const SubscriptFunctor& other) { vec = other.vec; return *this; }
- const VertexVector* vec;
-};
-typedef SubscriptFunctor VertexEvaluator;
-typedef std::vector<VertexVector> VertexVectorVector;
-typedef LSVineyard<Vertex, VertexEvaluator> PLVineyard;
-typedef PLVineyard::Simplex Smplx; // gotta start using namespaces
-
-void program_options(int argc, char* argv[], std::string& complex_fn,
- std::string& values_fn,
- std::string& output_prefix,
- bool& skip_infinite_vines,
- bool& save_vines,
- bool& explicit_events);
-void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices);
-void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices);
-
-
-int main(int argc, char** argv)
-{
-#ifdef LOGGING
- rlog::RLogInit(argc, argv);
-#endif
-
- std::string complex_fn, values_fn, output_prefix;
- bool skip_infinite_vines = false, explicit_events = false, save_vines = false;
- program_options(argc, argv, complex_fn, values_fn, output_prefix, skip_infinite_vines, save_vines, explicit_events);
-
-
- // Read in the complex
- PLVineyard::LSFiltration simplices;
- read_simplices(complex_fn, simplices);
- std::cout << "Complex read, size: " << simplices.size() << std::endl;
-
- // Read in vertex values
- VertexVectorVector vertices;
- read_vertices(values_fn, vertices);
-
- // Setup the vineyard
- VertexEvaluator veval(vertices[0]);
- PLVineyard::VertexComparison vcmp(veval);
- PLVineyard::SimplexComparison scmp(vcmp);
- simplices.sort(scmp);
- PLVineyard v(boost::counting_iterator<Vertex>(0),
- boost::counting_iterator<Vertex>(vertices[0].size()),
- simplices, veval);
- std::cout << "Pairing computed" << std::endl;
-
- // Compute vineyard
- for (size_t i = 1; i < vertices.size(); ++i)
- {
- veval = VertexEvaluator(vertices[i]);
- v.compute_vineyard(veval, explicit_events);
- std::cout << "Processed frame: " << i << std::endl;
- }
- std::cout << "Vineyard computed" << std::endl;
-
- if (save_vines)
- v.vineyard().save_vines(output_prefix, skip_infinite_vines);
- else
- v.vineyard().save_edges(output_prefix, skip_infinite_vines);
-}
-
-
-void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices)
-{
- std::ifstream in(complex_fn.c_str());
- std::string line;
- while (std::getline(in, line))
- {
- std::istringstream strin(line);
- simplices.push_back(Smplx(std::istream_iterator<Vertex>(strin), std::istream_iterator<Vertex>()));
- }
- std::cout << "Simplices read:" << std::endl;
- std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Smplx>(std::cout, "\n"));
-}
-
-void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices)
-{
- std::ifstream in(vertex_fn.c_str());
- std::string line;
- while (std::getline(in, line))
- {
- std::istringstream strin(line);
- vertices.push_back(VertexVector(std::istream_iterator<VertexValue>(strin), std::istream_iterator<VertexValue>()));
- }
- std::cout << "Vertex values read:" << std::endl;
- for (size_t i = 0; i < vertices.size(); ++i)
- {
- std::copy(vertices[i].begin(), vertices[i].end(), std::ostream_iterator<VertexValue>(std::cout, " "));
- std::cout << std::endl;
- }
-}
-
-void program_options(int argc, char* argv[], std::string& complex_fn,
- std::string& values_fn,
- std::string& output_prefix,
- bool& skip_infinite_vines,
- bool& save_vines,
- bool& explicit_events)
-{
- namespace po = boost::program_options;
-
- po::options_description hidden("Hidden options");
- hidden.add_options()
- ("complex-file", po::value<std::string>(&complex_fn), "file listing the simplices of the complex")
- ("values-file", po::value<std::string>(&values_fn), "file listing the values at the vertices")
- ("output-prefix", po::value<std::string>(&output_prefix), "output prefix");
-
- po::options_description visible("Allowed options", 100);
- visible.add_options()
- ("skip-infinite,s", po::bool_switch(&skip_infinite_vines), "skip infinite vines in the output")
- ("explicit-events,e", po::bool_switch(&explicit_events), "process kinetic sort events one by one")
- ("save-vines,v", po::bool_switch(&save_vines), "save vines instead of edges")
- ("help,h", "produce help message");
-#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("complex-file", 1);
- pos.add("values-file", 1);
- pos.add("output-prefix", 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("complex-file") || !vm.count("values-file") || !vm.count("output-prefix"))
- {
- std::cout << "Usage: " << argv[0] << " [options] complex-file values-file output-prefix" << std::endl;
- std::cout << visible << std::endl;
- std::abort();
- }
-}
--- a/examples/grid/test-grid2D-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,54 +0,0 @@
-#include <utilities/log.h>
-
-#include <iostream>
-#include <fstream>
-#include <algorithm>
-
-#include "grid2D.h"
-#include "lsvineyard.h"
-
-int main(int argc, char** argv)
-{
-#ifdef LOGGING
- rlog::RLogInit(argc, argv);
- // stdoutLog.subscribeTo(RLOG_CHANNEL("topology"));
- // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/persistence/transpositions"));
- // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
- stdoutLog.subscribeTo(RLOG_CHANNEL("lsvineyard"));
-#endif
-
- Grid2D g0(2, 2), g1(2, 2);
- g0(0,0) = 1; g0(0,1) = 2; g0(1,0) = 3; g0(1,1) = 0;
- g1(0,0) = 4; g1(0,1) = 2; g1(1,0) = 3; g1(1,1) = 5;
-
- // Generate the complex, initialize the vineyard (which also computes the pairing)
- typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
-
- Grid2DVineyard::LSFiltration simplices;
- g0.complex_generator(make_push_back_functor(simplices));
- Grid2DVineyard::VertexComparison vcmp(g0);
- Grid2DVineyard::SimplexComparison scmp(vcmp);
- simplices.sort(scmp);
- std::cout << "Complex generated, size: " << simplices.size() << std::endl;
- std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Grid2D::Smplx>(std::cout, "\n"));
-
- Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
- std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
- std::cout << "Pairing computed" << std::endl;
-
- // Simplex order before
- std::cout << "Simplex order:" << std::endl;
- for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
- std::cout << " " << v.filtration().simplex(cur) << std::endl;
-
- // Compute vineyard
- v.compute_vineyard(g1, true);
- std::cout << "Vineyard computed" << std::endl;
-
- // Simplex order after
- std::cout << "Simplex order:" << std::endl;
- for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
- std::cout << " " << v.filtration().simplex(cur) << std::endl;
-
- v.vineyard().save_edges("test-vineyard");
-}
--- a/examples/grid/test-grid2D.cpp Sun Dec 20 08:27:01 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,16 +0,0 @@
-#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;
-}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/CMakeLists.txt Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,24 @@
+set (targets
+ test-grid2D
+ test-grid2D-vineyard
+ combustion-vineyard)
+
+if (use_dsrpdb)
+ set (targets ${targets}
+ pdbdistance-vineyard)
+endif (use_dsrpdb)
+
+if (CGAL_FOUND)
+ # add_definitions (${CGAL_CXX_FLAGS_INIT})
+ add_definitions (-frounding-math)
+ include_directories (${CGAL_INCLUDE_DIRS})
+
+ foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${CGAL_LIBRARY} ${dsrpdb_LIBRARY})
+ endforeach (t ${targets})
+
+ # Extra special for program_options
+ add_executable (pl-vineyard pl-vineyard.cpp)
+ target_link_libraries (pl-vineyard ${libraries} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endif (CGAL_FOUND)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/combustion-vineyard.cpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,88 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005
+ */
+
+#include <utilities/log.h>
+
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+#include "grid2D.h"
+#include <topology/lsvineyard.h>
+
+const int xsize = 600;
+const int ysize = 600;
+const int var = 0; // which variable to use out of nc of them in each record in the file
+
+template<typename T>
+void read(std::ifstream& ifs, T& var)
+{
+ ifs.read(reinterpret_cast<char*>(&var), sizeof(T));
+}
+
+int main(int argc, char** argv)
+{
+ if (argc < 3)
+ {
+ std::cout << "Usage: combustion-vineyard FRAME1 FRAME2" << std::endl;
+ exit(0);
+ }
+
+ int size0, nc0;
+ int size1, nc1;
+
+ std::cout << "Reading: " << argv[1] << std::endl;
+ std::ifstream ifs0(argv[1], std::ios::binary);
+ std::cout << "Reading: " << argv[2] << std::endl;
+ std::ifstream ifs1(argv[2], std::ios::binary);
+
+ if (!ifs0 || !ifs1)
+ {
+ std::cout << "Could not open the frames" << std::endl;
+ exit(0);
+ }
+
+ read(ifs0, size0); read(ifs0, nc0);
+ read(ifs1, size1); read(ifs1, nc1);
+
+ assert(size0 == size1); assert(nc0 == nc1);
+ assert(size0 == xsize*ysize);
+
+ Grid2D g0(xsize, ysize), g1(xsize, ysize);
+
+ for (int y = 0; y < ysize; ++y)
+ for (int x = 0; x < xsize; ++x)
+ for (int d = 0; d < nc0; ++d)
+ {
+ float val0, val1;
+ read(ifs0, val0);
+ read(ifs1, val1);
+ if (d == var)
+ {
+ g0(x,y) = val0;
+ g1(x,y) = val1;
+ }
+ }
+ std::cout << "Grids read" << std::endl;
+
+ // Generate the complex, initialize the vineyard (which also computes the pairing)
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
+
+ Grid2DVineyard::LSFiltration simplices;
+ g0.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(g0);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+
+ Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
+ std::cout << "Pairing computed" << std::endl;
+
+ // Compute vineyard
+ v.compute_vineyard(g1, true);
+ std::cout << "Vineyard computed" << std::endl;
+
+ v.vineyard().save_edges("combustion");
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/grid2D.h Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,96 @@
+/*
+ * 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 <iostream>
+//#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 <boost/serialization/export.hpp>
+
+#include <boost/iterator/counting_iterator.hpp>
+
+#include "utilities/types.h"
+#include "topology/simplex.h"
+
+/**
+ * Grid2D stores a grid
+ */
+class Grid2D
+{
+ public:
+ typedef RealType ValueType;
+ typedef ValueType result_type;
+
+ typedef unsigned int CoordinateIndex;
+ typedef boost::counting_iterator<CoordinateIndex> iterator;
+
+ typedef Simplex<CoordinateIndex> Smplx;
+ 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; }
+
+ iterator begin() const { return iterator(0); }
+ iterator end() const { return iterator(size()); }
+
+ /* 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;
+
+ template<class Functor>
+ void complex_generator(const Functor& f);
+
+ 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/pl-functions/grid2D.hpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,121 @@
+#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
+
+template<class Functor>
+void
+Grid2D::
+complex_generator(const Functor& f)
+{
+ for (CoordinateIndex x = 0; x < xsize() - 1; ++x)
+ for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
+ {
+ CoordinateIndex v(seq(x,y));
+ CoordinateIndex vh(seq(x+1,y));
+ CoordinateIndex vv(seq(x,y+1));
+ CoordinateIndex vd(seq(x+1,y+1));
+
+ Smplx sh; sh.add(v);
+ f(sh);
+ sh.add(vh); f(sh); // Horizontal edge
+ sh.add(vd); f(sh); // "Horizontal" triangle
+
+ Smplx sv; sv.add(v);
+ sv.add(vv); f(sv); // Vertical edge
+ sv.add(vd); f(sv); // "Vertical" triangle
+
+ Smplx sd; sd.add(v);
+ sd.add(vd); f(sd); // Diagonal edge
+
+ if (y == ysize() - 2)
+ {
+ Smplx s; s.add(vv);
+ s.add(vd); f(s); // Top edge
+ }
+ if (x == xsize() - 2)
+ {
+ Smplx s; s.add(vh);
+ s.add(vd); f(s); // Right edge
+ }
+ }
+
+ // Last row
+ for (CoordinateIndex x = 0; x < xsize(); ++x)
+ {
+ std::cout << x << " " << ysize() - 1 << " " << seq(x, ysize() - 1) << std::endl;
+ CoordinateIndex v(seq(x, ysize() - 1));
+ Smplx s; s.add(v);
+ f(s);
+ }
+
+ // Last column
+ for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
+ {
+ std::cout << xsize() - 1 << " " << y << " " << seq(xsize() - 1, y) << std::endl;
+ CoordinateIndex v(seq(xsize() - 1, y));
+ Smplx s; s.add(v);
+ f(s);
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/pdbdistance-vineyard.cpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,90 @@
+//#include <boost/archive/binary_oarchive.hpp>
+#include "utilities/log.h"
+
+#include "pdbdistance.h"
+#include <topology/lsvineyard.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 LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
+#endif
+
+ if (argc < 5)
+ {
+ std::cout << "Usage: pdbdistance FILENAME LASTFRAME LASTSUBFRAME OUTFILENAME [CAs_ONLY]" << 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 << " CAs_ONLY - only use alpha carbons [1 = true, 0 = false, default: 1]" << 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];
+ bool cas_only = true;
+ if (argc > 5)
+ std::istringstream(argv[5]) >> cas_only;
+
+ // Compute initial filtration
+ int f = 0; int sf = 0;
+ std::ifstream in(frame_filename(infilename, f, sf++).c_str());
+ PDBDistanceGrid ginit(in, cas_only);
+ in.close();
+
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
+
+ Grid2DVineyard::LSFiltration simplices;
+ ginit.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(ginit);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ std::cout << "Complex generated, size: " << simplices.size() << std::endl;
+
+ Grid2DVineyard v(ginit.begin(), ginit.end(), simplices, ginit);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
+ 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(PDBDistanceGrid(in, cas_only));
+ 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/pl-functions/pdbdistance.h Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,87 @@
+/*
+ * 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 "utilities/types.h"
+#include "grid2D.h"
+
+#include <boost/serialization/export.hpp>
+
+
+class PDBDistanceGrid: public Grid2D
+{
+ public:
+ PDBDistanceGrid()
+ {}
+
+ PDBDistanceGrid(std::istream& in, bool ca_only = true)
+ {
+ load_stream(in, ca_only);
+ }
+
+ void load_stream(std::istream& in, bool ca_only = true)
+ {
+ dsrpdb::Protein p(in);
+ typedef std::vector<dsrpdb::Point> PointVector;
+ PointVector coordinates;
+ if (ca_only)
+ {
+ PointVector v(ca_coordinates_begin(p), ca_coordinates_end(p));
+ coordinates.swap(v);
+ }
+ else
+ {
+ PointVector v(backbone_coordinates_begin(p), backbone_coordinates_end(p));
+ coordinates.swap(v);
+ }
+
+ std::cout << "Coordinatess created, size: " << coordinates.size() << std::endl;
+
+ Grid2D::change_dimensions(coordinates.size(), coordinates.size());
+ for (Grid2D::CoordinateIndex i = 0; i < coordinates.size(); ++i)
+ for (Grid2D::CoordinateIndex j = 0; j < coordinates.size(); ++j)
+ {
+ if (i < j)
+ Grid2D::operator()(i,j) = distance(coordinates[i], coordinates[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/pl-functions/pl-vineyard.cpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,169 @@
+#include <utilities/log.h>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include <topology/lsvineyard.h>
+
+#include <boost/program_options.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
+
+typedef float VertexValue;
+typedef unsigned Vertex;
+typedef std::vector<VertexValue> VertexVector;
+struct SubscriptFunctor: public std::unary_function<Vertex, VertexValue>
+{
+ SubscriptFunctor(const VertexVector& v): vec(&v) {}
+ float operator()(Vertex i) const { return (*vec)[i]; }
+ SubscriptFunctor& operator=(const SubscriptFunctor& other) { vec = other.vec; return *this; }
+ const VertexVector* vec;
+};
+typedef SubscriptFunctor VertexEvaluator;
+typedef std::vector<VertexVector> VertexVectorVector;
+typedef LSVineyard<Vertex, VertexEvaluator> PLVineyard;
+typedef PLVineyard::Simplex Smplx; // gotta start using namespaces
+
+void program_options(int argc, char* argv[], std::string& complex_fn,
+ std::string& values_fn,
+ std::string& output_prefix,
+ bool& skip_infinite_vines,
+ bool& save_vines,
+ bool& explicit_events);
+void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices);
+void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices);
+
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+#endif
+
+ std::string complex_fn, values_fn, output_prefix;
+ bool skip_infinite_vines = false, explicit_events = false, save_vines = false;
+ program_options(argc, argv, complex_fn, values_fn, output_prefix, skip_infinite_vines, save_vines, explicit_events);
+
+
+ // Read in the complex
+ PLVineyard::LSFiltration simplices;
+ read_simplices(complex_fn, simplices);
+ std::cout << "Complex read, size: " << simplices.size() << std::endl;
+
+ // Read in vertex values
+ VertexVectorVector vertices;
+ read_vertices(values_fn, vertices);
+
+ // Setup the vineyard
+ VertexEvaluator veval(vertices[0]);
+ PLVineyard::VertexComparison vcmp(veval);
+ PLVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ PLVineyard v(boost::counting_iterator<Vertex>(0),
+ boost::counting_iterator<Vertex>(vertices[0].size()),
+ simplices, veval);
+ std::cout << "Pairing computed" << std::endl;
+
+ // Compute vineyard
+ for (size_t i = 1; i < vertices.size(); ++i)
+ {
+ veval = VertexEvaluator(vertices[i]);
+ v.compute_vineyard(veval, explicit_events);
+ std::cout << "Processed frame: " << i << std::endl;
+ }
+ std::cout << "Vineyard computed" << std::endl;
+
+ if (save_vines)
+ v.vineyard().save_vines(output_prefix, skip_infinite_vines);
+ else
+ v.vineyard().save_edges(output_prefix, skip_infinite_vines);
+}
+
+
+void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices)
+{
+ std::ifstream in(complex_fn.c_str());
+ std::string line;
+ while (std::getline(in, line))
+ {
+ std::istringstream strin(line);
+ simplices.push_back(Smplx(std::istream_iterator<Vertex>(strin), std::istream_iterator<Vertex>()));
+ }
+ std::cout << "Simplices read:" << std::endl;
+ std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Smplx>(std::cout, "\n"));
+}
+
+void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices)
+{
+ std::ifstream in(vertex_fn.c_str());
+ std::string line;
+ while (std::getline(in, line))
+ {
+ std::istringstream strin(line);
+ vertices.push_back(VertexVector(std::istream_iterator<VertexValue>(strin), std::istream_iterator<VertexValue>()));
+ }
+ std::cout << "Vertex values read:" << std::endl;
+ for (size_t i = 0; i < vertices.size(); ++i)
+ {
+ std::copy(vertices[i].begin(), vertices[i].end(), std::ostream_iterator<VertexValue>(std::cout, " "));
+ std::cout << std::endl;
+ }
+}
+
+void program_options(int argc, char* argv[], std::string& complex_fn,
+ std::string& values_fn,
+ std::string& output_prefix,
+ bool& skip_infinite_vines,
+ bool& save_vines,
+ bool& explicit_events)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("complex-file", po::value<std::string>(&complex_fn), "file listing the simplices of the complex")
+ ("values-file", po::value<std::string>(&values_fn), "file listing the values at the vertices")
+ ("output-prefix", po::value<std::string>(&output_prefix), "output prefix");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("skip-infinite,s", po::bool_switch(&skip_infinite_vines), "skip infinite vines in the output")
+ ("explicit-events,e", po::bool_switch(&explicit_events), "process kinetic sort events one by one")
+ ("save-vines,v", po::bool_switch(&save_vines), "save vines instead of edges")
+ ("help,h", "produce help message");
+#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("complex-file", 1);
+ pos.add("values-file", 1);
+ pos.add("output-prefix", 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("complex-file") || !vm.count("values-file") || !vm.count("output-prefix"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] complex-file values-file output-prefix" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/test-grid2D-vineyard.cpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,54 @@
+#include <utilities/log.h>
+
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+#include "grid2D.h"
+#include <topology/lsvineyard.h>
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology"));
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/persistence/transpositions"));
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
+ stdoutLog.subscribeTo(RLOG_CHANNEL("lsvineyard"));
+#endif
+
+ Grid2D g0(2, 2), g1(2, 2);
+ g0(0,0) = 1; g0(0,1) = 2; g0(1,0) = 3; g0(1,1) = 0;
+ g1(0,0) = 4; g1(0,1) = 2; g1(1,0) = 3; g1(1,1) = 5;
+
+ // Generate the complex, initialize the vineyard (which also computes the pairing)
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
+
+ Grid2DVineyard::LSFiltration simplices;
+ g0.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(g0);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ std::cout << "Complex generated, size: " << simplices.size() << std::endl;
+ std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Grid2D::Smplx>(std::cout, "\n"));
+
+ Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
+ std::cout << "Pairing computed" << std::endl;
+
+ // Simplex order before
+ std::cout << "Simplex order:" << std::endl;
+ for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
+ std::cout << " " << v.filtration().simplex(cur) << std::endl;
+
+ // Compute vineyard
+ v.compute_vineyard(g1, true);
+ std::cout << "Vineyard computed" << std::endl;
+
+ // Simplex order after
+ std::cout << "Simplex order:" << std::endl;
+ for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
+ std::cout << " " << v.filtration().simplex(cur) << std::endl;
+
+ v.vineyard().save_edges("test-vineyard");
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/pl-functions/test-grid2D.cpp Sun Dec 20 10:43:00 2009 -0800
@@ -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;
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/lsvineyard.h Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,245 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2009
+ */
+
+#ifndef __LSVINEYARD_H__
+#define __LSVINEYARD_H__
+
+#include <iostream>
+
+#include "topology/simplex.h"
+#include "topology/dynamic-persistence.h"
+#include "topology/lowerstarfiltration.h"
+#include "topology/vineyard.h"
+
+#include <utilities/indirect.h>
+
+#include <CGAL/Kinetic/Inexact_simulation_traits.h>
+#include <CGAL/Kinetic/Sort.h>
+#include <CGAL/Kinetic/Sort_visitor_base.h>
+
+#include <list>
+
+
+template<class Vertex_, class VertexEvaluator_, class Simplex_ = Simplex<Vertex_>, class Filtration_ = Filtration<Simplex_> >
+class LSVineyard
+{
+ public:
+ typedef LSVineyard Self;
+
+ typedef Vertex_ Vertex;
+ typedef VertexEvaluator_ VertexEvaluator;
+ typedef typename VertexEvaluator::result_type VertexValue;
+
+ typedef Simplex_ Simplex;
+ typedef Filtration_ LSFiltration;
+ typedef typename LSFiltration::Index LSFIndex;
+
+ class KineticVertexType;
+ class KineticVertexComparison;
+ typedef std::list<KineticVertexType> VertexContainer;
+ typedef typename VertexContainer::iterator VertexIndex;
+
+ struct AttachmentData: public VineData
+ {
+ void set_attachment(VertexIndex v) { attachment = v; }
+ VertexIndex attachment;
+ };
+ typedef DynamicPersistenceTrails<AttachmentData> Persistence;
+ typedef typename Persistence::OrderIndex Index;
+ typedef typename Persistence::iterator iterator;
+
+ typedef typename Persistence::template SimplexMap<LSFiltration>
+ PFMap;
+
+ class Evaluator;
+ class StaticEvaluator;
+ class KineticEvaluator;
+ class DimensionFromIterator;
+
+ typedef ThroughEvaluatorComparison<VertexEvaluator> VertexComparison;
+ typedef MaxVertexComparison<Simplex, VertexComparison> SimplexComparison;
+
+ class TranspositionVisitor;
+ friend class TranspositionVisitor;
+
+ typedef Vineyard<Index, iterator, Evaluator> Vnrd;
+
+ class SortVisitor;
+ typedef CGAL::Kinetic::Inexact_simulation_traits 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, VertexIndex> KeyVertexMap;
+
+ public:
+ template<class VertexIterator>
+ LSVineyard(VertexIterator begin, VertexIterator end,
+ LSFiltration& filtration,
+ const VertexEvaluator& veval = VertexEvaluator());
+ ~LSVineyard();
+
+ void compute_vineyard(const VertexEvaluator& veval, bool explicit_events = false);
+ bool transpose_vertices(VertexIndex vi);
+
+ const LSFiltration& filtration() const { return filtration_; }
+ const Vnrd& vineyard() const { return vineyard_; }
+ const Persistence& persistence() const { return persistence_; }
+ const VertexComparison& vertex_comparison() const { return vcmp_; }
+ const VertexEvaluator& vertex_evaluator() const { return veval_; }
+ const SimplexComparison& simplex_comparison() const { return scmp_; }
+
+ VertexValue vertex_value(const Vertex& v) const { return veval_(v); }
+ VertexValue simplex_value(const Simplex& s) const { return vertex_value(*std::max_element(s.vertices().begin(), s.vertices().end(), vcmp_)); }
+ const Simplex& pfmap(iterator i) const { return pfmap_[i]; }
+ const Simplex& pfmap(Index i) const { return pfmap_[i]; }
+
+ Index index(iterator i) const { return persistence_.index(i); }
+
+ public:
+ // For Kinetic Sort
+ void swap(Key a, Key b);
+
+ private:
+ void change_evaluator(Evaluator* eval);
+ void set_attachment(iterator i, VertexIndex vi) { persistence_.modifier()(i, boost::bind(&AttachmentData::set_attachment, bl::_1, vi)); }
+ void transpose_filtration(iterator i) { filtration_.transpose(filtration_.begin() + (i - persistence_.begin())); }
+
+ private:
+ VertexContainer vertices_;
+
+ VertexEvaluator veval_;
+ VertexComparison vcmp_;
+ SimplexComparison scmp_;
+
+ LSFiltration& filtration_;
+ Persistence persistence_;
+ PFMap pfmap_;
+
+ Vnrd vineyard_;
+ Evaluator* evaluator_;
+ unsigned time_count_;
+
+ KeyVertexMap kinetic_map_;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ LSVineyard() {}
+
+ 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(LSVineyard)
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticVertexType
+{
+ public:
+ KineticVertexType(const Vertex& v):
+ vertex_(v) {}
+
+ Key kinetic_key() const { return key_; }
+ void set_kinetic_key(Key k) { key_ = k; }
+
+ Vertex vertex() const { return vertex_; }
+ void set_vertex(Vertex v) { vertex_ = v; }
+
+ iterator simplex_index() const { return simplex_index_; }
+ void set_simplex_index(iterator i) { simplex_index_ = i; }
+
+ private:
+ Key key_;
+ Vertex vertex_;
+ iterator simplex_index_;
+};
+
+template<class V, class VE, class S, class C>
+std::ostream&
+operator<<(std::ostream& out, const typename LSVineyard<V,VE,S,C>::VertexIndex& vi)
+{ return out << vi->vertex(); }
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticVertexComparison: public std::binary_function<const KineticVertexType&, const KineticVertexType&, bool>
+{
+ public:
+ KineticVertexComparison(const VertexComparison& vcmp):
+ vcmp_(vcmp) {}
+
+ bool operator()(const KineticVertexType& v1, const KineticVertexType& v2) const
+ { return vcmp_(v1.vertex(), v2.vertex()); }
+
+ private:
+ VertexComparison vcmp_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::TranspositionVisitor: public Persistence::TranspositionVisitor
+{
+ public:
+ typedef typename Persistence::TranspositionVisitor Parent;
+ typedef typename LSVineyard<V,VE,S,C>::iterator iterator;
+ typedef typename LSVineyard<V,VE,S,C>::Index Index;
+
+ TranspositionVisitor(LSVineyard& v): lsvineyard_(v) {}
+
+ void transpose(iterator i) { lsvineyard_.transpose_filtration(i); }
+ void switched(iterator i, SwitchType type) { lsvineyard_.vineyard_.switched(index(i), index(boost::next(i))); }
+
+ private:
+ Index index(iterator i) { return lsvineyard_.index(i); }
+
+ LSVineyard& lsvineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::Evaluator: public std::unary_function<Index, RealType>
+{
+ public:
+ virtual RealType time() const =0;
+ virtual RealType operator()(Index i) const =0;
+ virtual Dimension dimension(Index i) const =0;
+ virtual RealType operator()(iterator i) const { return operator()(&*i); }
+ virtual Dimension dimension(iterator i) const { return dimension(&*i); }
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
+{
+ public:
+ SortVisitor(LSVineyard& v):
+ vineyard_(v) {}
+
+ template<class Vertex_handle>
+ void pre_swap(Vertex_handle a, Vertex_handle b) const { vineyard_.swap(*a,*b); }
+
+ private:
+ LSVineyard& vineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::DimensionFromIterator: std::unary_function<iterator, Dimension>
+{
+ public:
+ DimensionFromIterator(const PFMap& pfmap): pfmap_(pfmap) {}
+
+ Dimension operator()(iterator i) const { return pfmap_[i].dimension(); }
+
+ private:
+ const PFMap& pfmap_;
+};
+
+#include "lsvineyard.hpp"
+
+#endif // __LSVINEYARD_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/lsvineyard.hpp Sun Dec 20 10:43:00 2009 -0800
@@ -0,0 +1,254 @@
+#include <utilities/log.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlLSVineyard = DEF_CHANNEL("lsvineyard/info", rlog::Log_Debug);
+static rlog::RLogChannel* rlLSVineyardDebug = DEF_CHANNEL("lsvineyard/debug", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cVertexTransposition = GetCounter("lsfiltration/transposition"); // counts number of vertex transpositions
+static Counter* cAttachment = GetCounter("lsfiltration/attachment"); // counts the number of attachment changes
+#endif
+
+
+template<class V, class VE, class S, class F>
+template<class VertexIterator>
+LSVineyard<V,VE,S,F>::
+LSVineyard(VertexIterator begin, VertexIterator end,
+ LSFiltration& filtration,
+ const VertexEvaluator& veval):
+ vertices_(begin, end), filtration_(filtration),
+ persistence_(filtration_),
+ veval_(veval), vcmp_(veval_), scmp_(vcmp_),
+ pfmap_(persistence_.make_simplex_map(filtration_)),
+ time_count_(0)
+{
+ // TODO: LSVineyard really should sort the filtration_ itself
+ // filtration_.sort(scmp_);
+ // persistence_.init(filtration_);
+
+ rLog(rlLSVineyardDebug, "Initializing LSVineyard");
+ persistence_.pair_simplices();
+ rLog(rlLSVineyardDebug, "Simplices paired");
+
+ vertices_.sort(KineticVertexComparison(vcmp_)); // sort vertices w.r.t. vcmp_
+#if LOGGING
+ rLog(rlLSVineyardDebug, "Vertex order:");
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ rLog(rlLSVineyardDebug, " %d", cur->vertex());
+#endif
+
+ // Initialize simplex_index() and attachment
+ VertexIndex vi = boost::prior(vertices_.begin());
+ for (iterator i = persistence_.begin(); i != persistence_.end(); ++i)
+ {
+ const Simplex& s = pfmap_[i];
+ if (s.dimension() == 0)
+ {
+ ++vi;
+ AssertMsg(s.vertices().front() == vi->vertex(), "In constructor, simplices and vertices must match.");
+ vi->set_simplex_index(i);
+ }
+ set_attachment(i, vi);
+ rLog(rlLSVineyardDebug, "%s attached to %d", tostring(pfmap(i)).c_str(), vi->vertex());
+ }
+
+ evaluator_ = new StaticEvaluator(*this, time_count_);
+ vineyard_.set_evaluator(evaluator_);
+ vineyard_.start_vines(persistence_.begin(), persistence_.end());
+}
+
+template<class V, class VE, class S, class F>
+LSVineyard<V,VE,S,F>::
+~LSVineyard()
+{
+ delete evaluator_;
+}
+
+template<class V, class VE, class S, class F_>
+void
+LSVineyard<V,VE,S,F_>::
+compute_vineyard(const VertexEvaluator& veval, bool explicit_events)
+{
+ typedef Traits::Kinetic_kernel::Point_1 Point;
+ typedef Traits::Kinetic_kernel::Function_kernel::Construct_function CF;
+ typedef Traits::Kinetic_kernel::Motion_function F;
+
+ Traits tr(0,1);
+ Simulator::Handle sp = tr.simulator_handle();
+ ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
+ Sort sort(tr, SortVisitor(*this));
+
+ // Setup the (linear) trajectories
+ rLog(rlLSVineyard, "Setting up trajectories");
+ CF cf;
+ kinetic_map_.clear();
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ {
+ VertexValue val0 = veval_(cur->vertex());
+ VertexValue val1 = veval(cur->vertex());
+ rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
+ F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
+ Point p(x);
+ cur->set_kinetic_key(apt->insert(p));
+ kinetic_map_[cur->kinetic_key()] = cur;
+ rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
+ }
+
+ // Process all the events (compute the vineyard in the process)
+ change_evaluator(new KineticEvaluator(*this, sp, apt, time_count_));
+ if (explicit_events)
+ {
+ while (sp->next_event_time() < 1)
+ {
+ rLog(rlLSVineyardDebug, "Next event time: %f", sp->next_event_time());
+ sp->set_current_event_number(sp->current_event_number() + 1);
+ rLog(rlLSVineyardDebug, "Processed event");
+ }
+ } else
+ sp->set_current_time(1.0);
+ rLog(rlLSVineyard, "Processed %d events", sp->current_event_number());
+
+ veval_ = veval;
+ change_evaluator(new StaticEvaluator(*this, ++time_count_));
+ vineyard_.record_diagram(persistence_.begin(), persistence_.end());
+}
+
+template<class V, class VE, class S, class F>
+void
+LSVineyard<V,VE,S,F>::
+swap(Key a, Key b)
+{
+ rLog(rlLSVineyardDebug, "Entered swap");
+ VertexIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
+ AssertMsg(vcmp_(ao->vertex(), bo->vertex()), "In swap(a,b), a must precede b");
+ transpose_vertices(ao);
+ // AssertMsg(vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), b must precede a after the transposition");
+}
+
+template<class V, class VE, class S, class F>
+void
+LSVineyard<V,VE,S,F>::
+change_evaluator(Evaluator* eval)
+{
+ AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
+
+ delete evaluator_;
+ evaluator_ = eval;
+ vineyard_.set_evaluator(evaluator_);
+}
+
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
+template<class V, class VE, class S, class F>
+bool
+LSVineyard<V,VE,S,F>::
+transpose_vertices(VertexIndex vi)
+{
+ Count(cVertexTransposition);
+ rLog(rlLSVineyard, "Transposing vertices (%d, %d)", vi->vertex(), boost::next(vi)->vertex());
+
+ DimensionFromIterator dim(pfmap_);
+ TranspositionVisitor visitor(*this);
+
+ iterator i = vi->simplex_index();
+ iterator i_prev = boost::prior(i);
+ iterator i_next = boost::next(vi)->simplex_index();
+ iterator i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
+ iterator j = boost::next(i_next);
+
+ VertexIndex vi_next = boost::next(vi);
+ const Vertex& v = vi->vertex();
+
+ bool result = false; // has a switch in pairing occurred
+
+ // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
+ rLog(rlLSVineyardDebug, "Starting to move the vertex");
+ while (i_next_prev != i_prev)
+ {
+ rLog(rlLSVineyardDebug, " Transposing %s %s", tostring(pfmap(i_next_prev)).c_str(),
+ tostring(pfmap(boost::next(i_next_prev))).c_str());
+ result |= persistence_.transpose(i_next_prev, dim, visitor);
+ i_next_prev = boost::prior(i_next);
+ }
+ rLog(rlLSVineyardDebug, "Done moving the vertex");
+
+ // Second, move the simplices attached to it
+ rLog(rlLSVineyardDebug, "Moving attached simplices");
+ while (j != persistence_.end() && j->attachment == vi_next)
+ {
+ rLog(rlLSVineyardDebug, " Considering %s", tostring(pfmap(j)).c_str());
+ if (pfmap(j).contains(v)) // j becomes attached to v and does not move
+ {
+ Count(cAttachment);
+ rLog(rlLSVineyardDebug, " Attachment changed for ", tostring(pfmap(j)).c_str());
+ set_attachment(j, vi);
+ ++j;
+ continue;
+ }
+
+ iterator j_prev = j; ++j;
+ while ((--j_prev)->attachment != vi_next) // i.e., until we have reached vi_next (and the simplices that follow it) again
+ {
+ rLog(rlLSVineyardDebug, " Moving: %s, %s",
+ tostring(pfmap(j_prev)).c_str(),
+ tostring(pfmap(boost::next(j_prev))).c_str());
+ AssertMsg(j_prev->attachment == vi, "Simplex preceding the one being moved must be attached to v");
+ result |= persistence_.transpose(j_prev, dim, visitor);
+ --j_prev;
+ }
+ }
+ rLog(rlLSVineyard, "Done moving attached simplices");
+ vertices_.splice(vi, vertices_, vi_next); // swap vi and vi_next
+
+ return result;
+}
+
+
+/* Evaluators */
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::StaticEvaluator: public Evaluator
+{
+ public:
+ StaticEvaluator(const LSVineyard& v, RealType time):
+ time_(time), vineyard_(v) {}
+
+ virtual RealType time() const { return time_; }
+ virtual RealType operator()(Index i) const { return vineyard_.simplex_value(vineyard_.pfmap(i)); }
+ virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
+
+ private:
+ RealType time_;
+ const LSVineyard& vineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticEvaluator: public Evaluator
+{
+ public:
+ KineticEvaluator(const LSVineyard& v, Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset):
+ vineyard_(v), sp_(sp), apt_(apt), time_offset_(time_offset) {}
+
+ virtual RealType time() const { return time_offset_ + CGAL::to_double(get_time()); }
+ virtual RealType operator()(Index i) const
+ {
+ rLog(rlLSVineyard, "%s (attached to %d): %s(%f) = %f", tostring(vineyard_.pfmap(i)).c_str(),
+ i->attachment->vertex(),
+ tostring(apt_->at(i->attachment->kinetic_key()).x()).c_str(),
+ get_time(),
+ CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())));
+ return CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time()));
+ }
+ virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
+
+ private:
+ Simulator::Time get_time() const { return sp_->current_time(); }
+
+ const LSVineyard& vineyard_;
+ Simulator::Handle sp_;
+ ActivePointsTable::Handle apt_;
+ RealType time_offset_;
+};