Rips pairwise cohomology produces output necessary for 1-cocycle parametrization dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu Jul 09 02:42:47 2009 -0700 (2009-07-09)
branchdev
changeset 13896030f8d2f2c
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
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/examples/cohomology/cocycle.py	Thu Jul 09 02:42:47 2009 -0700
     1.3 @@ -0,0 +1,85 @@
     1.4 +#!/usr/bin/env python
     1.5 +
     1.6 +from    cvxopt          import spmatrix, matrix
     1.7 +from    cvxopt.blas     import copy
     1.8 +from    lsqr            import lsqr
     1.9 +from    sys             import argv, exit
    1.10 +import  os.path
    1.11 +
    1.12 +def smooth(boundary_filename, cocycle_filename, vertexmap_filename):
    1.13 +    boundary_list = []
    1.14 +    with open(boundary_filename) as fp:
    1.15 +        for line in fp.xreadlines():
    1.16 +            if line.startswith('#'): continue
    1.17 +            boundary_list.append(map(int, line.split()))
    1.18 +
    1.19 +    cocycle_list = []
    1.20 +    with open(cocycle_filename) as fp:
    1.21 +        for line in fp.xreadlines():
    1.22 +            if line.startswith('#'): continue
    1.23 +            cocycle_list.append(map(int, line.split()))
    1.24 +
    1.25 +    vertices = []
    1.26 +    with open(vertexmap_filename) as fp:
    1.27 +        for line in fp.xreadlines():
    1.28 +            if line.startswith('#'): continue
    1.29 +            vertices.append(map(int, line.split())[1])
    1.30 +
    1.31 +    dimension = max((max(d[1], d[2]) for d in boundary_list))
    1.32 +    dimension += 1
    1.33 +
    1.34 +    # NB: D is a coboundary matrix; 1 and 2 below are transposed
    1.35 +    D = spmatrix([d[0] for d in boundary_list],
    1.36 +                 [d[2] for d in boundary_list],
    1.37 +                 [d[1] for d in boundary_list], (dimension, dimension))
    1.38 +
    1.39 +           
    1.40 +    z = spmatrix([zz[0] for zz in cocycle_list],
    1.41 +                 [zz[1] for zz in cocycle_list],
    1.42 +                 [0     for zz in cocycle_list], (dimension, 1))
    1.43 +
    1.44 +    v1 = D * z
    1.45 +    # print "D^2 is zero:", not bool(D*D)
    1.46 +    # print "D*z is zero:", not bool(v1)
    1.47 +    z = matrix(z)
    1.48 +
    1.49 +    def Dfun(x,y,trans = 'N'):
    1.50 +        if trans == 'N':
    1.51 +            copy(D * x, y)
    1.52 +        elif trans == 'T':
    1.53 +            copy(D.T * x, y)
    1.54 +        else:
    1.55 +            assert False, "Unexpected trans parameter"
    1.56 +
    1.57 +    tol = 1e-10
    1.58 +    show = False
    1.59 +    maxit = None
    1.60 +    solution = lsqr(Dfun, matrix(z), show = show, atol = tol, btol = tol, itnlim = maxit)
    1.61 +    
    1.62 +    v = z - D*solution[0]
    1.63 +
    1.64 +    # print sum(v**2)
    1.65 +    # assert sum((D*v)**2) < tol and sum((D.T*v)**2) < tol, "Expected a harmonic cocycle"
    1.66 +    if not (sum((D*v)**2) < tol and sum((D.T*v)**2) < tol):
    1.67 +        print "Expected a harmonic cocycle:", sum((D*v)**2), sum((D.T*v)**2) 
    1.68 +
    1.69 +    values = [None]*len(vertices)
    1.70 +    for i in xrange(len(vertices)):
    1.71 +        values[vertices[i]] = solution[0][i]
    1.72 +    return values
    1.73 +
    1.74 +
    1.75 +if __name__ == '__main__':
    1.76 +    if len(argv) < 4:
    1.77 +        print "Usage: %s BOUNDARY COCYCLE VERTEXMAP" % argv[0]
    1.78 +        exit()
    1.79 +
    1.80 +    boundary_filename = argv[1]
    1.81 +    cocycle_filename = argv[2]
    1.82 +    vertexmap_filename = argv[3]
    1.83 +    values = smooth(boundary_filename, cocycle_filename, vertexmap_filename)
    1.84 +
    1.85 +    outfn = os.path.splitext(cocycle_filename)[0] + '.val'
    1.86 +    with open(outfn, 'w') as f:
    1.87 +        for i,v in enumerate(values):
    1.88 +            f.write('%d %f\n' % (i,v))
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/examples/cohomology/lsqr.py	Thu Jul 09 02:42:47 2009 -0700
     2.3 @@ -0,0 +1,408 @@
     2.4 +# LSQR solver from http://pages.cs.wisc.edu/~kline/cvxopt/
     2.5 +
     2.6 +from cvxopt import matrix
     2.7 +from cvxopt.lapack import *
     2.8 +from cvxopt.blas import *
     2.9 +from math import sqrt
    2.10 +
    2.11 +"""
    2.12 +a,b are scalars
    2.13 +
    2.14 +On exit, returns scalars c,s,r
    2.15 +"""
    2.16 +def SymOrtho(a,b):
    2.17 +    aa=abs(a)
    2.18 +    ab=abs(b)
    2.19 +    if b==0.:
    2.20 +        s=0.
    2.21 +        r=aa
    2.22 +        if aa==0.:
    2.23 +            c=1.
    2.24 +        else:
    2.25 +            c=a/aa
    2.26 +    elif a==0.:
    2.27 +        c=0.
    2.28 +        s=b/ab
    2.29 +        r=ab
    2.30 +    elif ab>aa:
    2.31 +        sb=1
    2.32 +        if b<0: sb=-1
    2.33 +        tau=a/b
    2.34 +        s=sb*(1+tau**2)**-0.5
    2.35 +        c=s*tau
    2.36 +        r=b/s
    2.37 +    elif aa>ab:
    2.38 +        sa=1
    2.39 +        if a<0: sa=-1
    2.40 +        tau=b/a
    2.41 +        c=sa*(1+tau**2)**-0.5
    2.42 +        s=c*tau
    2.43 +        r=a/c
    2.44 +        
    2.45 +    return c,s,r
    2.46 +
    2.47 +"""
    2.48 +
    2.49 +It is usually recommended to use SYMMLQ for symmetric matrices
    2.50 +
    2.51 +Requires the syntax
    2.52 +                   A(x,y)   == y:=[A]*x
    2.53 +and
    2.54 +           A(x,y,trans='T') == y:=[A.T]*x
    2.55 +
    2.56 +comments with '###' are followed by the intent of the original matlab
    2.57 +code. This may be useful for debugging.
    2.58 +
    2.59 +"""
    2.60 +
    2.61 +def lsqr(  A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8, itnlim=None, show=False, wantvar=False):
    2.62 +    """
    2.63 +    
    2.64 +    [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
    2.65 +     = lsqr( m, n,  'aprod',  iw, rw, b, damp, atol, btol, conlim, itnlim, show );
    2.66 +    
    2.67 +     LSQR solves  Ax = b  or  min ||b - Ax||_2  if damp = 0,
    2.68 +     or   min || (b)  -  (  A   )x ||   otherwise.
    2.69 +              || (0)     (damp I)  ||2
    2.70 +     A  is an m by n matrix defined by  y = aprod( mode,m,n,x,iw,rw ),
    2.71 +     where the parameter 'aprodname' refers to a function 'aprod' that
    2.72 +     performs the matrix-vector operations.
    2.73 +     If mode = 1,   aprod  must return  y = Ax   without altering x.
    2.74 +     If mode = 2,   aprod  must return  y = A'x  without altering x.
    2.75 +     WARNING:   The file containing the function 'aprod'
    2.76 +                must not be called aprodname.m !!!!
    2.77 +
    2.78 +    -----------------------------------------------------------------------
    2.79 +     LSQR uses an iterative (conjugate-gradient-like) method.
    2.80 +     For further information, see 
    2.81 +     1. C. C. Paige and M. A. Saunders (1982a).
    2.82 +        LSQR: An algorithm for sparse linear equations and sparse least squares,
    2.83 +        ACM TOMS 8(1), 43-71.
    2.84 +     2. C. C. Paige and M. A. Saunders (1982b).
    2.85 +        Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
    2.86 +        ACM TOMS 8(2), 195-209.
    2.87 +     3. M. A. Saunders (1995).  Solution of sparse rectangular systems using
    2.88 +        LSQR and CRAIG, BIT 35, 588-604.
    2.89 +    
    2.90 +     Input parameters:
    2.91 +     iw, rw      are not used by lsqr, but are passed to aprod.
    2.92 +     atol, btol  are stopping tolerances.  If both are 1.0e-9 (say),
    2.93 +                 the final residual norm should be accurate to about 9 digits.
    2.94 +                 (The final x will usually have fewer correct digits,
    2.95 +                 depending on cond(A) and the size of damp.)
    2.96 +     conlim      is also a stopping tolerance.  lsqr terminates if an estimate
    2.97 +                 of cond(A) exceeds conlim.  For compatible systems Ax = b,
    2.98 +                 conlim could be as large as 1.0e+12 (say).  For least-squares
    2.99 +                 problems, conlim should be less than 1.0e+8.
   2.100 +                 Maximum precision can be obtained by setting
   2.101 +                 atol = btol = conlim = zero, but the number of iterations
   2.102 +                 may then be excessive.
   2.103 +     itnlim      is an explicit limit on iterations (for safety).
   2.104 +     show = 1    gives an iteration log,
   2.105 +     show = 0    suppresses output.
   2.106 +    
   2.107 +     Output parameters:
   2.108 +     x           is the final solution.
   2.109 +     istop       gives the reason for termination.
   2.110 +     istop       = 1 means x is an approximate solution to Ax = b.
   2.111 +                 = 2 means x approximately solves the least-squares problem.
   2.112 +     r1norm      = norm(r), where r = b - Ax.
   2.113 +     r2norm      = sqrt( norm(r)^2  +  damp^2 * norm(x)^2 )
   2.114 +                 = r1norm if damp = 0.
   2.115 +     anorm       = estimate of Frobenius norm of Abar = [  A   ].
   2.116 +                                                        [damp*I]
   2.117 +     acond       = estimate of cond(Abar).
   2.118 +     arnorm      = estimate of norm(A'*r - damp^2*x).
   2.119 +     xnorm       = norm(x).
   2.120 +     var         (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
   2.121 +                 or more generally (A'A + damp^2*I)^{-1}.
   2.122 +                 This is well defined if A has full column rank or damp > 0.
   2.123 +                 (Not sure what var means if rank(A) < n and damp = 0.)
   2.124 +                 
   2.125 +    
   2.126 +            1990: Derived from Fortran 77 version of LSQR.
   2.127 +     22 May 1992: bbnorm was used incorrectly.  Replaced by anorm.
   2.128 +     26 Oct 1992: More input and output parameters added.
   2.129 +     01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
   2.130 +                  Print log reformatted.
   2.131 +     14 Jun 1997: show  added to allow printing or not.
   2.132 +     30 Jun 1997: var   added as an optional output parameter.
   2.133 +     07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
   2.134 +                  Michael Saunders, Systems Optimization Laboratory,
   2.135 +                  Dept of MS&E, Stanford University.
   2.136 +    -----------------------------------------------------------------------
   2.137 +    """
   2.138 +    """
   2.139 +         Initialize.
   2.140 +    """
   2.141 +    n=len(b)
   2.142 +    m=n
   2.143 +    if itnlim is None: itnlim=2*n    
   2.144 +
   2.145 +    msg=('The exact solution is  x = 0                              ',
   2.146 +         'Ax - b is small enough, given atol, btol                  ',
   2.147 +         'The least-squares solution is good enough, given atol     ',
   2.148 +         'The estimate of cond(Abar) has exceeded conlim            ',
   2.149 +         'Ax - b is small enough for this machine                   ',
   2.150 +         'The least-squares solution is good enough for this machine',
   2.151 +         'Cond(Abar) seems to be too large for this machine         ',
   2.152 +         'The iteration limit has been reached                      ');
   2.153 +
   2.154 +    var = matrix(0.,(n,1));
   2.155 +
   2.156 +    if show:
   2.157 +        print ' '
   2.158 +        print 'LSQR            Least-squares solution of  Ax = b'
   2.159 +        str1 = 'The matrix A has %8g rows  and %8g cols' % (m, n)
   2.160 +        str2 = 'damp = %20.14e    wantvar = %8g' %( damp,wantvar)
   2.161 +        str3 = 'atol = %8.2e                 conlim = %8.2e'%( atol, conlim)
   2.162 +        str4 = 'btol = %8.2e                 itnlim = %8g'  %( btol, itnlim)
   2.163 +        print str1
   2.164 +        print str2
   2.165 +        print str3
   2.166 +        print str4
   2.167 +
   2.168 +    itn    = 0;		istop  = 0;		nstop  = 0;
   2.169 +    ctol   = 0;
   2.170 +    if conlim > 0: ctol = 1/conlim
   2.171 +    anorm  = 0;		acond  = 0;
   2.172 +    dampsq = damp**2;	ddnorm = 0;		res2   = 0;
   2.173 +    xnorm  = 0;		xxnorm = 0;		z      = 0;
   2.174 +    cs2    = -1;		sn2    = 0;
   2.175 +    
   2.176 +    """
   2.177 +    Set up the first vectors u and v for the bidiagonalization.
   2.178 +     These satisfy  beta*u = b,  alfa*v = A'u.
   2.179 +    """
   2.180 +    __x    = matrix(0., (n,1)) # a matrix for temporary holding
   2.181 +    v      = matrix(0., (n,1))
   2.182 +    u      = +b;	
   2.183 +    x      = matrix(0., (n,1))
   2.184 +    alfa   = 0;
   2.185 +    beta = nrm2( u );
   2.186 +    w      = matrix(0., (n,1))
   2.187 +    
   2.188 +    if beta > 0:
   2.189 +        ### u = (1/beta) * u;
   2.190 +        ### v = feval( aprodname, 2, m, n, u, iw, rw );
   2.191 +        scal(1/beta,u)
   2.192 +	A(u,v,trans='T'); #v = feval( aprodname, 2, m, n, u, iw, rw );
   2.193 +        alfa = nrm2( v );
   2.194 +
   2.195 +    if alfa > 0:
   2.196 +        ### v = (1/alfa) * v;
   2.197 +        scal(1/alfa,v)
   2.198 +        copy(v,w)
   2.199 +
   2.200 +
   2.201 +    rhobar = alfa;		phibar = beta;		bnorm  = beta;
   2.202 +    rnorm  = beta;
   2.203 +    r1norm = rnorm;
   2.204 +    r2norm = rnorm;
   2.205 +
   2.206 +    # reverse the order here from the original matlab code because
   2.207 +    # there was an error on return when arnorm==0
   2.208 +    arnorm = alfa * beta;
   2.209 +    if arnorm == 0:
   2.210 +        print msg[0];
   2.211 +        return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var 
   2.212 +
   2.213 +    head1  = '   Itn      x[0]       r1norm     r2norm ';
   2.214 +    head2  = ' Compatible   LS      Norm A   Cond A';
   2.215 +    
   2.216 +    if show:
   2.217 +        print ' '
   2.218 +        print head1, head2
   2.219 +        test1  = 1;		test2  = alfa / beta;
   2.220 +        str1   = '%6g %12.5e'    %(    itn,   x[0] );
   2.221 +        str2   = ' %10.3e %10.3e'%( r1norm, r2norm );
   2.222 +        str3   = '  %8.1e %8.1e' %(  test1,  test2 );
   2.223 +        print str1, str2, str3
   2.224 +        
   2.225 +    """
   2.226 +    %------------------------------------------------------------------
   2.227 +    %     Main iteration loop.
   2.228 +    %------------------------------------------------------------------
   2.229 +    """
   2.230 +    while itn < itnlim:
   2.231 +        itn = itn + 1;
   2.232 +        """
   2.233 +        %     Perform the next step of the bidiagonalization to obtain the
   2.234 +        %     next  beta, u, alfa, v.  These satisfy the relations
   2.235 +        %                beta*u  =  a*v   -  alfa*u,
   2.236 +        %                alfa*v  =  A'*u  -  beta*v.
   2.237 +        """
   2.238 +        ### u    = feval( aprodname, 1, m, n, v, iw, rw )  -  alfa*u;
   2.239 +        copy(u, __x)
   2.240 +        A(v,u)
   2.241 +        axpy(__x,u,-alfa)
   2.242 +
   2.243 +        beta = nrm2( u );
   2.244 +        if beta > 0:
   2.245 +            ### u     = (1/beta) * u;
   2.246 +            scal(1/beta,u)
   2.247 +            anorm = sqrt(anorm**2 + alfa**2 + beta**2 + damp**2);
   2.248 +            ### v     = feval( aprodname, 2, m, n, u, iw, rw )  -  beta*v;
   2.249 +            copy(v,__x)
   2.250 +            A(u,v,trans='T')
   2.251 +            axpy(__x,v,-beta)
   2.252 +
   2.253 +            alfa  = nrm2( v );
   2.254 +            if alfa > 0:
   2.255 +                ### v = (1/alfa) * v;
   2.256 +                scal(1/alfa, v)
   2.257 +
   2.258 +        """
   2.259 +        %     Use a plane rotation to eliminate the damping parameter.
   2.260 +        %     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
   2.261 +        """
   2.262 +
   2.263 +        rhobar1 = sqrt(rhobar**2 + damp**2);
   2.264 +        cs1     = rhobar / rhobar1;
   2.265 +        sn1     = damp   / rhobar1;
   2.266 +        psi     = sn1 * phibar;
   2.267 +        phibar  = cs1 * phibar;
   2.268 +        """
   2.269 +        %     Use a plane rotation to eliminate the subdiagonal element (beta)
   2.270 +        %     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
   2.271 +        """
   2.272 +
   2.273 +
   2.274 +        ###cs      =   rhobar1/ rho;
   2.275 +        ###sn      =   beta   / rho;
   2.276 +        cs,sn,rho = SymOrtho(rhobar1,beta)
   2.277 +        
   2.278 +        theta   =   sn * alfa;
   2.279 +        rhobar  = - cs * alfa;
   2.280 +        phi     =   cs * phibar;
   2.281 +        phibar  =   sn * phibar;
   2.282 +        tau     =   sn * phi;
   2.283 +        """
   2.284 +        %     Update x and w.
   2.285 +        """
   2.286 +        t1      =   phi  /rho;
   2.287 +        t2      = - theta/rho;
   2.288 +        dk      =   (1/rho)*w;
   2.289 +
   2.290 +        ### x       = x      +  t1*w;
   2.291 +        axpy(w,x,t1)
   2.292 +        ### w       = v      +  t2*w;
   2.293 +        scal(t2,w)
   2.294 +        axpy(v,w)
   2.295 +        ddnorm  = ddnorm +  nrm2(dk)**2;
   2.296 +        if wantvar:
   2.297 +            ### var = var  +  dk.*dk; 
   2.298 +            axpy(dk**2, var)
   2.299 +        """
   2.300 +        %     Use a plane rotation on the right to eliminate the
   2.301 +        %     super-diagonal element (theta) of the upper-bidiagonal matrix.
   2.302 +        %     Then use the result to estimate  norm(x).
   2.303 +        """
   2.304 +
   2.305 +        delta   =   sn2 * rho;
   2.306 +        gambar  = - cs2 * rho;
   2.307 +        rhs     =   phi  -  delta * z;
   2.308 +        zbar    =   rhs / gambar;
   2.309 +        xnorm   =   sqrt(xxnorm + zbar**2);
   2.310 +        gamma   =   sqrt(gambar**2 +theta**2);
   2.311 +        cs2     =   gambar / gamma;
   2.312 +        sn2     =   theta  / gamma;
   2.313 +        z       =   rhs    / gamma;
   2.314 +        xxnorm  =   xxnorm  +  z**2;
   2.315 +        """
   2.316 +        %     Test for convergence.
   2.317 +        %     First, estimate the condition of the matrix  Abar,
   2.318 +        %     and the norms of  rbar  and  Abar'rbar.
   2.319 +        """
   2.320 +        acond   =   anorm * sqrt(ddnorm);
   2.321 +        res1    =   phibar**2;
   2.322 +        res2    =   res2  +  psi**2;
   2.323 +        rnorm   =   sqrt( res1 + res2 );
   2.324 +        arnorm  =   alfa * abs( tau );
   2.325 +        """
   2.326 +        %     07 Aug 2002:
   2.327 +        %     Distinguish between
   2.328 +        %        r1norm = ||b - Ax|| and
   2.329 +        %        r2norm = rnorm in current code
   2.330 +        %               = sqrt(r1norm^2 + damp^2*||x||^2).
   2.331 +        %        Estimate r1norm from
   2.332 +        %        r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
   2.333 +        %     Although there is cancellation, it might be accurate enough.
   2.334 +        """
   2.335 +        r1sq    =   rnorm**2  -  dampsq * xxnorm;
   2.336 +        r1norm  =   sqrt( abs(r1sq) );
   2.337 +        if r1sq < 0: r1norm = - r1norm; 
   2.338 +        r2norm  =   rnorm;
   2.339 +        """
   2.340 +        %     Now use these norms to estimate certain other quantities,
   2.341 +        %     some of which will be small near a solution.
   2.342 +        """
   2.343 +        test1   =   rnorm / bnorm;
   2.344 +        test2   =   arnorm/( anorm * rnorm );
   2.345 +        test3   =       1 / acond;
   2.346 +        t1      =   test1 / (1    +  anorm * xnorm / bnorm);
   2.347 +        rtol    =   btol  +  atol *  anorm * xnorm / bnorm;
   2.348 +        """
   2.349 +        %     The following tests guard against extremely small values of
   2.350 +        %     atol, btol  or  ctol.  (The user may have set any or all of
   2.351 +        %     the parameters  atol, btol, conlim  to 0.)
   2.352 +        %     The effect is equivalent to the normal tests using
   2.353 +        %     atol = eps,  btol = eps,  conlim = 1/eps.
   2.354 +        """
   2.355 +        if itn >= itnlim  : istop = 7; 
   2.356 +        if 1 + test3  <= 1: istop = 6; 
   2.357 +        if 1 + test2  <= 1: istop = 5; 
   2.358 +        if 1 + t1     <= 1: istop = 4; 
   2.359 +        """
   2.360 +        %     Allow for tolerances set by the user.
   2.361 +        """
   2.362 +        if  test3 <= ctol:  istop = 3;
   2.363 +        if  test2 <= atol:  istop = 2;
   2.364 +        if  test1 <= rtol:  istop = 1;
   2.365 +        """
   2.366 +        %     See if it is time to print something.
   2.367 +        """
   2.368 +        prnt = False;
   2.369 +        if n     <= 40       : prnt = True;
   2.370 +        if itn   <= 10       : prnt = True;
   2.371 +        if itn   >= itnlim-10: prnt = True;
   2.372 +        # if itn%10 == 0       : prnt = True;
   2.373 +        if test3 <=  2*ctol  : prnt = True;
   2.374 +        if test2 <= 10*atol  : prnt = True;
   2.375 +        if test1 <= 10*rtol  : prnt = True;
   2.376 +        if istop !=  0       : prnt = True;
   2.377 +
   2.378 +        if prnt:
   2.379 +            if show:
   2.380 +                str1 = '%6g %12.5e'%        (itn,   x[0] );
   2.381 +                str2 = ' %10.3e %10.3e'% (r1norm, r2norm );
   2.382 +                str3 = '  %8.1e %8.1e'%  ( test1,  test2 );
   2.383 +                str4 = ' %8.1e %8.1e'%   ( anorm,  acond );
   2.384 +                print str1, str2, str3, str4
   2.385 +
   2.386 +        if istop != 0: break
   2.387 +
   2.388 +    """
   2.389 +    %     End of iteration loop.
   2.390 +    %     Print the stopping condition.
   2.391 +    """
   2.392 +    if show:
   2.393 +        print ' '
   2.394 +        print 'LSQR finished'
   2.395 +        print msg[istop]
   2.396 +        print ' '
   2.397 +        str1 = 'istop =%8g   r1norm =%8.1e'%  ( istop, r1norm );
   2.398 +        str2 = 'anorm =%8.1e   arnorm =%8.1e'%( anorm, arnorm );
   2.399 +        str3 = 'itn   =%8g   r2norm =%8.1e'%  (   itn, r2norm );
   2.400 +        str4 = 'acond =%8.1e   xnorm  =%8.1e'%( acond, xnorm  );
   2.401 +        print str1+ '   ' +str2
   2.402 +        print str3+ '   ' +str4
   2.403 +        print ' '
   2.404 +
   2.405 +    return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var 
   2.406 +            
   2.407 +    """
   2.408 +    %-----------------------------------------------------------------------
   2.409 +    % End of lsqr.m
   2.410 +    %-----------------------------------------------------------------------
   2.411 +    """
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/examples/cohomology/output.h	Thu Jul 09 02:42:47 2009 -0700
     3.3 @@ -0,0 +1,64 @@
     3.4 +#include <iostream>
     3.5 +#include <sstream>
     3.6 +#include <fstream>
     3.7 +
     3.8 +#include <boost/tuple/tuple.hpp>
     3.9 +#include <boost/tuple/tuple_comparison.hpp>
    3.10 +
    3.11 +#include <cassert>
    3.12 +
    3.13 +bool neq(const Smplx& s1, const Smplx& s2)               
    3.14 +{ 
    3.15 +    Smplx::VertexComparison cmp;
    3.16 +    return cmp(s1, s2) || cmp(s2, s1);
    3.17 +}
    3.18 +
    3.19 +unsigned index(const SimplexVector& v, const Smplx& s, const Generator::Comparison& cmp)
    3.20 +{
    3.21 +    SimplexVector::const_iterator it = std::lower_bound(v.begin(), v.end(), s, cmp);
    3.22 +    while (neq(*it, s)) ++it;
    3.23 +    return it - v.begin();
    3.24 +}
    3.25 +
    3.26 +void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Generator::Comparison& cmp)
    3.27 +{
    3.28 +    unsigned i = 0;
    3.29 +    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
    3.30 +    {
    3.31 +        // std::cout << "Simplex: " << *cur << std::endl;
    3.32 +        bool                sign = true;
    3.33 +        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
    3.34 +        {
    3.35 +            // std::cout << "  " << *bcur << std::endl;
    3.36 +            out << (sign ? 1 : -1) << " ";
    3.37 +            out << index(v, *bcur, cmp) << " " << i << "\n";
    3.38 +            sign = !sign;
    3.39 +        }
    3.40 +        ++i;
    3.41 +    }
    3.42 +}
    3.43 +
    3.44 +void output_vertex_indices(std::ostream& out, const SimplexVector& v)
    3.45 +{
    3.46 +    unsigned i = 0;
    3.47 +    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
    3.48 +    {
    3.49 +        if (cur->dimension() == 0)
    3.50 +            out << i << " " << cur->vertices()[0] << std::endl;
    3.51 +        ++i;
    3.52 +    }
    3.53 +}
    3.54 +
    3.55 +void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime, const Generator::Comparison& cmp)
    3.56 +{
    3.57 +    std::ostringstream istr; istr << '-' << i;
    3.58 +    std::string filename = cocycle_prefix + istr.str() + ".ccl";
    3.59 +    std::ofstream out(filename.c_str());
    3.60 +    out << "# Cocycle born at " << c.birth.get<1>() << std::endl;
    3.61 +    for (Persistence::ZColumn::const_iterator zcur = c.zcolumn.begin(); zcur != c.zcolumn.end(); ++zcur)
    3.62 +    {
    3.63 +        const Smplx& s = **(zcur->si);
    3.64 +        out << (zcur->coefficient <= prime/2 ? zcur->coefficient : zcur->coefficient - prime) << " ";
    3.65 +        out << index(v, s, cmp) << "\n";
    3.66 +    }
    3.67 +}
     4.1 --- a/examples/cohomology/rips-pairwise-cohomology.cpp	Thu Jul 09 00:59:32 2009 -0700
     4.2 +++ b/examples/cohomology/rips-pairwise-cohomology.cpp	Thu Jul 09 02:42:47 2009 -0700
     4.3 @@ -28,7 +28,9 @@
     4.4  typedef     Persistence::Death                                      Death;
     4.5  typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
     4.6  
     4.7 -void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime);
     4.8 +#include "output.h"         // for output_*()
     4.9 +
    4.10 +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);
    4.11  
    4.12  int main(int argc, char* argv[])
    4.13  {
    4.14 @@ -41,9 +43,16 @@
    4.15      Dimension               skeleton;
    4.16      DistanceType            max_distance;
    4.17      ZpField::Element        prime;
    4.18 -    std::string             infilename;
    4.19 +    std::string             infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
    4.20  
    4.21 -    program_options(argc, argv, infilename, skeleton, max_distance, prime);
    4.22 +    program_options(argc, argv, infilename, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
    4.23 +    std::ofstream           bdry_out(boundary_name.c_str());
    4.24 +    std::ofstream           vertices_out(vertices_name.c_str());
    4.25 +    std::ofstream           diagram_out(diagram_name.c_str());
    4.26 +    std::cout << "Boundary matrix: " << boundary_name << std::endl;
    4.27 +    std::cout << "Cocycles:        " << cocycle_prefix << "*.ccl" << std::endl;
    4.28 +    std::cout << "Vertices:        " << vertices_name << std::endl;
    4.29 +    std::cout << "Diagram:         " << diagram_name << std::endl;
    4.30  
    4.31      Timer total_timer; total_timer.start();
    4.32      PointContainer          points;
    4.33 @@ -52,15 +61,19 @@
    4.34      PairDistances           distances(points);
    4.35      Generator               rips(distances);
    4.36      Generator::Evaluator    size(distances);
    4.37 +    Generator::Comparison   cmp(distances);
    4.38      SimplexVector           v;
    4.39      Complex                 c;
    4.40      
    4.41      Timer rips_timer; rips_timer.start();
    4.42      rips.generate(skeleton, max_distance, make_push_back_functor(v));
    4.43 -    std::sort(v.begin(), v.end(), Generator::Comparison(distances));
    4.44 +    std::sort(v.begin(), v.end(), cmp);
    4.45      rips_timer.stop();
    4.46      std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
    4.47  
    4.48 +    output_boundary_matrix(bdry_out, v, cmp);
    4.49 +    output_vertex_indices(vertices_out, v);
    4.50 +
    4.51      Timer persistence_timer; persistence_timer.start();
    4.52      ZpField                 zp(prime);
    4.53      Persistence             p(zp);
    4.54 @@ -81,28 +94,24 @@
    4.55          if (d && (size(*cur) - d->get<1>()) > 0)
    4.56          {
    4.57              AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
    4.58 -            std::cout << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
    4.59 +            diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
    4.60          }
    4.61      }
    4.62 -    // output infinte persistence cocycles
    4.63 +    // output infinte persistence pairs 
    4.64      for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
    4.65 -        std::cout << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
    4.66 +        diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
    4.67      persistence_timer.stop();
    4.68  
    4.69  
    4.70      // p.show_cocycles();
    4.71 -    // Output alive cocycles
    4.72 +    // Output alive cocycles of dimension 1
    4.73 +    unsigned i = 0;
    4.74      for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
    4.75      {
    4.76 -        std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
    4.77 -        for (Persistence::ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
    4.78 -        {
    4.79 -            const Smplx& s = **(zcur->si);
    4.80 -            std::cout << zcur->coefficient << " ";
    4.81 -            for (Smplx::VertexContainer::const_iterator vcur = s.vertices().begin(); vcur != s.vertices().end(); ++vcur)
    4.82 -                std::cout << *vcur << " ";
    4.83 -            std::cout << std::endl;
    4.84 -        }
    4.85 +        if (cur->birth.get<0>() != 1) continue;
    4.86 +        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
    4.87 +        // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
    4.88 +        ++i;
    4.89      }
    4.90      total_timer.stop();
    4.91      rips_timer.check("Rips timer");
    4.92 @@ -110,7 +119,7 @@
    4.93      total_timer.check("Total timer");
    4.94  }
    4.95  
    4.96 -void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime)
    4.97 +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)
    4.98  {
    4.99      namespace po = boost::program_options;
   4.100  
   4.101 @@ -122,8 +131,12 @@
   4.102      visible.add_options()
   4.103          ("help,h",                                                                                  "produce help message")
   4.104          ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2),                  "Dimension of the Rips complex we want to compute")
   4.105 -        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(2),              "Prime p for the field F_p")
   4.106 -        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction");
   4.107 +        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(11),             "Prime p for the field F_p")
   4.108 +        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction")
   4.109 +        ("boundary,b",          po::value<std::string>(&boundary_name),                             "Filename where to output the boundary matrix")
   4.110 +        ("cocycle,c",           po::value<std::string>(&cocycle_prefix),                            "Prefix of the filename where to output the 1-dimensional cocycles")
   4.111 +        ("vertices,v",          po::value<std::string>(&vertices_name),                             "Filename where to output the simplex-vertex mapping")
   4.112 +        ("diagram,d",           po::value<std::string>(&diagram_name),                              "Filename where to output the persistence diagram");
   4.113  #if LOGGING
   4.114      std::vector<std::string>    log_channels;
   4.115      visible.add_options()
   4.116 @@ -132,7 +145,6 @@
   4.117  
   4.118      po::positional_options_description pos;
   4.119      pos.add("input-file", 1);
   4.120 -    pos.add("output-file", 2);
   4.121      
   4.122      po::options_description all; all.add(visible).add(hidden);
   4.123  
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/tools/plot-values/plot.py	Thu Jul 09 02:42:47 2009 -0700
     5.3 @@ -0,0 +1,53 @@
     5.4 +#!/usr/bin/env python
     5.5 +
     5.6 +from    pylab           import scatter, show, cm, colorbar, savefig, axis, \
     5.7 +                               figure, xlim, axes, hsv, subplots_adjust as adjust
     5.8 +from    itertools       import izip
     5.9 +from    sys             import argv, exit
    5.10 +import  os.path         as     osp
    5.11 +
    5.12 +
    5.13 +def plot(val_fn, pts_fn, output_fn):
    5.14 +    points = []
    5.15 +    with open(pts_fn) as fp:
    5.16 +        for line in fp.xreadlines():
    5.17 +            points.append(map(float, line.split()))
    5.18 +    
    5.19 +    values = []
    5.20 +    with open(val_fn) as fp:
    5.21 +        for line in fp.xreadlines():
    5.22 +            values.append(float(line.split()[1]))
    5.23 +
    5.24 +    xx = [pt[0] for pt in points]
    5.25 +    yy = [pt[1] for pt in points]
    5.26 +    print "X:", min(xx), max(xx)
    5.27 +    print "Y:", min(yy), max(yy)
    5.28 +
    5.29 +    m = min(values)
    5.30 +    values = [(v-m) % 1. for v in values]
    5.31 +    print "V:", min(values), max(values)
    5.32 +
    5.33 +    aspect = (max(yy) - min(yy))/(max(xx) - min(xx)) + .1
    5.34 +    # aspect = .5
    5.35 +
    5.36 +    # hsv()
    5.37 +    # fig = figure(figsize = (3,3*aspect))
    5.38 +    # ax = fig.add_axes([-.05,-.1,1.1,1.1])
    5.39 +    ax = axes()
    5.40 +    ax.set_axis_off()
    5.41 +    ax.set_aspect('equal', 'box')
    5.42 +    ax.scatter(xx,yy,s=10,c=values)
    5.43 +    # adjust(0,0,1,1,0,0)
    5.44 +    colorbar()
    5.45 +    fig.savefig(output_fn)
    5.46 +
    5.47 +if __name__ == '__main__':
    5.48 +    if len(argv) < 3:
    5.49 +        print "Usage: %s VALUES POINTS" % argv[0]
    5.50 +        exit()
    5.51 +
    5.52 +    val_fn = argv[1]
    5.53 +    pts_fn  = argv[2]
    5.54 +    output_fn, ext = osp.splitext(val_fn)
    5.55 +    output_fn += '.pdf'
    5.56 +    plot(val_fn, pts_fn, output_fn)
     6.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     6.2 +++ b/tools/plot-values/scatter.py	Thu Jul 09 02:42:47 2009 -0700
     6.3 @@ -0,0 +1,35 @@
     6.4 +#!/usr/bin/env python
     6.5 +
     6.6 +from    pylab           import scatter, show, cm, colorbar, axes, savefig
     6.7 +from    itertools       import izip
     6.8 +from    sys             import argv, exit
     6.9 +import  os.path         as     osp
    6.10 +
    6.11 +
    6.12 +def plot(val1_fn, val2_fn):
    6.13 +    values1 = []
    6.14 +    with open(val1_fn) as fp:
    6.15 +        for line in fp.xreadlines():
    6.16 +            values1.append(float(line.split()[1]))
    6.17 +    
    6.18 +    values2 = []
    6.19 +    with open(val2_fn) as fp:
    6.20 +        for line in fp.xreadlines():
    6.21 +            values2.append(float(line.split()[1]))
    6.22 +    
    6.23 +    values1 = [v % 1. for v in values1]
    6.24 +    values2 = [v % 1. for v in values2]
    6.25 +    print min(values1), max(values2), min(values1), min(values2)
    6.26 +
    6.27 +    scatter(values1, values2)
    6.28 +    axes().set_aspect('equal')
    6.29 +    show()
    6.30 +
    6.31 +if __name__ == '__main__':
    6.32 +    if len(argv) < 3:
    6.33 +        print "Usage: %s VALUES1 VALUES2" % argv[0]
    6.34 +        exit()
    6.35 +
    6.36 +    val1_fn = argv[1]
    6.37 +    val2_fn  = argv[2]
    6.38 +    plot(val1_fn, val2_fn)