--- a/doc/examples/index.rst	Thu Jul 09 07:47:22 2009 -0700
+++ b/doc/examples/index.rst	Fri Jul 24 11:08:48 2009 -0700
@@ -42,3 +42,6 @@
     :maxdepth: 1
 
     cohomology
+
+A simple example of computing persistence of a lower-star filtration is in
+:sfile:`examples/lsfiltration.py`.    

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/consistency/make-zigzag-subsamples.py	Fri Jul 24 11:08:48 2009 -0700
@@ -0,0 +1,47 @@
+#!/usr/bin/env python 
+
+# Creates POINTS and SUBSAMPLES files from a list of points file in a format
+# suitable for rips-consistency-zigzag.
+
+
+from    sys     import argv, exit
+
+
+def create_subsamples(points_fn, subsamples_fn, points_list):
+    points = []
+    count = []
+    for pfn in points_list:
+        count.append(0)
+        with open(pfn) as f:
+            for line in f:
+                if line.startswith('#'): continue
+                points.append(line)
+                count[-1] += 1
+
+    with open(points_fn, 'w') as f:
+        for line in points:
+            f.write(line)
+
+    cur = 0
+    counts = []
+    for c in count:
+        counts.append(' '.join(map(str, xrange(cur, cur+c))) + '\n')
+        cur += c
+        counts.append(' '.join(map(str, xrange(cur-c, cur+c))) + '\n')
+    
+    with open(subsamples_fn, 'w') as f:
+        f.writelines(counts[:-1])
+
+
+if __name__ == '__main__':
+    if len(argv) < 4:
+        print "Usage: %s POINTS SUBSAMPLES POINTS1 [POINTS2 [POINTS3 [...]]]" % argv[0]
+        print
+        print "Creates a file POINTS with the union of POINTS* and SUBSAMPLES which lists"
+        print "the indices of the points and their pairwise unions, one per line"
+        exit()
+
+    points_fn = argv[1]
+    subsamples_fn = argv[2]
+    points_list = argv[3:]
+    create_subsamples(points_fn, subsamples_fn, points_list)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/lsfiltration.py	Fri Jul 24 11:08:48 2009 -0700
@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+
+from    dionysus        import Simplex, Filtration, StaticPersistence, vertex_cmp
+from    sys             import argv, exit
+
+
+def max_vertex(s, vertices):
+    return max((vertices[v] for v in s.vertices))
+
+def max_vertex_cmp(s1, s2, vertices):
+    m1 = max_vertex(s1, vertices)
+    m2 = max_vertex(s2, vertices)
+    return cmp(m1, m2) or cmp(s1.dimension(), s2.dimension())
+
+def lsf(values_filename, simplices_filename):
+    # Read vertices
+    vertices = []
+    with open(values_filename) as f:
+        for line in f:
+            if line.startswith('#'): continue
+            vertices.append(float(line.split()[1]))
+
+    # Read simplices
+    simplices = []            
+    with open(simplices_filename) as f:
+        for line in f:
+            if line.startswith('#'): continue
+            simplices.append(map(int, line.split()))
+
+    # Setup the complex
+    complex = [Simplex(s) for s in simplices]
+    complex.sort(vertex_cmp)
+
+    # Compute persistence
+    f = Filtration(complex, lambda x,y: max_vertex_cmp(x,y,vertices))
+    p = StaticPersistence(f)
+    p.pair_simplices()
+    
+    # Output the persistence diagram
+    for i in p:
+        if not i.sign(): continue
+
+        b = complex[f[p(i)]]
+        d = complex[f[p(i.pair)]]
+
+        if i == i.pair:
+            print b.dimension(), max_vertex(b, vertices), "inf"
+            continue
+
+        print b.dimension(), max_vertex(b, vertices), max_vertex(d, vertices)
+
+
+if __name__ == '__main__':
+    if len(argv) < 3:
+        print "Usage: %s VERTICES SIMPLICES" % argv[0]
+        print 
+        print "Computes persistence of the lower star filtration of the simplicial "
+        print "complex explicitly listed out in SIMPLICES with vertex values given in VERTICES."
+        exit()
+
+    values = argv[1]
+    simplices = argv[2]
+
+    lsf(values, simplices)

--- a/examples/rips/rips-pairwise.cpp	Thu Jul 09 07:47:22 2009 -0700
+++ b/examples/rips/rips-pairwise.cpp	Fri Jul 24 11:08:48 2009 -0700
@@ -48,7 +48,7 @@
     // Generate 2-skeleton of the Rips complex for epsilon = 50
     rips.generate(skeleton, max_distance, make_push_back_functor(complex));
     std::sort(complex.begin(), complex.end(), Smplx::VertexComparison());       // unnecessary
-    std::cout << "Generated complex of size: " << complex.size() << std::endl;
+    std::cout << "# Generated complex of size: " << complex.size() << std::endl;
 
     // Generate filtration with respect to distance and compute its persistence
     Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
@@ -70,30 +70,34 @@
             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;
+            // if (b.dimension() != 1) continue;
+            // std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
+            if (b.dimension() >= skeleton) continue;
+            std::cout << b.dimension() << " " << 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;
+            // if (b.dimension() != 1) continue;
             
-            std::cout << "Unpaired birth: " << size(b) << std::endl;
-            cycle = cur->chain;
+            // std::cout << "Unpaired birth: " << size(b) << std::endl;
+            // cycle = cur->chain;
+            if (b.dimension() >= skeleton) continue;
+            std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
         }
 
         // Iterate over the cycle
-        for (Persistence::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;
-        }
+        // for (Persistence::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");
+    persistence_timer.check("# Persistence timer");
 }
 
 void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)

--- a/include/geometry/l2distance.h	Thu Jul 09 07:47:22 2009 -0700
+++ b/include/geometry/l2distance.h	Fri Jul 24 11:08:48 2009 -0700
@@ -2,6 +2,7 @@
 #define __L2_DISTANCE_H__
 
 #include <utilities/types.h>
+#include <utilities/log.h>
 
 #include <vector>
 #include <fstream>