Add diagram_point::get_persistence(order, func)

diff --git a/include/libstick/persistence.h b/include/libstick/persistence.h
index 481d6480d2dd50821d82845672fdf0dd45873673..72b07b813944d758fe4927484d4fb03f96651926 100644 (file)
--- a/include/libstick/persistence.h
+++ b/include/libstick/persistence.h
@@ -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. */
@@ -302,12 +330,8 @@ class persistence {
                     if (death > 0) {
                         os << birth << "\033[1;31m -> " << death << "\033[0;m";
 
-                        if (f != NULL) {
-                            const IT birth_idx = order.order_to_complex(birth);
-                            const IT death_idx = order.order_to_complex(death);
-                            const VT pers = f->get_value(death_idx) - f->get_value(birth_idx);
-                            os << " \tpers: " << pers;
-                        }
+                        if (f != NULL)
+                            os << " \tpers: " << dia.births[i].get_persistence(get_order(), *f);
                     } else {
                         os << "\033[1;32m" << birth << "\033[0;m";
                         if (f != NULL)

diff --git a/tests/persistence.h b/tests/persistence.h
index 186039c53e2aef8e0ee21c087837bc959175dbf3..b7678bfd76993fb3ed03bb45046f3bbcc3b2eb5e 100644 (file)
--- a/tests/persistence.h
+++ b/tests/persistence.h
@@ -34,6 +34,7 @@ class persistence_TestSuite: public Test::Suite {
             TEST_ADD(persistence_TestSuite::test_betti_numbers);
             TEST_ADD(persistence_TestSuite::test_lowestones);
             TEST_ADD(persistence_TestSuite::test_diagram);
+            TEST_ADD(persistence_TestSuite::test_diagrampoint_value);
         }
 
     protected:
@@ -221,6 +222,28 @@ class persistence_TestSuite: public Test::Suite {
 #endif
         }
 
+        void test_diagrampoint_value() {
+            pers torusp(otorus);
+            torusp.compute_diagrams();
+
+            // We make a simplicial function that assigns to each simplex twice
+            // the index as value.
+            sfunction ftorus(torus);
+            ftorus.make_complete();
+            for (unsigned i=1; i < torus.size(); ++i)
+                ftorus.set_value(i, 2*i);
+
+            for (unsigned d=0; d <= 3; ++d) {
+                const pers::diagram &dia = torusp.get_persistence_diagram(d);
+                for (unsigned i=0; i < dia.births.size(); ++i) {
+                    TEST_ASSERT(2 * dia.births[i].birth == dia.births[i].get_birth_value(otorus, ftorus));
+                    TEST_ASSERT(2 * dia.births[i].death == dia.births[i].get_death_value(otorus, ftorus));
+                    TEST_ASSERT(2 * (dia.births[i].death - dia.births[i].birth) ==
+                            dia.births[i].get_persistence(otorus, ftorus));
+                }
+            }
+        }
+
         void test_betti_numbers() {
             pers ringp(oring);
             ringp.compute_diagrams();