Add diagram_point::get_persistence(order, func)

[libstick.git] / include / libstick / persistence.h

diff --git a/include/libstick/persistence.h b/include/libstick/persistence.h
index b9fc080ee589c951f9a16a0e354dcc774888028e..72b07b813944d758fe4927484d4fb03f96651926 100644 (file)
--- a/include/libstick/persistence.h
+++ b/include/libstick/persistence.h
@@ -6,7 +6,7 @@
 #include <vector>
 #include <cassert>
 
-#include "simplicialcomplex.h"
+#include "simplicialfunction.h"
 #include "booleanmatrix.h"
 
 
@@ -15,11 +15,11 @@ namespace libstick {
 
 /** Persistence computers various persistent homology information on a
  * simplex_order of a simplicial_complex. */
-template<int MAXDIM, class IT, class VT>
+template<int MAXDIM, class IT>
 class persistence {
 
     public:
-        typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
+        typedef simplicial_complex<MAXDIM, IT> scomplex;
         typedef typename scomplex::simplex_order simplex_order;
         typedef boolean_colmatrix<IT> boundary_matrix;
         typedef boolean_colmatrix<IT> transformation_matrix;
@@ -35,6 +35,34 @@ class persistence {
              * zero, the cycle never dies. Otherwise, 'death' is always larger
              * than 'birth'. */
             IT death;
+
+
+            /** Get the value of the birth index w.r.t. to a simplicial order o
+             * and a simplicial function f. Requires that o and f are defined
+             * on the same simplicial complex. */
+            template<class VT>
+            VT get_birth_value(const simplex_order &o, const simplicial_function<MAXDIM, IT, VT> &f) const {
+                assert(&f.get_complex() == &o.get_complex());
+                return f.get_value(o.order_to_complex(birth));
+            }
+
+            /** Get the value of the death index w.r.t. to a simplicial order o
+             * and a simplicial function f. Requires that o and f are defined
+             * on the same simplicial complex. */
+            template<class VT>
+            VT get_death_value(const simplex_order &o, const simplicial_function<MAXDIM, IT, VT> &f) const {
+                assert(&f.get_complex() == &o.get_complex());
+                return f.get_value(o.order_to_complex(death));
+            }
+
+            /** Get the persistence of this point w.r.t. to a simplicial order
+             * o and a simplicial function f. Requires that o and f are defined
+             * on the same simplicial complex. */
+            template<class VT>
+            VT get_persistence(const simplex_order &o, const simplicial_function<MAXDIM, IT, VT> &f) const {
+                assert(&f.get_complex() == &o.get_complex());
+                return get_death_value<VT>(o, f) - get_birth_value<VT>(o, f);
+            }
         };
 
         /** Save the births and deaths of simplices of a fixed dimension, say p. */
@@ -172,8 +200,11 @@ class persistence {
             return tm;
         }
 
-        /** Print the sequence of births and deaths of cycles (homology classes). */
-        void interpret_reduction(std::ostream &os) const {
+        /** Print the sequence of births and deaths of cycles (homology
+         * classes). If a non-NULL simplicial function is given, also the
+         * values of the simplices are printed.  */
+        template<class VT>
+        void interpret_reduction(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
             assert(done_matrices);
 
             for (unsigned c=0; c < bm.width(); ++c) {
@@ -199,21 +230,28 @@ class persistence {
                         os << "simplex";
                         break;
                 }
+                os << std::endl;
 
-                os << " \t\tvalue: " << order.get_simplex(c).value << std::endl;
+                if (f != NULL)
+                    os << "  value: " <<  f->get_value(order.order_to_complex(c)) << std::endl;
                 os << "  boundary: " << bm.get_column(c) << std::endl;
 
                 typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
                 if (col.isempty()) {
                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
                 } else {
-                    os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
-                    os << "  boundary of: " << tm.get_column(c) << std::endl;
+                    os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c);
+                    os << ", boundary of: " << tm.get_column(c) << std::endl;
                 }
                 os << std::endl;
             }
         }
 
+        /** Print the sequence of births and deaths of cycles (homology classes). */
+        void interpret_reduction(std::ostream &os) const {
+            interpret_reduction<float>(os, NULL);
+        }
+
         /** Get the 'dim'-dimensional peristence diagram. */
         const diagram& get_persistence_diagram(int dim) const{
             assert(-1 <= dim);
@@ -273,8 +311,11 @@ class persistence {
 #endif
         }
 
-        /** Print birth and death information in the persistence diagrams. */
-        void interpret_persistence(std::ostream &os) const {
+        /** Print birth and death information in the persistence diagrams. If a
+         * non-NULL simplicial function is given, also the persistence w.r.t.
+         * this function is printed. */
+        template<class VT>
+        void interpret_persistence(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
             assert(done_diagrams);
 
             for (int d=-1; d <= MAXDIM; ++d) {
@@ -287,17 +328,25 @@ class persistence {
 
                     os << "  ";
                     if (death > 0) {
-                        const scomplex& c = get_order().get_complex();
-                        const VT pers = c.simplices[death].value - c.simplices[birth].value;
-                        os << birth << "\033[1;31m -> " << death << "\033[0;m" << " \tpers: " << pers;
+                        os << birth << "\033[1;31m -> " << death << "\033[0;m";
+
+                        if (f != NULL)
+                            os << " \tpers: " << dia.births[i].get_persistence(get_order(), *f);
                     } else {
-                        os << "\033[1;32m" << birth << "\033[0;m" << " \t\t[essential]";
+                        os << "\033[1;32m" << birth << "\033[0;m";
+                        if (f != NULL)
+                            os << " \t\t[essential]";
                     }
                     os << std::endl;
                 }
             }
         }
 
+        /** Print birth and death information in the persistence diagrams. */
+        void interpret_persistence(std::ostream &os) const {
+            interpret_persistence<float>(os, NULL);
+        }
+
     private:
 
         diagram& _get_persistence_diagram(int dim) {