Rips pairwise cohomology produces output necessary for 1-cocycle parametrization dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu, 09 Jul 2009 02:42:47 -0700
branchdev
changeset 138 96030f8d2f2c
parent 137 069596c71902
child 139 656ec2838fb8
Rips pairwise cohomology produces output necessary for 1-cocycle parametrization
examples/cohomology/cocycle.py
examples/cohomology/lsqr.py
examples/cohomology/output.h
examples/cohomology/rips-pairwise-cohomology.cpp
tools/plot-values/plot.py
tools/plot-values/scatter.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/cocycle.py	Thu Jul 09 02:42:47 2009 -0700
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+
+from    cvxopt          import spmatrix, matrix
+from    cvxopt.blas     import copy
+from    lsqr            import lsqr
+from    sys             import argv, exit
+import  os.path
+
+def smooth(boundary_filename, cocycle_filename, vertexmap_filename):
+    boundary_list = []
+    with open(boundary_filename) as fp:
+        for line in fp.xreadlines():
+            if line.startswith('#'): continue
+            boundary_list.append(map(int, line.split()))
+
+    cocycle_list = []
+    with open(cocycle_filename) as fp:
+        for line in fp.xreadlines():
+            if line.startswith('#'): continue
+            cocycle_list.append(map(int, line.split()))
+
+    vertices = []
+    with open(vertexmap_filename) as fp:
+        for line in fp.xreadlines():
+            if line.startswith('#'): continue
+            vertices.append(map(int, line.split())[1])
+
+    dimension = max((max(d[1], d[2]) for d in boundary_list))
+    dimension += 1
+
+    # NB: D is a coboundary matrix; 1 and 2 below are transposed
+    D = spmatrix([d[0] for d in boundary_list],
+                 [d[2] for d in boundary_list],
+                 [d[1] for d in boundary_list], (dimension, dimension))
+
+           
+    z = spmatrix([zz[0] for zz in cocycle_list],
+                 [zz[1] for zz in cocycle_list],
+                 [0     for zz in cocycle_list], (dimension, 1))
+
+    v1 = D * z
+    # print "D^2 is zero:", not bool(D*D)
+    # print "D*z is zero:", not bool(v1)
+    z = matrix(z)
+
+    def Dfun(x,y,trans = 'N'):
+        if trans == 'N':
+            copy(D * x, y)
+        elif trans == 'T':
+            copy(D.T * x, y)
+        else:
+            assert False, "Unexpected trans parameter"
+
+    tol = 1e-10
+    show = False
+    maxit = None
+    solution = lsqr(Dfun, matrix(z), show = show, atol = tol, btol = tol, itnlim = maxit)
+    
+    v = z - D*solution[0]
+
+    # print sum(v**2)
+    # assert sum((D*v)**2) < tol and sum((D.T*v)**2) < tol, "Expected a harmonic cocycle"
+    if not (sum((D*v)**2) < tol and sum((D.T*v)**2) < tol):
+        print "Expected a harmonic cocycle:", sum((D*v)**2), sum((D.T*v)**2) 
+
+    values = [None]*len(vertices)
+    for i in xrange(len(vertices)):
+        values[vertices[i]] = solution[0][i]
+    return values
+
+
+if __name__ == '__main__':
+    if len(argv) < 4:
+        print "Usage: %s BOUNDARY COCYCLE VERTEXMAP" % argv[0]
+        exit()
+
+    boundary_filename = argv[1]
+    cocycle_filename = argv[2]
+    vertexmap_filename = argv[3]
+    values = smooth(boundary_filename, cocycle_filename, vertexmap_filename)
+
+    outfn = os.path.splitext(cocycle_filename)[0] + '.val'
+    with open(outfn, 'w') as f:
+        for i,v in enumerate(values):
+            f.write('%d %f\n' % (i,v))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/lsqr.py	Thu Jul 09 02:42:47 2009 -0700
@@ -0,0 +1,408 @@
+# LSQR solver from http://pages.cs.wisc.edu/~kline/cvxopt/
+
+from cvxopt import matrix
+from cvxopt.lapack import *
+from cvxopt.blas import *
+from math import sqrt
+
+"""
+a,b are scalars
+
+On exit, returns scalars c,s,r
+"""
+def SymOrtho(a,b):
+    aa=abs(a)
+    ab=abs(b)
+    if b==0.:
+        s=0.
+        r=aa
+        if aa==0.:
+            c=1.
+        else:
+            c=a/aa
+    elif a==0.:
+        c=0.
+        s=b/ab
+        r=ab
+    elif ab>aa:
+        sb=1
+        if b<0: sb=-1
+        tau=a/b
+        s=sb*(1+tau**2)**-0.5
+        c=s*tau
+        r=b/s
+    elif aa>ab:
+        sa=1
+        if a<0: sa=-1
+        tau=b/a
+        c=sa*(1+tau**2)**-0.5
+        s=c*tau
+        r=a/c
+        
+    return c,s,r
+
+"""
+
+It is usually recommended to use SYMMLQ for symmetric matrices
+
+Requires the syntax
+                   A(x,y)   == y:=[A]*x
+and
+           A(x,y,trans='T') == y:=[A.T]*x
+
+comments with '###' are followed by the intent of the original matlab
+code. This may be useful for debugging.
+
+"""
+
+def lsqr(  A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8, itnlim=None, show=False, wantvar=False):
+    """
+    
+    [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
+     = lsqr( m, n,  'aprod',  iw, rw, b, damp, atol, btol, conlim, itnlim, show );
+    
+     LSQR solves  Ax = b  or  min ||b - Ax||_2  if damp = 0,
+     or   min || (b)  -  (  A   )x ||   otherwise.
+              || (0)     (damp I)  ||2
+     A  is an m by n matrix defined by  y = aprod( mode,m,n,x,iw,rw ),
+     where the parameter 'aprodname' refers to a function 'aprod' that
+     performs the matrix-vector operations.
+     If mode = 1,   aprod  must return  y = Ax   without altering x.
+     If mode = 2,   aprod  must return  y = A'x  without altering x.
+     WARNING:   The file containing the function 'aprod'
+                must not be called aprodname.m !!!!
+
+    -----------------------------------------------------------------------
+     LSQR uses an iterative (conjugate-gradient-like) method.
+     For further information, see 
+     1. C. C. Paige and M. A. Saunders (1982a).
+        LSQR: An algorithm for sparse linear equations and sparse least squares,
+        ACM TOMS 8(1), 43-71.
+     2. C. C. Paige and M. A. Saunders (1982b).
+        Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
+        ACM TOMS 8(2), 195-209.
+     3. M. A. Saunders (1995).  Solution of sparse rectangular systems using
+        LSQR and CRAIG, BIT 35, 588-604.
+    
+     Input parameters:
+     iw, rw      are not used by lsqr, but are passed to aprod.
+     atol, btol  are stopping tolerances.  If both are 1.0e-9 (say),
+                 the final residual norm should be accurate to about 9 digits.
+                 (The final x will usually have fewer correct digits,
+                 depending on cond(A) and the size of damp.)
+     conlim      is also a stopping tolerance.  lsqr terminates if an estimate
+                 of cond(A) exceeds conlim.  For compatible systems Ax = b,
+                 conlim could be as large as 1.0e+12 (say).  For least-squares
+                 problems, conlim should be less than 1.0e+8.
+                 Maximum precision can be obtained by setting
+                 atol = btol = conlim = zero, but the number of iterations
+                 may then be excessive.
+     itnlim      is an explicit limit on iterations (for safety).
+     show = 1    gives an iteration log,
+     show = 0    suppresses output.
+    
+     Output parameters:
+     x           is the final solution.
+     istop       gives the reason for termination.
+     istop       = 1 means x is an approximate solution to Ax = b.
+                 = 2 means x approximately solves the least-squares problem.
+     r1norm      = norm(r), where r = b - Ax.
+     r2norm      = sqrt( norm(r)^2  +  damp^2 * norm(x)^2 )
+                 = r1norm if damp = 0.
+     anorm       = estimate of Frobenius norm of Abar = [  A   ].
+                                                        [damp*I]
+     acond       = estimate of cond(Abar).
+     arnorm      = estimate of norm(A'*r - damp^2*x).
+     xnorm       = norm(x).
+     var         (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
+                 or more generally (A'A + damp^2*I)^{-1}.
+                 This is well defined if A has full column rank or damp > 0.
+                 (Not sure what var means if rank(A) < n and damp = 0.)
+                 
+    
+            1990: Derived from Fortran 77 version of LSQR.
+     22 May 1992: bbnorm was used incorrectly.  Replaced by anorm.
+     26 Oct 1992: More input and output parameters added.
+     01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
+                  Print log reformatted.
+     14 Jun 1997: show  added to allow printing or not.
+     30 Jun 1997: var   added as an optional output parameter.
+     07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
+                  Michael Saunders, Systems Optimization Laboratory,
+                  Dept of MS&E, Stanford University.
+    -----------------------------------------------------------------------
+    """
+    """
+         Initialize.
+    """
+    n=len(b)
+    m=n
+    if itnlim is None: itnlim=2*n    
+
+    msg=('The exact solution is  x = 0                              ',
+         'Ax - b is small enough, given atol, btol                  ',
+         'The least-squares solution is good enough, given atol     ',
+         'The estimate of cond(Abar) has exceeded conlim            ',
+         'Ax - b is small enough for this machine                   ',
+         'The least-squares solution is good enough for this machine',
+         'Cond(Abar) seems to be too large for this machine         ',
+         'The iteration limit has been reached                      ');
+
+    var = matrix(0.,(n,1));
+
+    if show:
+        print ' '
+        print 'LSQR            Least-squares solution of  Ax = b'
+        str1 = 'The matrix A has %8g rows  and %8g cols' % (m, n)
+        str2 = 'damp = %20.14e    wantvar = %8g' %( damp,wantvar)
+        str3 = 'atol = %8.2e                 conlim = %8.2e'%( atol, conlim)
+        str4 = 'btol = %8.2e                 itnlim = %8g'  %( btol, itnlim)
+        print str1
+        print str2
+        print str3
+        print str4
+
+    itn    = 0;		istop  = 0;		nstop  = 0;
+    ctol   = 0;
+    if conlim > 0: ctol = 1/conlim
+    anorm  = 0;		acond  = 0;
+    dampsq = damp**2;	ddnorm = 0;		res2   = 0;
+    xnorm  = 0;		xxnorm = 0;		z      = 0;
+    cs2    = -1;		sn2    = 0;
+    
+    """
+    Set up the first vectors u and v for the bidiagonalization.
+     These satisfy  beta*u = b,  alfa*v = A'u.
+    """
+    __x    = matrix(0., (n,1)) # a matrix for temporary holding
+    v      = matrix(0., (n,1))
+    u      = +b;	
+    x      = matrix(0., (n,1))
+    alfa   = 0;
+    beta = nrm2( u );
+    w      = matrix(0., (n,1))
+    
+    if beta > 0:
+        ### u = (1/beta) * u;
+        ### v = feval( aprodname, 2, m, n, u, iw, rw );
+        scal(1/beta,u)
+	A(u,v,trans='T'); #v = feval( aprodname, 2, m, n, u, iw, rw );
+        alfa = nrm2( v );
+
+    if alfa > 0:
+        ### v = (1/alfa) * v;
+        scal(1/alfa,v)
+        copy(v,w)
+
+
+    rhobar = alfa;		phibar = beta;		bnorm  = beta;
+    rnorm  = beta;
+    r1norm = rnorm;
+    r2norm = rnorm;
+
+    # reverse the order here from the original matlab code because
+    # there was an error on return when arnorm==0
+    arnorm = alfa * beta;
+    if arnorm == 0:
+        print msg[0];
+        return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var 
+
+    head1  = '   Itn      x[0]       r1norm     r2norm ';
+    head2  = ' Compatible   LS      Norm A   Cond A';
+    
+    if show:
+        print ' '
+        print head1, head2
+        test1  = 1;		test2  = alfa / beta;
+        str1   = '%6g %12.5e'    %(    itn,   x[0] );
+        str2   = ' %10.3e %10.3e'%( r1norm, r2norm );
+        str3   = '  %8.1e %8.1e' %(  test1,  test2 );
+        print str1, str2, str3
+        
+    """
+    %------------------------------------------------------------------
+    %     Main iteration loop.
+    %------------------------------------------------------------------
+    """
+    while itn < itnlim:
+        itn = itn + 1;
+        """
+        %     Perform the next step of the bidiagonalization to obtain the
+        %     next  beta, u, alfa, v.  These satisfy the relations
+        %                beta*u  =  a*v   -  alfa*u,
+        %                alfa*v  =  A'*u  -  beta*v.
+        """
+        ### u    = feval( aprodname, 1, m, n, v, iw, rw )  -  alfa*u;
+        copy(u, __x)
+        A(v,u)
+        axpy(__x,u,-alfa)
+
+        beta = nrm2( u );
+        if beta > 0:
+            ### u     = (1/beta) * u;
+            scal(1/beta,u)
+            anorm = sqrt(anorm**2 + alfa**2 + beta**2 + damp**2);
+            ### v     = feval( aprodname, 2, m, n, u, iw, rw )  -  beta*v;
+            copy(v,__x)
+            A(u,v,trans='T')
+            axpy(__x,v,-beta)
+
+            alfa  = nrm2( v );
+            if alfa > 0:
+                ### v = (1/alfa) * v;
+                scal(1/alfa, v)
+
+        """
+        %     Use a plane rotation to eliminate the damping parameter.
+        %     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
+        """
+
+        rhobar1 = sqrt(rhobar**2 + damp**2);
+        cs1     = rhobar / rhobar1;
+        sn1     = damp   / rhobar1;
+        psi     = sn1 * phibar;
+        phibar  = cs1 * phibar;
+        """
+        %     Use a plane rotation to eliminate the subdiagonal element (beta)
+        %     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
+        """
+
+
+        ###cs      =   rhobar1/ rho;
+        ###sn      =   beta   / rho;
+        cs,sn,rho = SymOrtho(rhobar1,beta)
+        
+        theta   =   sn * alfa;
+        rhobar  = - cs * alfa;
+        phi     =   cs * phibar;
+        phibar  =   sn * phibar;
+        tau     =   sn * phi;
+        """
+        %     Update x and w.
+        """
+        t1      =   phi  /rho;
+        t2      = - theta/rho;
+        dk      =   (1/rho)*w;
+
+        ### x       = x      +  t1*w;
+        axpy(w,x,t1)
+        ### w       = v      +  t2*w;
+        scal(t2,w)
+        axpy(v,w)
+        ddnorm  = ddnorm +  nrm2(dk)**2;
+        if wantvar:
+            ### var = var  +  dk.*dk; 
+            axpy(dk**2, var)
+        """
+        %     Use a plane rotation on the right to eliminate the
+        %     super-diagonal element (theta) of the upper-bidiagonal matrix.
+        %     Then use the result to estimate  norm(x).
+        """
+
+        delta   =   sn2 * rho;
+        gambar  = - cs2 * rho;
+        rhs     =   phi  -  delta * z;
+        zbar    =   rhs / gambar;
+        xnorm   =   sqrt(xxnorm + zbar**2);
+        gamma   =   sqrt(gambar**2 +theta**2);
+        cs2     =   gambar / gamma;
+        sn2     =   theta  / gamma;
+        z       =   rhs    / gamma;
+        xxnorm  =   xxnorm  +  z**2;
+        """
+        %     Test for convergence.
+        %     First, estimate the condition of the matrix  Abar,
+        %     and the norms of  rbar  and  Abar'rbar.
+        """
+        acond   =   anorm * sqrt(ddnorm);
+        res1    =   phibar**2;
+        res2    =   res2  +  psi**2;
+        rnorm   =   sqrt( res1 + res2 );
+        arnorm  =   alfa * abs( tau );
+        """
+        %     07 Aug 2002:
+        %     Distinguish between
+        %        r1norm = ||b - Ax|| and
+        %        r2norm = rnorm in current code
+        %               = sqrt(r1norm^2 + damp^2*||x||^2).
+        %        Estimate r1norm from
+        %        r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
+        %     Although there is cancellation, it might be accurate enough.
+        """
+        r1sq    =   rnorm**2  -  dampsq * xxnorm;
+        r1norm  =   sqrt( abs(r1sq) );
+        if r1sq < 0: r1norm = - r1norm; 
+        r2norm  =   rnorm;
+        """
+        %     Now use these norms to estimate certain other quantities,
+        %     some of which will be small near a solution.
+        """
+        test1   =   rnorm / bnorm;
+        test2   =   arnorm/( anorm * rnorm );
+        test3   =       1 / acond;
+        t1      =   test1 / (1    +  anorm * xnorm / bnorm);
+        rtol    =   btol  +  atol *  anorm * xnorm / bnorm;
+        """
+        %     The following tests guard against extremely small values of
+        %     atol, btol  or  ctol.  (The user may have set any or all of
+        %     the parameters  atol, btol, conlim  to 0.)
+        %     The effect is equivalent to the normal tests using
+        %     atol = eps,  btol = eps,  conlim = 1/eps.
+        """
+        if itn >= itnlim  : istop = 7; 
+        if 1 + test3  <= 1: istop = 6; 
+        if 1 + test2  <= 1: istop = 5; 
+        if 1 + t1     <= 1: istop = 4; 
+        """
+        %     Allow for tolerances set by the user.
+        """
+        if  test3 <= ctol:  istop = 3;
+        if  test2 <= atol:  istop = 2;
+        if  test1 <= rtol:  istop = 1;
+        """
+        %     See if it is time to print something.
+        """
+        prnt = False;
+        if n     <= 40       : prnt = True;
+        if itn   <= 10       : prnt = True;
+        if itn   >= itnlim-10: prnt = True;
+        # if itn%10 == 0       : prnt = True;
+        if test3 <=  2*ctol  : prnt = True;
+        if test2 <= 10*atol  : prnt = True;
+        if test1 <= 10*rtol  : prnt = True;
+        if istop !=  0       : prnt = True;
+
+        if prnt:
+            if show:
+                str1 = '%6g %12.5e'%        (itn,   x[0] );
+                str2 = ' %10.3e %10.3e'% (r1norm, r2norm );
+                str3 = '  %8.1e %8.1e'%  ( test1,  test2 );
+                str4 = ' %8.1e %8.1e'%   ( anorm,  acond );
+                print str1, str2, str3, str4
+
+        if istop != 0: break
+
+    """
+    %     End of iteration loop.
+    %     Print the stopping condition.
+    """
+    if show:
+        print ' '
+        print 'LSQR finished'
+        print msg[istop]
+        print ' '
+        str1 = 'istop =%8g   r1norm =%8.1e'%  ( istop, r1norm );
+        str2 = 'anorm =%8.1e   arnorm =%8.1e'%( anorm, arnorm );
+        str3 = 'itn   =%8g   r2norm =%8.1e'%  (   itn, r2norm );
+        str4 = 'acond =%8.1e   xnorm  =%8.1e'%( acond, xnorm  );
+        print str1+ '   ' +str2
+        print str3+ '   ' +str4
+        print ' '
+
+    return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var 
+            
+    """
+    %-----------------------------------------------------------------------
+    % End of lsqr.m
+    %-----------------------------------------------------------------------
+    """
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/output.h	Thu Jul 09 02:42:47 2009 -0700
@@ -0,0 +1,64 @@
+#include <iostream>
+#include <sstream>
+#include <fstream>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/tuple/tuple_comparison.hpp>
+
+#include <cassert>
+
+bool neq(const Smplx& s1, const Smplx& s2)               
+{ 
+    Smplx::VertexComparison cmp;
+    return cmp(s1, s2) || cmp(s2, s1);
+}
+
+unsigned index(const SimplexVector& v, const Smplx& s, const Generator::Comparison& cmp)
+{
+    SimplexVector::const_iterator it = std::lower_bound(v.begin(), v.end(), s, cmp);
+    while (neq(*it, s)) ++it;
+    return it - v.begin();
+}
+
+void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Generator::Comparison& cmp)
+{
+    unsigned i = 0;
+    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+    {
+        // std::cout << "Simplex: " << *cur << std::endl;
+        bool                sign = true;
+        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
+        {
+            // std::cout << "  " << *bcur << std::endl;
+            out << (sign ? 1 : -1) << " ";
+            out << index(v, *bcur, cmp) << " " << i << "\n";
+            sign = !sign;
+        }
+        ++i;
+    }
+}
+
+void output_vertex_indices(std::ostream& out, const SimplexVector& v)
+{
+    unsigned i = 0;
+    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+    {
+        if (cur->dimension() == 0)
+            out << i << " " << cur->vertices()[0] << std::endl;
+        ++i;
+    }
+}
+
+void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime, const Generator::Comparison& cmp)
+{
+    std::ostringstream istr; istr << '-' << i;
+    std::string filename = cocycle_prefix + istr.str() + ".ccl";
+    std::ofstream out(filename.c_str());
+    out << "# Cocycle born at " << c.birth.get<1>() << std::endl;
+    for (Persistence::ZColumn::const_iterator zcur = c.zcolumn.begin(); zcur != c.zcolumn.end(); ++zcur)
+    {
+        const Smplx& s = **(zcur->si);
+        out << (zcur->coefficient <= prime/2 ? zcur->coefficient : zcur->coefficient - prime) << " ";
+        out << index(v, s, cmp) << "\n";
+    }
+}
--- a/examples/cohomology/rips-pairwise-cohomology.cpp	Thu Jul 09 00:59:32 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp	Thu Jul 09 02:42:47 2009 -0700
@@ -28,7 +28,9 @@
 typedef     Persistence::Death                                      Death;
 typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
 
