--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/cohomology.rst Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,87 @@
+.. _cohomology-parametrization:
+
+Parametrizing a point set using circle valued functions
+=======================================================
+
+The procedure described below is explained in detail in [dSVJ09]_.
+
+.. program:: rips-pairwise-cohomology
+
+One can use :sfile:`examples/cohomology/rips-pairwise-cohomology.cpp` to compute
+persistent pairing of the Rips filtration using the persistent cohomology
+algorithm. It takes as input a file containing a point set in Euclidean space
+(one per line) as well as the following command-line flags:
+
+.. cmdoption:: -p, --prime
+
+ The prime to use in the computation (defaults to 11).
+
+.. cmdoption:: -m, --max-distance
+
+ Maximum cutoff parameter up to which to compute the complex.
+
+.. cmdoption:: -s, --skeleton-dimension
+
+ Skeleton to compute; persistent pairs output will be this number minus 1
+ (defaults to 2).
+
+.. cmdoption:: -b, --boundary
+
+ Filename where to output the boundary matrix.
+
+.. cmdoption:: -c, --cocycle
+
+ Prefix of the filenames where to output the 1-dimensional cocycles.
+
+.. cmdoption:: -v, --vertices
+
+ Filename where to output the simplex vertex mapping.
+
+.. cmdoption:: -d, --diagram
+
+ Filename where to output the persistence diagram.
+
+
+For example::
+
+ rips-pairwise-cohomology points.txt -m 1 -b points.bdry -c points -v points.vrt -d points.dgm
+
+Assuming that at the threshold value of 1 (``-m 1`` above) Rips complex contains
+1-dimensional cocycles, they will be output into filenames of the form
+``points-0.ccl``, ``points-1.ccl``, etc.
+
+Subsequently one can use :sfile:`examples/cohomology/cocycle.py` to assign to
+each vertex of the input point set a circle-valued function. It takes the
+boundary matrix, cocycle, and simplex-vertex map as an input (all produced at
+the previous step)::
+
+ cocycle.py points.bdry points-0.ccl points.vrt
+
+The above command outputs a file ``points-0.val`` which contains values assigned
+to the input points (the lines match the lines of the input file
+``points.txt``, but also contains the indices).
+
+
+Plotting
+--------
+
+Two auxilliary tools allow one to visualize the values assigned to the points
+(using Matplotlib_): :sfile:`tools/plot-values/plot.py` and
+:sfile:`tools/plot-values/scatter.py`::
+
+ plot.py points-0.val points.txt scatter.py points-0.val points-1.val
+
+.. _Matplotlib: http://matplotlib.sourceforge.net/
+
+
+Dependency
+----------
+
+The Python `LSQR code`_ (ported from the `Stanford MATLAB implementation`_ to
+Python by `Jeffery Kline`_) included with Dionysus, and used in
+:sfile:`examples/cohomology/cocycle.py`, requires CVXOPT_.
+
+.. _`LSQR code`: http://pages.cs.wisc.edu/~kline/cvxopt/
+.. _CVXOPT: http://abel.ee.ucla.edu/cvxopt/
+.. _`Stanford MATLAB implementation`: http://www.stanford.edu/group/SOL/software/lsqr.html
+.. _`Jeffery Kline`: http://pages.cs.wisc.edu/~kline/
--- a/doc/examples/index.rst Tue Jul 28 10:58:53 2009 -0700
+++ b/doc/examples/index.rst Tue Jul 28 11:42:03 2009 -0700
@@ -32,3 +32,16 @@
:maxdepth: 1
rips
+
+One may use persistent cohomology algorithm to extract persistent cocycles,
+turn them into harmonic cocycles, and use them to parametrize the input point
+set; for details see [dSVJ09]_. The explanation of how to use Dionysus to
+achieve this is available.
+
+.. toctree::
+ :maxdepth: 1
+
+ cohomology
+
+A simple example of computing persistence of a lower-star filtration is in
+:sfile:`examples/lsfiltration.py`.
--- a/doc/index.rst Tue Jul 28 10:58:53 2009 -0700
+++ b/doc/index.rst Tue Jul 28 11:42:03 2009 -0700
@@ -4,6 +4,17 @@
Dionysus is a C++ library for computing persistent homology. It provides
implementations of the following algorithms:
+.. sidebar:: Contents
+
+ .. toctree::
+ :maxdepth: 1
+
+ get-build-install
+ tutorial
+ examples/index
+ python/overview
+ bibliography
+
* Persistent homology computation [ELZ02]_ [ZC05]_
* Vineyards [CEM06]_ |cpp-only|
* Persistent cohomology computation (described in [dSVJ09]_) |cpp-only|
@@ -13,6 +24,7 @@
* :ref:`Alpha shape construction <alpha-shape-example>` in 2D and 3D
* :ref:`Rips complex construction <rips-example>`
* Cech complex construction |cpp-only|
+ * :ref:`Circle-valued parametrization using persistent cohomology <cohomology-parametrization>`
.. todo::
Document more examples.
@@ -27,17 +39,6 @@
simplicity, elegance, and interactivity of Python. Since they mimick the C++
functionality, their documentation may be a helpful resource for the latter.
-Contents:
-
-.. toctree::
- :maxdepth: 1
-
- get-build-install
- tutorial
- examples/index
- python/overview
- bibliography
-
.. include:: substitutions.aux
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/cocycle.py Tue Jul 28 11:42:03 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 Tue Jul 28 11:42:03 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 Tue Jul 28 11:42:03 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 Tue Jul 28 10:58:53 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Tue Jul 28 11:42:03 2009 -0700
@@ -6,6 +6,7 @@
#include <utilities/containers.h> // for BackInsertFunctor
#include <utilities/timer.h>
+#include <utilities/log.h>
#include <string>
@@ -15,61 +16,110 @@
typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
typedef PairDistances::DistanceType DistanceType;
typedef PairDistances::IndexType Vertex;
-
-typedef CohomologyPersistence<DistanceType> Persistence;
-typedef Persistence::SimplexIndex Index;
-typedef Persistence::Death Death;
-
+
typedef Rips<PairDistances> Generator;
typedef Generator::Simplex Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef SimplexVector::const_iterator SV_const_iterator;
-typedef std::map<Smplx, Index,
- Smplx::VertexComparison> Complex;
-typedef std::vector<Smplx> SimplexVector;
+typedef boost::tuple<Dimension, DistanceType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo, SV_const_iterator> Persistence;
+typedef Persistence::SimplexIndex Index;
+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);
+#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[])
{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
Dimension skeleton;
DistanceType max_distance;
- std::string infilename;
+ ZpField::Element prime;
+ std::string infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
- program_options(argc, argv, infilename, skeleton, max_distance);
+ 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;
read_points(infilename, points);
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();
- Persistence p;
+ ZpField zp(prime);
+ Persistence p(zp);
for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
{
std::vector<Index> boundary;
- for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
- bcur != cur->boundary_end(); ++bcur)
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
boundary.push_back(c[*bcur]);
Index idx; Death d;
- boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), size(*cur));
- c[*cur] = idx;
- if (d && (size(*cur) - *d) > 0)
- std::cout << (cur->dimension() - 1) << " " << *d << " " << size(*cur) << std::endl;
+ bool store = cur->dimension() < skeleton;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
+
+ // c[*cur] = idx;
+ if (store)
+ c.insert(std::make_pair(*cur, idx));
+
+ if (d && (size(*cur) - d->get<1>()) > 0)
+ {
+ AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+ diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+ }
}
+ // output infinte persistence pairs
+ for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+ diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
persistence_timer.stop();
+
+
+ // p.show_cocycles();
+ // Output alive cocycles of dimension 1
+ unsigned i = 0;
+ for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ 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");
persistence_timer.check("Persistence timer");
+ total_timer.check("Total timer");
}
-void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+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;
@@ -81,11 +131,20 @@
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")
- ("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()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
po::positional_options_description pos;
pos.add("input-file", 1);
- pos.add("output-file", 2);
po::options_description all; all.add(visible).add(hidden);
@@ -94,6 +153,11 @@
options(all).positional(pos).run(), vm);
po::notify(vm);
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
if (vm.count("help") || !vm.count("input-file"))
{
std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
--- a/examples/cohomology/triangle-cohomology.cpp Tue Jul 28 10:58:53 2009 -0700
+++ b/examples/cohomology/triangle-cohomology.cpp Tue Jul 28 11:42:03 2009 -0700
@@ -54,7 +54,7 @@
stderrLog.subscribeTo(RLOG_CHANNEL("info"));
stderrLog.subscribeTo(RLOG_CHANNEL("error"));
- //stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
+ stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
#endif
SimplexVector v;
@@ -66,10 +66,13 @@
// Compute persistence
Complex c;
- Persistence p;
+ ZpField zp(11);
+ Persistence p(zp);
unsigned i = 0;
for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
{
+ std::cout << "-------" << std::endl;
+
std::vector<Index> boundary;
for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
bcur != cur->boundary_end(); ++bcur)
@@ -84,7 +87,7 @@
// (i.e. when a 1 class kills 0, it's really that in cohomology forward 0 kills 1,
// in cohomology backward 1 kills 0, and in homology 1 kills 0)
- //p.show_cycles();
+ p.show_cocycles();
}
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/consistency/make-zigzag-subsamples.py Tue Jul 28 11:42:03 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 Tue Jul 28 11:42:03 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 Tue Jul 28 10:58:53 2009 -0700
+++ b/examples/rips/rips-pairwise.cpp Tue Jul 28 11:42:03 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 Tue Jul 28 10:58:53 2009 -0700
+++ b/include/geometry/l2distance.h Tue Jul 28 11:42:03 2009 -0700
@@ -2,6 +2,7 @@
#define __L2_DISTANCE_H__
#include <utilities/types.h>
+#include <utilities/log.h>
#include <vector>
#include <fstream>
--- a/include/topology/cohomology-persistence.h Tue Jul 28 10:58:53 2009 -0700
+++ b/include/topology/cohomology-persistence.h Tue Jul 28 11:42:03 2009 -0700
@@ -16,6 +16,7 @@
#include <list>
#include <utility>
+#include <topology/field-arithmetic.h>
#include "utilities/types.h"
#include <boost/optional.hpp>
@@ -23,86 +24,105 @@
namespace bi = boost::intrusive;
-template<class BirthInfo_, class SimplexData_ = Empty<> >
+template<class BirthInfo_, class SimplexData_ = Empty<>, class Field_ = ZpField>
class CohomologyPersistence
{
public:
typedef BirthInfo_ BirthInfo;
typedef SimplexData_ SimplexData;
+ typedef Field_ Field;
+
+ typedef typename Field::Element FieldElement;
- struct SNode;
- typedef bi::list<SNode, bi::constant_time_size<false> > ZRow;
+ CohomologyPersistence(const Field& field = Field()):
+ field_(field) {}
+
- // Simplex representation
- struct SHead: public SimplexData
- {
- SHead(const SHead& other):
- SimplexData(other), order(other.order) {} // don't copy row since we can't
- SHead(const SimplexData& sd, unsigned o):
- SimplexData(sd), order(o) {}
+ // An entry in a cocycle column
+ struct SNode; // members: si, coefficient, ci
+ typedef s::vector<SNode> ZColumn;
+ typedef bi::list<SNode, bi::constant_time_size<false> > ZRow;
+ class CompareSNode;
- // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
- ZRow row;
- unsigned order;
- };
-
+ struct SHead; // members: row, order
typedef s::list<SHead> Simplices;
typedef typename Simplices::iterator SimplexIndex;
- struct Cocycle;
+ struct Cocycle; // members: zcolumn, birth, order
typedef s::list<Cocycle> Cocycles;
typedef typename Cocycles::iterator CocycleIndex;
-
- // An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
- typedef bi::list_base_hook<bi::link_mode<bi::auto_unlink> > auto_unlink_hook;
- struct SNode: public auto_unlink_hook
- {
- SNode() {}
- SNode(SimplexIndex sidx): si(sidx) {}
-
- // eventually store a field element
-
- SimplexIndex si;
- CocycleIndex cocycle; // TODO: is there no way to get rid of this overhead?
-
- void unlink() { auto_unlink_hook::unlink(); }
- };
- class CompareSNode;
+ typedef std::pair<CocycleIndex, FieldElement> CocycleCoefficientPair;
- typedef s::vector<SNode> ZColumn;
- struct Cocycle
- {
- Cocycle(const BirthInfo& b, unsigned o):
- birth(b), order(o) {}
-
- ZColumn cocycle;
- BirthInfo birth;
- unsigned order;
-
- bool operator<(const Cocycle& other) const { return order > other.order; }
- bool operator==(const Cocycle& other) const { return order == other.order; }
- };
-
- typedef boost::optional<BirthInfo> Death;
- typedef std::pair<SimplexIndex,
- Death> IndexDeathPair;
+ typedef boost::optional<BirthInfo> Death;
+ typedef std::pair<SimplexIndex, Death> IndexDeathPair;
// return either a SimplexIndex or a Death
// BI = BoundaryIterator; it should dereference to a SimplexIndex
template<class BI>
- IndexDeathPair add(BI begin, BI end, BirthInfo b, const SimplexData& sd = SimplexData());
+ IndexDeathPair add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData());
- void show_cycles() const;
-
+ void show_cocycles() const;
+ CocycleIndex begin() { return cocycles_.begin(); }
+ CocycleIndex end() { return cocycles_.end(); }
private:
- void add_cocycle(Cocycle& z1, Cocycle& z2);
+ void add_cocycle(CocycleCoefficientPair& z1, CocycleCoefficientPair& z2);
private:
Simplices simplices_;
Cocycles cocycles_;
+ Field field_;
};
+
+// Simplex representation
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SHead: public SimplexData
+{
+ SHead(const SHead& other):
+ SimplexData(other), order(other.order) {} // don't copy row since we can't
+ SHead(const SimplexData& sd, unsigned o):
+ SimplexData(sd), order(o) {}
+
+ // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
+ ZRow row;
+ unsigned order;
+};
+
+// An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
+typedef bi::list_base_hook<bi::link_mode<bi::auto_unlink> > auto_unlink_hook;
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::SNode: public auto_unlink_hook
+{
+ SNode(const SNode& other):
+ si(other.si), coefficient(other.coefficient),
+ ci(other.ci) {}
+
+ SNode(SimplexIndex sidx, FieldElement coef, CocycleIndex cidx):
+ si(sidx), coefficient(coef), ci(cidx) {}
+
+ SimplexIndex si;
+ FieldElement coefficient;
+
+ CocycleIndex ci; // TODO: is there no way to get rid of this overhead?
+
+ void unlink() { auto_unlink_hook::unlink(); }
+};
+
+template<class BirthInfo_, class SimplexData_, class Field_>
+struct CohomologyPersistence<BirthInfo_, SimplexData_, Field_>::Cocycle
+{
+ Cocycle(const BirthInfo& b, unsigned o):
+ birth(b), order(o) {}
+
+ ZColumn zcolumn;
+ BirthInfo birth;
+ unsigned order;
+
+ bool operator<(const Cocycle& other) const { return order > other.order; }
+ bool operator==(const Cocycle& other) const { return order == other.order; }
+};
+
#include "cohomology-persistence.hpp"
--- a/include/topology/cohomology-persistence.hpp Tue Jul 28 10:58:53 2009 -0700
+++ b/include/topology/cohomology-persistence.hpp Tue Jul 28 11:42:03 2009 -0700
@@ -9,52 +9,62 @@
static rlog::RLogChannel* rlCohomology = DEF_CHANNEL("topology/cohomology", rlog::Log_Debug);
#endif
-template<class BirthInfo, class SimplexData>
-class CohomologyPersistence<BirthInfo, SimplexData>::CompareSNode
+template<class BirthInfo, class SimplexData, class Field>
+class CohomologyPersistence<BirthInfo, SimplexData, Field>::CompareSNode
{
public:
bool operator()(const SNode& s1, const SNode& s2) const { return s1.si->order < s2.si->order; }
};
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
template<class BI>
-typename CohomologyPersistence<BirthInfo, SimplexData>::IndexDeathPair
-CohomologyPersistence<BirthInfo, SimplexData>::
-add(BI begin, BI end, BirthInfo birth, const SimplexData& sd)
+typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd)
{
// Create simplex representation
simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
SimplexIndex si = boost::prior(simplices_.end());
// Find out if there are cocycles that evaluate to non-zero on the new simplex
- typedef std::list<CocycleIndex> Candidates;
+ typedef std::list<CocycleCoefficientPair> Candidates;
Candidates candidates, candidates_bulk;
rLog(rlCohomology, "Boundary");
+ bool sign = true; // TODO: this is very crude, fix later
for (BI cur = begin; cur != end; ++cur)
{
- rLog(rlCohomology, " %d", (*cur)->order);
+ rLog(rlCohomology, " %d %d", (*cur)->order, sign);
for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
- candidates_bulk.push_back(zcur->cocycle);
+ candidates_bulk.push_back(std::make_pair(zcur->ci,
+ (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
+ sign = !sign;
}
- candidates_bulk.sort(make_indirect_comparison(std::less<Cocycle>()));
+ candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
+#if LOGGING
rLog(rlCohomology, " Candidates bulk");
for (typename Candidates::iterator cur = candidates_bulk.begin();
cur != candidates_bulk.end(); ++cur)
- rLog(rlCohomology, " %d", (*cur)->order);
+ rLog(rlCohomology, " %d %d", cur->first->order, cur->second);
+#endif
- // Remove duplicates --- this is really Z_2, we need a more sophisticated
+ // Remove duplicates
{
typename Candidates::const_iterator cur = candidates_bulk.begin();
while (cur != candidates_bulk.end())
{
typename Candidates::const_iterator next = cur;
- unsigned count = 0;
- while (next != candidates_bulk.end() && *next == *cur) { ++next; ++count; }
+ FieldElement sum = field_.zero();
+ while (next != candidates_bulk.end() && next->first == cur->first)
+ {
+ sum = field_.add(sum, next->second);
+ ++next;
+ }
- if (count % 2)
- candidates.push_back(*cur);
+ rLog(rlCohomology, " Bulk candidate %d sum %d", cur->first->order, sum);
+ if (!field_.is_zero(sum))
+ candidates.push_back(CocycleCoefficientPair(cur->first, sum));
cur = next;
}
@@ -63,16 +73,22 @@
// Birth
if (candidates.empty())
{
- rLog(rlCohomology, "Birth");
+ if (!store)
+ {
+ simplices_.pop_back();
+ return std::make_pair(simplices_.begin(), Death()); // TODO: shouldn't return front
+ }
unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
cocycles_.push_front(Cocycle(birth, order));
+
+ rLog(rlCohomology, "Birth: %d", cocycles_.front().order);
// set up the cocycle
- ZColumn& cocycle = cocycles_.front().cocycle;
- cocycle.push_back(si);
- cocycle.front().cocycle = cocycles_.begin();
- si->row.push_back(cocycles_.front().cocycle.front());
+ ZColumn& cocycle = cocycles_.front().zcolumn;
+ cocycle.push_back(SNode(si, field_.id(), cocycles_.begin()));
+ si->row.push_back(cocycles_.front().zcolumn.front());
+ rLog(rlCohomology, " Cocyle: %d", si->order);
return std::make_pair(si, Death());
}
@@ -80,67 +96,97 @@
// Death
rLog(rlCohomology, "Death");
-#if 0
+#if LOGGING
// Debug only, output candidates
rLog(rlCohomology, " Candidates");
- for (typename Candidates::iterator cur = candidates.begin();
- cur != candidates.end(); ++cur)
- rLog(rlCohomology, " %d", (*cur)->order);
+ for (typename Candidates::iterator cur = candidates.begin(); cur != candidates.end(); ++cur)
+ rLog(rlCohomology, " %d %d", cur->first->order, cur->second);
#endif
- Cocycle& z = *candidates.front();
- Death d = z.birth;
+ CocycleCoefficientPair& z = candidates.front();
+ Death d = z.first->birth;
// add z to everything else in candidates
for (typename Candidates::iterator cur = boost::next(candidates.begin());
cur != candidates.end(); ++cur)
- add_cocycle(**cur, z);
+ add_cocycle(*cur, z);
- for (typename ZColumn::iterator cur = z.cocycle.begin(); cur != z.cocycle.end(); ++cur)
+ for (typename ZColumn::iterator cur = z.first->zcolumn.begin(); cur != z.first->zcolumn.end(); ++cur)
cur->unlink();
- cocycles_.erase(candidates.front());
+ cocycles_.erase(candidates.front().first);
return std::make_pair(si, d);
}
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
void
-CohomologyPersistence<BirthInfo, SimplexData>::
-show_cycles() const
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+show_cocycles() const
{
- std::cout << "Cocycles" << std::endl;
+ std::cout << "Cocycles: " << cocycles_.size() << std::endl;
for (typename Cocycles::const_iterator cur = cocycles_.begin(); cur != cocycles_.end(); ++cur)
{
- std::cout << cur->order << ": ";
- for (typename ZColumn::const_iterator zcur = cur->cocycle.begin(); zcur != cur->cocycle.end(); ++zcur)
- std::cout << zcur->si->order << ", ";
+ // std::cout << cur->order << " (" << cur->birth << "): ";
+ for (typename ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
+ std::cout << zcur->coefficient << " * " << zcur->si->order << ", ";
std::cout << std::endl;
}
}
-template<class BirthInfo, class SimplexData>
+template<class BirthInfo, class SimplexData, class Field>
void
-CohomologyPersistence<BirthInfo, SimplexData>::
-add_cocycle(Cocycle& to, Cocycle& from)
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add_cocycle(CocycleCoefficientPair& to, CocycleCoefficientPair& from)
{
- rLog(rlCohomology, "Adding cocycle %d to %d", from.order, to.order);
+ rLog(rlCohomology, "Adding cocycle %d to %d", from.first->order, to.first->order);
+
+ ZColumn nw;
+ FieldElement multiplier = field_.neg(field_.div(to.second, from.second));
+ CocycleIndex ci = to.first;
- ZColumn nw;
+ typename ZColumn::iterator tcur = to.first->zcolumn.begin();
+ typename ZColumn::iterator fcur = from.first->zcolumn.begin();
+ CompareSNode cmp;
+ while (tcur != to.first->zcolumn.end() && fcur != from.first->zcolumn.end())
+ {
+ rLog(rlCohomology, " %d %d", tcur->si->order, fcur->si->order);
+ if (cmp(*tcur, *fcur))
+ {
+ nw.push_back(*tcur);
+ ++tcur;
+ }
+ else if (cmp(*fcur, *tcur))
+ {
+ nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+ ++fcur;
+ }
+ else // equality
+ {
+ FieldElement res = field_.mul(multiplier, tcur->coefficient);
+ res = field_.add(fcur->coefficient, res);
+ if (!field_.is_zero(res))
+ nw.push_back(SNode(fcur->si, res, ci));
+ ++tcur; ++fcur;
+ }
+ }
+ for (; tcur != to.first->zcolumn.end(); ++tcur)
+ {
+ rLog(rlCohomology, " %d", tcur->si->order);
+ nw.push_back(SNode(*tcur));
+ }
+ for (; fcur != from.first->zcolumn.end(); ++fcur)
+ {
+ rLog(rlCohomology, " %d", fcur->si->order);
+ nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
+ }
- std::set_symmetric_difference(to.cocycle.begin(), to.cocycle.end(), // this is catered to Z_2
- from.cocycle.begin(), from.cocycle.end(),
- std::back_insert_iterator<ZColumn>(nw),
- CompareSNode());
- for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
+ for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
cur->unlink();
- to.cocycle.swap(nw);
+ to.first->zcolumn.swap(nw);
- for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
- {
+ for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
cur->si->row.push_back(*cur);
- cur->cocycle = nw[0].cocycle;
- }
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/field-arithmetic.h Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,64 @@
+#ifndef __FIELD_ARITHMETIC_H__
+#define __FIELD_ARITHMETIC_H__
+
+#include <vector>
+#include <gmpxx.h>
+
+class ZpField
+{
+ public:
+ typedef int Element;
+
+ ZpField(Element p = 2);
+
+ Element id() const { return 1; }
+ Element zero() const { return 0; }
+
+ Element neg(Element a) const { return p_ - a; }
+ Element add(Element a, Element b) const { return (a+b) % p_; }
+
+ Element inv(Element a) const { return inverses_[a]; }
+ Element mul(Element a, Element b) const { return (a*b) % p_; }
+ Element div(Element a, Element b) const { return mul(a, inv(b)); }
+
+ bool is_zero(Element a) const { return (a % p_) == 0; }
+
+ private:
+ Element p_;
+ std::vector<Element> inverses_;
+};
+
+ZpField::
+ZpField(Element p):
+ p_(p), inverses_(p_)
+{
+ for (Element i = 1; i < p_; ++i)
+ for (Element j = 1; j < p_; ++j)
+ if (mul(i,j) == 1)
+ {
+ inverses_[i] = j;
+ break;
+ }
+}
+
+class QField
+{
+ public:
+ typedef mpq_class Element;
+
+ QField() {}
+
+ Element id() const { return 1; }
+ Element zero() const { return 0; }
+
+ Element neg(Element a) const { return -a; }
+ Element add(Element a, Element b) const { return (a+b); }
+
+ Element inv(Element a) const { return id()/a; }
+ Element mul(Element a, Element b) const { return (a*b); }
+ Element div(Element a, Element b) const { return a/b; }
+
+ bool is_zero(Element a) const { return a == 0; }
+};
+
+#endif // __FIELD_ARITHMETIC_H__
--- a/include/utilities/indirect.h Tue Jul 28 10:58:53 2009 -0700
+++ b/include/utilities/indirect.h Tue Jul 28 11:42:03 2009 -0700
@@ -28,6 +28,27 @@
{ return IndirectComparison<Comparison>(cmp); }
+template<class Comparison_>
+class FirstComparison
+{
+ public:
+ typedef Comparison_ Comparison;
+
+ FirstComparison(Comparison cmp):
+ cmp_(cmp) {}
+
+ template<class Pair>
+ bool operator()(const Pair& a, const Pair& b) const
+ { return cmp_(a.first, b.first); }
+
+ private:
+ Comparison cmp_;
+};
+template<class Comparison>
+FirstComparison<Comparison> make_first_comparison(const Comparison& cmp)
+{ return FirstComparison<Comparison>(cmp); }
+
+
template<class Comparison>
struct ThreeOutcomeCompare: public Comparison
{
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/plot-values/plot.py Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,52 @@
+#!/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)
+
+ # hsv()
+ fig = figure()
+ scatter(xx,yy,s=10,c=values)
+ colorbar()
+
+ # ax = fig.add_axes([-.05,-.1,1.1,1.1])
+ ax = axes()
+ ax.set_axis_off()
+ ax.set_aspect('equal', 'box')
+ # adjust(0,0,1,1,0,0)
+
+ 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 Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,43 @@
+#!/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, outfn = None):
+ 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, s=10)
+ axes().set_aspect('equal')
+ if not outfn:
+ show()
+ else:
+ savefig(outfn)
+
+if __name__ == '__main__':
+ if len(argv) < 3:
+ print "Usage: %s VALUES1 VALUES2 [OUTPUT]" % argv[0]
+ exit()
+
+ val1_fn = argv[1]
+ val2_fn = argv[2]
+
+ outfn = None
+ if len(argv) > 3:
+ outfn = argv[3]
+
+ plot(val1_fn, val2_fn, outfn)