Merged in cohomology branch dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Mon, 23 Mar 2009 11:40:29 -0700
branchdev
changeset 117 06a4361bddaa
parent 115 a3410b6ba79c (diff)
parent 116 5095771715ab (current diff)
child 118 c4e25fb4082c
Merged in cohomology branch
examples/rips/rips-pairwise.cpp
--- a/README	Fri Feb 27 10:14:26 2009 -0800
+++ b/README	Mon Mar 23 11:40:29 2009 -0700
@@ -19,10 +19,6 @@
 [rlog]:         http://www.arg0.net/rlog
 [SYNAPS]:       http://www-sop.inria.fr/galaad/synaps/
 
-## Configuration
-  The path to CGAL's Makefile is expected to be set in `$CGAL_MAKEFILE`, the rest
-  is just the usual CMake configuration
-
 ## Building
   To build examples, create a directory build (to keep everything in one place),
   go to that directory and run cmake and make:
--- a/examples/rips/rips-pairwise.cpp	Fri Feb 27 10:14:26 2009 -0800
+++ b/examples/rips/rips-pairwise.cpp	Mon Mar 23 11:40:29 2009 -0700
@@ -1,6 +1,7 @@
 #include <topology/rips.h>
 #include <topology/filtration.h>
 #include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
 #include <topology/persistence-diagram.h>
 
 #include <utilities/containers.h>           // for BackInsertFunctor
@@ -20,7 +21,9 @@
 typedef         Generator::Simplex                                      Smplx;
 typedef         std::vector<Smplx>                                      SimplexVector;
 typedef         Filtration<SimplexVector, unsigned>                     Fltr;
-typedef         StaticPersistence<>                                     Persistence;
+//typedef         StaticPersistence<>                                     Persistence;
+typedef         DynamicPersistenceChains<>                              Persistence;
+typedef         PersistenceDiagram<>                                    PDgm;
 
 void            program_options(int argc, char* argv[], std::string& infilename, Dimension& ambient, Dimension& skeleton, DistanceType& max_distance);
 
@@ -55,15 +58,38 @@
 
     // Output cycles
     for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+    {
+        Persistence::OrderDescriptor::Cycle& cycle = cur->cycle;
+
         if (!cur->sign())        // only negative simplices have non-empty cycles
         {
             Persistence::OrderIndex birth = cur->pair;      // the cycle that cur killed was born when we added birth (another simplex)
 
             const Smplx& b = f.simplex(f.begin() + (birth - p.begin()));        // eventually this will be encapsulated behind an interface
             const Smplx& d = f.simplex(f.begin() + (cur   - p.begin()));
-            if (size(d) - size(b) > 0)
-                std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
+            
+            if (b.dimension() != 1) continue;
+            std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
+        } else if (cur->pair == cur)    // positive could be unpaired
+        {
+            const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+            if (b.dimension() != 1) continue;
+            
+            std::cout << "Unpaired birth: " << size(b) << std::endl;
+            cycle = cur->chain;
         }
+
+        // Iterate over the cycle
+        for (Persistence::OrderDescriptor::Cycle::const_iterator si =  cycle.begin();
+                                                                 si != cycle.end();     ++si)
+        {
+            const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
+            //std::cout << s.dimension() << std::endl;
+            const Smplx::VertexContainer& vertices = s.vertices();          // std::vector<Vertex> where Vertex = Distances::IndexType
+            AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
+            std::cout << vertices[0] << " " << vertices[1] << std::endl;
+        }
+    }
     
     persistence_timer.check("Persistence timer");
 }
--- a/examples/rips/rips.cpp	Fri Feb 27 10:14:26 2009 -0800
+++ b/examples/rips/rips.cpp	Mon Mar 23 11:40:29 2009 -0700
@@ -1,6 +1,7 @@
 #include <topology/rips.h>
 #include <topology/filtration.h>
 #include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
 #include <topology/persistence-diagram.h>
 #include <utilities/containers.h>           // for BackInsertFunctor
 
@@ -26,8 +27,9 @@
 typedef         Generator::Simplex                                      Smplx;
 typedef         std::vector<Smplx>                                      SimplexVector;
 typedef         Filtration<SimplexVector, unsigned>                     Fltr;
-typedef         StaticPersistence<>                                     Persistence;
-typedef         PersistenceDiagram<>            PDgm;
+//typedef         StaticPersistence<>                                     Persistence;
+typedef         DynamicPersistenceChains<>                              Persistence;
+typedef         PersistenceDiagram<>                                    PDgm;
 
 
 int main(int argc, char* argv[])
@@ -51,7 +53,7 @@
     SimplexVector           complex;
     
     // Generate 2-skeleton of the Rips complex for epsilon = 50
-    rips.generate(2, 50, make_push_back_functor(complex));
+    rips.generate(2, 10, make_push_back_functor(complex));
     std::sort(complex.begin(), complex.end(), Smplx::VertexComparison());       // unnecessary
     rInfo("Generated complex of size: %d",  complex.size());
 