-void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime);
+#include "output.h"         // for output_*()
+
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name);
 
 int main(int argc, char* argv[])
 {
@@ -41,9 +43,16 @@
     Dimension               skeleton;
     DistanceType            max_distance;
     ZpField::Element        prime;
-    std::string             infilename;
+    std::string             infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
 
-    program_options(argc, argv, infilename, skeleton, max_distance, prime);
+    program_options(argc, argv, infilename, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
+    std::ofstream           bdry_out(boundary_name.c_str());
+    std::ofstream           vertices_out(vertices_name.c_str());
+    std::ofstream           diagram_out(diagram_name.c_str());
+    std::cout << "Boundary matrix: " << boundary_name << std::endl;
+    std::cout << "Cocycles:        " << cocycle_prefix << "*.ccl" << std::endl;
+    std::cout << "Vertices:        " << vertices_name << std::endl;
+    std::cout << "Diagram:         " << diagram_name << std::endl;
 
     Timer total_timer; total_timer.start();
     PointContainer          points;
@@ -52,15 +61,19 @@
     PairDistances           distances(points);
     Generator               rips(distances);
     Generator::Evaluator    size(distances);
+    Generator::Comparison   cmp(distances);
     SimplexVector           v;
     Complex                 c;
     
     Timer rips_timer; rips_timer.start();
     rips.generate(skeleton, max_distance, make_push_back_functor(v));
-    std::sort(v.begin(), v.end(), Generator::Comparison(distances));
+    std::sort(v.begin(), v.end(), cmp);
     rips_timer.stop();
     std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
 
+    output_boundary_matrix(bdry_out, v, cmp);
+    output_vertex_indices(vertices_out, v);
+
     Timer persistence_timer; persistence_timer.start();
     ZpField                 zp(prime);
     Persistence             p(zp);
@@ -81,28 +94,24 @@
         if (d && (size(*cur) - d->get<1>()) > 0)
         {
             AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
-            std::cout << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+            diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
         }
     }
-    // output infinte persistence cocycles
+    // output infinte persistence pairs 
     for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
-        std::cout << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
+        diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
     persistence_timer.stop();
 
 
     // p.show_cocycles();
-    // Output alive cocycles
+    // Output alive cocycles of dimension 1
+    unsigned i = 0;
     for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
     {
-        std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
-        for (Persistence::ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
-        {
-            const Smplx& s = **(zcur->si);
-            std::cout << zcur->coefficient << " ";
-            for (Smplx::VertexContainer::const_iterator vcur = s.vertices().begin(); vcur != s.vertices().end(); ++vcur)
-                std::cout << *vcur << " ";
-            std::cout << std::endl;
-        }
+        if (cur->birth.get<0>() != 1) continue;
+        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
+        // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+        ++i;
     }
     total_timer.stop();
     rips_timer.check("Rips timer");
@@ -110,7 +119,7 @@
     total_timer.check("Total timer");
 }
 
-void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime)
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name)
 {
     namespace po = boost::program_options;
 
@@ -122,8 +131,12 @@
     visible.add_options()
         ("help,h",                                                                                  "produce help message")
         ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2),                  "Dimension of the Rips complex we want to compute")
