X-Git-Url: https://git.sthu.org/?a=blobdiff_plain;f=include%2Flibstick-0.1%2Fpersistence.h;h=b9fc080ee589c951f9a16a0e354dcc774888028e;hb=909909ce5de57f1454686cf6251d0bd9652a96d9;hp=6df70abc7c8345b0f11a86d168fdf3bcd9f4f0f2;hpb=799cd8e53b06b4f9ec7d4dde13ab600717982099;p=libstick.git

diff --git a/include/libstick-0.1/persistence.h b/include/libstick-0.1/persistence.h
index 6df70ab..b9fc080 100644
--- a/include/libstick-0.1/persistence.h
+++ b/include/libstick-0.1/persistence.h
@@ -94,6 +94,11 @@ class persistence {
             done_diagrams = false;
         }
 
+        /** Get simplicial order of this persistence object. */
+        const simplex_order& get_order() const {
+            return order;
+        }
+
         /** Get boundary matrix 'bm' of 'order', compute reduces boundary
          * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
          * tm. */
@@ -112,7 +117,7 @@ class persistence {
                 //if (c % 100 == 0)
                 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
                 // Reduce as long as we need to reduce
-                while (rbm.get_column(c).size() > 0) {
+                while (!rbm.get_column(c).isempty()) {
                     // (r, c) is the lowest one of col
                     const IT r = rbm.get_column(c).back();
                     // This column contains its lowest one on the same row
@@ -132,7 +137,7 @@ class persistence {
                 }
 
                 // A lowest one remained, recall it
-                if (rbm.get_column(c).size() > 0) {
+                if (!rbm.get_column(c).isempty()) {
                     const IT r = rbm.get_column(c).back();
                     assert(lowestones[r] == 0);
                     assert(r < c);
@@ -195,11 +200,11 @@ class persistence {
                         break;
                 }
 
-                os << std::endl;
+                os << " \t\tvalue: " << order.get_simplex(c).value << std::endl;
                 os << "  boundary: " << bm.get_column(c) << std::endl;
 
                 typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
-                if (col.size() == 0) {
+                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;
@@ -235,7 +240,7 @@ class persistence {
 
             for (unsigned c=0; c < lowestones.size(); ++c) {
                 // A cycle was born
-                if (rbm.get_column(c).size() == 0) {
+                if (rbm.get_column(c).isempty()) {
                     const int dim = order.get_simplex(c).dim;
 
                     // Create a diagram point
@@ -260,7 +265,7 @@ class persistence {
             // or not be, tertium non datur.
             for (unsigned c=0; c < lowestones.size(); ++c) {
                 // A class died with c
-                if (rbm.get_column(c).size() != 0)
+                if (!rbm.get_column(c).isempty())
                     assert(c == 0 || deaths.count(c) > 0);
                 else
                     assert(births.count(c) > 0);
@@ -275,13 +280,20 @@ class persistence {
             for (int d=-1; d <= MAXDIM; ++d) {
                 os << "Dimension " << d << std::endl;
                 const diagram &dia = get_persistence_diagram(d);
+
                 for (unsigned i=0; i < dia.births.size(); ++i) {
+                    const IT birth = dia.births[i].birth;
+                    const IT death = dia.births[i].death;
+
                     os << "  ";
-                    if (dia.births[i].death > 0)
-                        os << dia.births[i].birth << "\033[1;31m -> " << dia.births[i].death;
-                    else
-                        os << "\033[1;32m" << dia.births[i].birth;
-                    os << "\033[0;m" << std::endl;
+                    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;
+                    } else {
+                        os << "\033[1;32m" << birth << "\033[0;m" << " \t\t[essential]";
+                    }
+                    os << std::endl;
                 }
             }
         }
@@ -319,54 +331,6 @@ class persistence {
 };
 
 
-/** Takes the boundary matrix b and transforms it inplace into a reduced matrix
- * b. The matrix v is required to be a unit matrix having the same dimensions
- * as b. It is maintained in order to leave the product b*inverse(v) as an
- * invariant.  Hence, the resulting b is equal to the product of the boundary
- * matrix times v. */
-template<class IT, class D>
-void reduce_boundary_matrix(boolean_colmatrix_base<IT, D> &b, boolean_colmatrix_base<IT, D> &v) {
-    assert(b.width() == v.width());
-
-    // lowestones[i] gives the column whose lowest 1 is at row i. If it contains
-    // zero, there is no such column.
-    std::vector<IT> lowestones(b.width());
-
-    // Make every column reduced, i.e., it is a zero-column or contains only one 1.
-    for (unsigned c=0; c < b.width(); ++c) {
-        //if (c % 100 == 0)
-            //std::cout << "c = " << c << " (" << (100.0*float(c)/b.width()) << " %)" << std::endl;
-        // Reduce as long as we need to reduce
-        while (b.get_column(c).size() > 0) {
-            // (r, c) is the lowest one of col
-            const IT r = b.get_column(c).back();
-            // This column contains its lowest one on the same row
-            const IT c2 = lowestones[r];
-            assert(c2 < c);
-
-            // There is no valid c2, hence we have completely reduced
-            if (c2 == 0)
-                break;
-            // Reduce col by column c2
-            else {
-                assert(b.get(r, c));
-                b.add_column(c, b.get_column(c2));
-                v.add_column(c, v.get_column(c2));
-                assert(!b.get(r, c));
-            }
-        }
-
-        // A lowest one remained, recall it
-        if (b.get_column(c).size() > 0) {
-            const IT r = b.get_column(c).back();
-            assert(lowestones[r] == 0);
-            assert(r < c);
-            lowestones[r] = c;
-        }
-    }
-}
-
-
 } // namespace libstick
 
 #endif