@@ -73,4 +75,39 @@
     std::ofstream ofs("rips-diagrams");
     boost::archive::binary_oarchive oa(ofs);
     oa << dgms;
+
+    // Output cycles
+    for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+    {
+        Persistence::OrderDescriptor::Cycle& cycle = cur->cycle;
+
+        if (!cur->sign())        // only negative simplices have non-empty cycles
+        {
+            Persistence::OrderIndex birth = cur->pair;      // the cycle that cur killed was born when we added birth (another simplex)
+
+            const Smplx& b = f.simplex(f.begin() + (birth - p.begin()));        // eventually this will be encapsulated behind an interface
+            const Smplx& d = f.simplex(f.begin() + (cur   - p.begin()));
+            
+            if (b.dimension() != 1) continue;
+            std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
+        } else if (cur->pair == cur)    // positive could be unpaired
+        {
+            const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+            if (b.dimension() != 1) continue;
+            
+            std::cout << "Unpaired birth: " << size(b) << std::endl;
+            cycle = cur->chain;
+        }
+
+        // Iterate over the cycle
+        for (Persistence::OrderDescriptor::Cycle::const_iterator si =  cycle.begin();
+                                                                 si != cycle.end();     ++si)
+        {
+            const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
+            //std::cout << s.dimension() << std::endl;
+            const Smplx::VertexContainer& vertices = s.vertices();          // std::vector<Vertex> where Vertex = Distances::IndexType
+            AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
+            std::cout << vertices[0] << " " << vertices[1] << std::endl;
+        }
+    }
 }
--- a/include/topology/dynamic-persistence.h	Fri Feb 27 10:14:26 2009 -0800
+++ b/include/topology/dynamic-persistence.h	Mon Mar 23 11:40:29 2009 -0700
@@ -4,8 +4,11 @@
 #include "static-persistence.h"
 #include <utilities/types.h>
 
+#include <boost/progress.hpp>
+
 #ifdef COUNTERS
 static Counter*  cTrailLength =             GetCounter("persistence/pair/traillength");     // the size of matrix U in RU decomposition
+static Counter*  cChainLength =             GetCounter("persistence/pair/chainlength");     // the size of matrix V in R=DV decomposition
 #endif // COUNTERS
 
 template<class Data, class OrderDescriptor_, class ConsistencyIndex>
@@ -97,6 +100,7 @@
 
         using                           Parent::begin;
         using                           Parent::end;
+        using                           Parent::size;
 
         // Struct: TranspositionVisitor
         //
@@ -139,6 +143,140 @@
         ConsistencyComparison           ccmp_;
 };
 