-        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(2),              "Prime p for the field F_p")
-        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction");
+        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(11),             "Prime p for the field F_p")
+        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction")
+        ("boundary,b",          po::value<std::string>(&boundary_name),                             "Filename where to output the boundary matrix")
+        ("cocycle,c",           po::value<std::string>(&cocycle_prefix),                            "Prefix of the filename where to output the 1-dimensional cocycles")
+        ("vertices,v",          po::value<std::string>(&vertices_name),                             "Filename where to output the simplex-vertex mapping")
+        ("diagram,d",           po::value<std::string>(&diagram_name),                              "Filename where to output the persistence diagram");
 #if LOGGING
     std::vector<std::string>    log_channels;
     visible.add_options()
@@ -132,7 +145,6 @@
 
     po::positional_options_description pos;
     pos.add("input-file", 1);
-    pos.add("output-file", 2);
     
     po::options_description all; all.add(visible).add(hidden);
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/plot-values/plot.py	Thu Jul 09 02:42:47 2009 -0700
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+from    pylab           import scatter, show, cm, colorbar, savefig, axis, \
+                               figure, xlim, axes, hsv, subplots_adjust as adjust
+from    itertools       import izip
+from    sys             import argv, exit
+import  os.path         as     osp
+
+
+def plot(val_fn, pts_fn, output_fn):
+    points = []
+    with open(pts_fn) as fp:
+        for line in fp.xreadlines():
+            points.append(map(float, line.split()))
+    
+    values = []
+    with open(val_fn) as fp:
+        for line in fp.xreadlines():
+            values.append(float(line.split()[1]))
+
+    xx = [pt[0] for pt in points]
+    yy = [pt[1] for pt in points]
+    print "X:", min(xx), max(xx)
+    print "Y:", min(yy), max(yy)
+
+    m = min(values)
+    values = [(v-m) % 1. for v in values]
+    print "V:", min(values), max(values)
+
+    aspect = (max(yy) - min(yy))/(max(xx) - min(xx)) + .1
+    # aspect = .5
+
+    # hsv()
+    # fig = figure(figsize = (3,3*aspect))
+    # ax = fig.add_axes([-.05,-.1,1.1,1.1])
+    ax = axes()
+    ax.set_axis_off()
+    ax.set_aspect('equal', 'box')
+    ax.scatter(xx,yy,s=10,c=values)
+    # adjust(0,0,1,1,0,0)
+    colorbar()
+    fig.savefig(output_fn)
+
+if __name__ == '__main__':
+    if len(argv) < 3:
+        print "Usage: %s VALUES POINTS" % argv[0]
+        exit()
+
+    val_fn = argv[1]
+    pts_fn  = argv[2]
+    output_fn, ext = osp.splitext(val_fn)
+    output_fn += '.pdf'
+    plot(val_fn, pts_fn, output_fn)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/plot-values/scatter.py	Thu Jul 09 02:42:47 2009 -0700
@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+
+from    pylab           import scatter, show, cm, colorbar, axes, savefig
+from    itertools       import izip
+from    sys             import argv, exit
+import  os.path         as     osp
+
+
+def plot(val1_fn, val2_fn):
+    values1 = []
+    with open(val1_fn) as fp:
+        for line in fp.xreadlines():
+            values1.append(float(line.split()[1]))
+    
+    values2 = []
+    with open(val2_fn) as fp:
+        for line in fp.xreadlines():
+            values2.append(float(line.split()[1]))
+    
+    values1 = [v % 1. for v in values1]
+    values2 = [v % 1. for v in values2]
+    print min(values1), max(values2), min(values1), min(values2)
+
+    scatter(values1, values2)
+    axes().set_aspect('equal')
+    show()
+
+if __name__ == '__main__':
+    if len(argv) < 3:
+        print "Usage: %s VALUES1 VALUES2" % argv[0]
+        exit()
+
+    val1_fn = argv[1]
+    val2_fn  = argv[2]
+    plot(val1_fn, val2_fn)