+template<class Data, class OrderDescriptor_, class ConsistencyIndex>
+struct ChainData_: public Data
+{
+    typedef ChainData_<Data, OrderDescriptor_, ConsistencyIndex>                Self;
+
+    typedef typename OrderTraits<typename OrderDescriptor_::
+                                          template RebindData<Self>::
+                                          other>::Index                         OrderIndex;
+    typedef typename OrderDescriptor_::
+                     Chains::template rebind<OrderIndex>::other::Chain          Chain;
+    
+    template<class Comparison>
+    struct ConsistencyComparison: public std::binary_function<const OrderIndex&, const OrderIndex&, bool>
+    {
+        ConsistencyComparison(const Comparison& cmp = Comparison()): cmp_(cmp)  {}
+
+        bool operator()(const OrderIndex& a, const OrderIndex& b) const         { return cmp_(a->consistency, b->consistency); }
+
+        Comparison      cmp_;
+    };
+
+    Chain                                                                       chain;
+    ConsistencyIndex                                                            consistency;
+};
+
+/**
+ * Class: DynamicPersistenceChains
+ *
+ * TODO: below comment is incorrect; nothing dynamic about this yet.
+ * Derives from StaticPersistence and allows one to update persistence 
+ * after a transposition of two contiguous simplices in a filtration. 
+ * In addition to reduced cycles, it stores each OrderElement's chains,
+ * i.e. in addition to matrix R, it stores matrix V in vineyard notation.
+ *
+ * Template parameters:
+ *   Data_ -                auxilliary contents to store with each OrderElement
+ *   OrderDescriptor_ -     class describing how the order is stored; it defaults to <VectorOrderDescriptor> 
+ *                          which serves as a prototypical class
+ */
+template<class Data_ =                  Empty<>, 
+         class OrderDescriptor_ =       VectorOrderDescriptor<>,
+         class ConsistencyIndex_ =      size_t,
+         class ConsistencyComparison_ = std::less<ConsistencyIndex_> >
+class DynamicPersistenceChains: 
+    public StaticPersistence<ChainData_<Data_, OrderDescriptor_, ConsistencyIndex_>, 
+                             OrderDescriptor_>
+{
+    public:
+        typedef         Data_                                                           Data;
+        typedef         ChainData_<Data_, OrderDescriptor_, ConsistencyIndex_>          ChainData;
+        typedef         StaticPersistence<ChainData, OrderDescriptor_>                  Parent;
+ 
+        typedef         typename Parent::Traits                                         Traits;
+        typedef         typename Parent::OrderDescriptor                                OrderDescriptor;
+        typedef         typename Parent::OrderComparison                                OrderComparison;
+        typedef         typename Parent::OrderIndex                                     OrderIndex;
+        typedef         ConsistencyIndex_                                               ConsistencyIndex;
+        typedef         ThreeOutcomeCompare<
+                            typename ChainData::
+                            template ConsistencyComparison<ConsistencyComparison_> >    ConsistencyComparison;
+
+        /**
+         * Constructor: DynamicPersistenceChains()
+         * TODO: write a description
+         *
+         * Template parameters:
+         *   Filtration -           filtration of the complex whose persistence we are computing
+         */
+        template<class Filtration>      DynamicPersistenceChains(const Filtration&              f, 
+                                                                 const OrderComparison&         ocmp =  OrderComparison(),
+                                                                 const ConsistencyComparison_&  ccmp =  ConsistencyComparison_());
+        
+        void                            pair_simplices();
+
+        // Function: transpose(i)
+        // Tranpose i and the next element. 
+        // Returns: true iff the pairing switched.
+        // TODO
+        //bool                            transpose(OrderIndex i)                         { return transpose(i, TranspositionVisitor()); }
+        
+        // TODO: the main missing piece to be dynamic
+        //template<class Visitor>
+        //bool                            transpose(OrderIndex i, const Visitor& visitor = Visitor());
+
+        using                           Parent::begin;
+        using                           Parent::end;
+        using                           Parent::size;
+
+        // Struct: TranspositionVisitor
+        //
+        // For example, a VineardVisitor could implement this archetype.
+        struct TranspositionVisitor
+        {
+            // Function: transpose(i)
+            // This function is called before transposition is processed 
+            // (at the very beginning of <transpose(i, visitor)>). It is meant to update any structures 
+            // that may need to be updated, but perhaps it has other uses as well.
+            void                        transpose(OrderIndex i) const                   {}
+
+            // Function: switched(i, type)
+            // This function is called after the transposition if the switch in pairing has occured.
+            // `i` is the index of the preceding simplex after the transposition. 
+            // `type` indicates the <SwitchType>.
+            void                        switched(OrderIndex i, SwitchType type) const   {}
+        };
+
+    protected:
+        using                           Parent::order;
+
+    private:
+        void                            swap(OrderIndex i, OrderIndex j);
+        void                            pairing_switch(OrderIndex i, OrderIndex j);
+
+        struct PairingChainsVisitor: public Parent::PairVisitor 
+        {
+            // TODO: this is specialized for std::vector
+                                        PairingChainsVisitor(OrderIndex bg, ConsistencyComparison ccmp, unsigned size): 
+                                            bg_(bg), ccmp_(ccmp), show_progress(size)   {}
+
+            void                        init(OrderIndex i) const                        { i->consistency = i - bg_; i->chain.append(i, ccmp_); }
+            void                        update(OrderIndex j, OrderIndex i) const        { j->chain.add(i->pair->chain, ccmp_); }
+            void                        finished(OrderIndex i) const                    { CountBy(cChainLength, i->chain.size()); ++show_progress; }
+
+            OrderIndex                  bg_;
+            ConsistencyComparison       ccmp_;
+
+            mutable boost::progress_display     
+                                        show_progress;
+        };
+
+        ConsistencyComparison           ccmp_;
+};
+
+
 #include "dynamic-persistence.hpp"
 
 #endif  // __DYNAMIC_PERSISTENCE_H__
--- a/include/topology/dynamic-persistence.hpp	Fri Feb 27 10:14:26 2009 -0800
+++ b/include/topology/dynamic-persistence.hpp	Mon Mar 23 11:40:29 2009 -0700
@@ -20,6 +20,8 @@
 #endif // COUNTERS
 
 
+/* Trails */
+
 template<class D, class OD, class CI, class CC>
 template<class Filtration>
 DynamicPersistenceTrails<D,OD,CI,CC>::
@@ -256,3 +258,22 @@
         j_pair->pair = i;
     }
 }
+
+
+/* Chains */
+
+template<class D, class OD, class CI, class CC>
+template<class Filtration>
+DynamicPersistenceChains<D,OD,CI,CC>::
+DynamicPersistenceChains(const Filtration& f, const OrderComparison& ocmp, const CC& ccmp):
+    Parent(f, ocmp), ccmp_(ccmp)
+{}
+        
+template<class D, class OD, class CI, class CC>
+void
+DynamicPersistenceChains<D,OD,CI,CC>::
+pair_simplices()
+{ 
+    Parent::pair_simplices(begin(), end(), PairingChainsVisitor(begin(), ccmp_, size()));
+}
+