--- a/.hgignore Thu May 07 14:53:13 2009 -0700
+++ b/.hgignore Mon May 11 12:45:49 2009 -0700
@@ -3,3 +3,4 @@
docs
aux
build
+doc/.build
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/Makefile Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,75 @@
+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS =
+SPHINXBUILD = sphinx-build
+PAPER =
+
+# Internal variables.
+PAPEROPT_a4 = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS = -d .build/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+
+.PHONY: help clean html web pickle htmlhelp latex changes linkcheck
+
+help:
+ @echo "Please use \`make <target>' where <target> is one of"
+ @echo " html to make standalone HTML files"
+ @echo " pickle to make pickle files"
+ @echo " json to make JSON files"
+ @echo " htmlhelp to make HTML files and a HTML help project"
+ @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+ @echo " changes to make an overview over all changed/added/deprecated items"
+ @echo " linkcheck to check all external links for integrity"
+
+clean:
+ -rm -rf .build/*
+
+html:
+ mkdir -p .build/html .build/doctrees
+ $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) .build/html
+ @echo
+ @echo "Build finished. The HTML pages are in .build/html."
+
+pickle:
+ mkdir -p .build/pickle .build/doctrees
+ $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) .build/pickle
+ @echo
+ @echo "Build finished; now you can process the pickle files."
+
+web: pickle
+
+json:
+ mkdir -p .build/json .build/doctrees
+ $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) .build/json
+ @echo
+ @echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+ mkdir -p .build/htmlhelp .build/doctrees
+ $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) .build/htmlhelp
+ @echo
+ @echo "Build finished; now you can run HTML Help Workshop with the" \
+ ".hhp project file in .build/htmlhelp."
+
+latex:
+ mkdir -p .build/latex .build/doctrees
+ $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) .build/latex
+ @echo
+ @echo "Build finished; the LaTeX files are in .build/latex."
+ @echo "Run \`make all-pdf' or \`make all-ps' in that directory to" \
+ "run these through (pdf)latex."
+
+changes:
+ mkdir -p .build/changes .build/doctrees
+ $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) .build/changes
+ @echo
+ @echo "The overview file is in .build/changes."
+
+linkcheck:
+ mkdir -p .build/linkcheck .build/doctrees
+ $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) .build/linkcheck
+ @echo
+ @echo "Link check complete; look for any errors in the above output " \
+ "or in .build/linkcheck/output.txt."
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/bibliography.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,34 @@
+Bibliography
+============
+
+.. |edels| replace:: Herbert Edelsbrunner
+.. .. _edels: http://www.cs.duke.edu/~edels/
+.. |david| replace:: David Cohen-Steiner
+.. .. _david: http://www-sop.inria.fr/geometrica/team/David.Cohen-Steiner/
+.. |morozov| replace:: Dmitriy Morozov
+.. .. _morozov: http://www.mrzv.org/
+.. |afra| replace:: Afra Zomorodian
+.. .. _afra: http://www.cs.dartmouth.edu/~afra/
+.. |letscher| replace:: David Letscher
+.. |gunnar| replace:: Gunnar Carlsson
+.. .. _gunnar: http://math.stanford.edu/~gunnar/
+.. |vin| replace:: Vin de Silva
+.. .. _vin: http://pages.pomona.edu/~vds04747/
+.. |mikael| replace:: Mikael Vejdemo Johansson
+
+.. [CdSM09] |gunnar|, |vin|, and |morozov|.
+ `Zigzag Persistent Homology and Real-valued Functions
+ <http://www.mrzv.org/publications/zigzags/>`__.
+
+.. [CEM06] |david|, |edels|, and |morozov|.
+ `Vines and Vineyards by Updating Persistence in Linear Time
+ <http://www.mrzv.org/publications/vineyards/>`__.
+
+.. [dSVJ09] |vin|, |mikael|.
+ Persistent Cohomology and Circular Coordinates.
+
+.. [ELZ02] |edels|, |letscher|, and |afra|.
+ Topological Persistence and Simplification.
+
+.. [ZC05] |afra| and |gunnar|.
+ Computing Persistent Homology.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/conf.py Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,205 @@
+# -*- coding: utf-8 -*-
+#
+# Dionysus documentation build configuration file, created by
+# sphinx-quickstart on Tue Dec 9 10:21:40 2008.
+#
+# This file is execfile()d with the current directory set to its containing dir.
+#
+# The contents of this file are pickled, so don't put values in the namespace
+# that aren't pickleable (module imports are okay, they're removed automatically).
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+import sys, os
+
+# If your extensions are in another directory, add it here. If the directory
+# is relative to the documentation root, use os.path.abspath to make it
+# absolute, like shown here.
+#sys.path.append(os.path.abspath('.'))
+
+# General configuration
+# ---------------------
+
+# Add any Sphinx extension module names here, as strings. They can be extensions
+# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+sys.path.append('.')
+
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.todo',
+ #'sphinx.ext.jsmath', 'ext.sprite_jsmath',
+ 'sphinx.ext.pngmath',
+ 'ext.sfile', 'ext.cppdocs']
+todo_include_todos = True
+
+jsmath_path = 'jsMath/easy/load.js'
+
+# Base URI for sfile extension
+sfile_base_uri = 'http://hg.mrzv.org/Dionysus/file/tip/'
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['.templates']
+
+# The suffix of source filenames.
+source_suffix = '.rst'
+
+# The encoding of source files.
+#source_encoding = 'utf-8'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'Dionysus'
+copyright = u'2005--2009, Dmitriy Morozov'
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = '3'
+# The full version, including alpha/beta/rc tags.
+release = '3'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+#today = ''
+# Else, today_fmt is used as the format for a strftime call.
+#today_fmt = '%B %d, %Y'
+
+# List of documents that shouldn't be included in the build.
+#unused_docs = []
+
+# List of directories, relative to source directory, that shouldn't be searched
+# for source files.
+exclude_trees = ['.build']
+
+# The reST default role (used for this markup: `text`) to use for all documents.
+#default_role = None
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+#add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+#add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+#show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+
+# Options for HTML output
+# -----------------------
+
+# The theme to use for HTML and HTML Help pages. Major themes that come with
+# Sphinx are currently 'default' and 'sphinxdoc'.
+html_theme = 'default'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+#html_theme_path = []
+
+# The name for this set of Sphinx documents. If None, it defaults to
+# "<project> v<release> documentation".
+html_title = 'Dionysus'
+
+# A shorter title for the navigation bar. Default is the same as html_title.
+#html_short_title = None
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+#html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+#html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['.static']
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+html_last_updated_fmt = '%b %d, %Y'
+
+# If true, SmartyPants will be used to convert quotes and dashes to
+# typographically correct entities.
+html_use_smartypants = True
+
+# Custom sidebar templates, maps document names to template names.
+#html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+#html_additional_pages = {}
+
+# If false, no module index is generated.
+#html_use_modindex = True
+
+# If false, no index is generated.
+#html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+#html_split_index = False
+
+# If true, the reST sources are included in the HTML build as _sources/<name>.
+#html_copy_source = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it. The value of this option must be the
+# base URL from which the finished HTML is served.
+#html_use_opensearch = ''
+
+# If nonempty, this is the file name suffix for HTML files (e.g. ".xhtml").
+#html_file_suffix = ''
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'Dionysusdoc'
+
+
+# Options for LaTeX output
+# ------------------------
+
+# The paper size ('letter' or 'a4').
+#latex_paper_size = 'letter'
+
+# The font size ('10pt', '11pt' or '12pt').
+#latex_font_size = '10pt'
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title, author, document class [howto/manual]).
+latex_documents = [
+ ('index', 'Dionysus.tex', ur'Dionysus Documentation',
+ ur'Dmitriy Morozov', 'manual'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+#latex_logo = None
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+#latex_use_parts = False
+
+# Additional stuff for the LaTeX preamble.
+#latex_preamble = ''
+
+# Documents to append as an appendix to all manuals.
+#latex_appendices = []
+
+# If false, no module index is generated.
+#latex_use_modindex = True
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/alphashape.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,36 @@
+.. _alpha-shape-example:
+
+Alpha shape example
+===================
+
+The example given in :sfile:`examples/alphashapes/alphashapes.py` takes a
+filename containing points in 2D or 3D on the command line. It generates the
+alpha shape filtration of those points, and computes its persistence. It then
+outputs the persistence diagram in the format of a point (dimension, birth,
+death) per line.
+
+.. literalinclude:: ../../examples/alphashapes/alphashapes.py
+ :language: python
+
+After the points are read into the list ``points``, the functions
+:ref:`fill_alpha*_complex <alphashapes>` fill the list simplices with the
+simplices of the Delaunay triangulation. Each one has its :attr:`~Simplex.data`
+attribute set to its alpha shape value (the minimum value of the squared
+distance function on its dual Voronoi cell).
+
+The simplices are put into lexicographic order (required for
+:class:`Filtration`), and then a filtration is created that sorts simplices with
+respect to their data and dimension (via :func:`data_dim_cmp`)::
+
+ simplices.sort(vertex_cmp)
+ f = Filtration(simplices, data_dim_cmp)
+
+We initialize :class:`StaticPersistence`, and pair the simplices::
+
+ p = StaticPersistence(f)
+ p.pair_simplices()
+
+Iterating over the :class:`StaticPersistence`, we output the points of the
+persistence diagram (dimension, birth, death) in the last for loop. The ``i ==
+i.pair`` condition indicates that the positive simplex is unpaired (i.e. the
+class it creates survives till infinity).
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/index.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,34 @@
+.. _examples:
+
+Examples
+========
+
+The most basic example and therefore a good place to start getting acquainted
+with the library is the :ref:`triangle-example`. It adds simplices of a triangle
+one by one, and then (in case of a :ref:`triangle-zigzag-example`), removes them
+one by one.
+
+.. toctree::
+
+ triangle
+ triangle-zigzag
+
+The simplest example that instead of specifying the complex explicitly,
+constructs it from the input point set is the :ref:`alpha-shape-example`. The
+example reads points from a file, determines their dimension dynamically (based
+on the number of coordinates in the first line of the file), and then constructs
+an alpha shape and outputs its persistence diagram.
+
+.. toctree::
+
+ alphashape
+
+Another example that follows a similar strategy is the computation of the
+Vietoris-Rips complex. Since only pairwise distances are required it works with
+points in arbitrary dimension. (Of course, in dimensions 2 and 3 the complexes
+are much larger than those for the :ref:`alpha-shape-example`).
+
+.. toctree::
+ :maxdepth: 1
+
+ rips
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/rips.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,71 @@
+.. _rips-example:
+
+Rips complex example
+====================
+
+.. todo::
+ Explain `Vietoris-Rips complex`_.
+
+There is an elementary example in :sfile:`examples/rips/rips.py` that computes a
+Rips complex of a few points with integer coordinates on a line. It illustrates
+the use of Rips complexes and in particular of defining your own notion of a
+:ref:`Distance <Distances>` based on which the Rips complex is constructed.
+
+A more useful example is given in :sfile:`examples/rips/rips-pairwise.py` (and
+its C++ counterpart in :sfile:`examples/rips/rips-pairwise.cpp`). The example
+takes on the command line the filename of a file with points in Euclidean space
+(one point per line), and a cut off parameters for the skeleton and the
+:math:`\epsilon` parameter for the Rips complex construction. It then constructs
+the Rips complex up to these cutoff parameters, computes its persistence, and
+outputs the persistence diagram (one point per line).
+
+.. literalinclude:: ../../examples/rips/rips-pairwise.py
+
+The bit that sets up the Rips complex is::
+
+ distances = PairwiseDistances(points)
+ rips = Rips(distances)
+ simplices = []
+ rips.generate(skeleton, max, simplices.append)
+
+The computation of persistence and output of the persistence diagram are the
+same as in the :ref:`alpha-shape-example`. The example also incorporates
+the :ref:`speed-up-suggestions` given in the :ref:`tutorial`.
+
+
+
+C++ sketch
+----------
+.. highlight:: cpp
+
+.. warning:: This section is not finished.
+
+The example given in :sfile:`examples/rips/rips.cpp` illustrates how one can use
+the library to compute persistence of a `Vietoris-Rips complex`_ for a given set of
+distances. At the top of the file a `struct Distances` is defined. The
+particular distances in the example are trivial (they are points with integer
+coordinates on a real line), however, the `struct` illustrates the basic
+requirements of any such class to be passed to the `Rips<Distances>` class.
+
+.. _`Vietoris-Rips complex`: http://en.wikipedia.org/wiki/Vietoris-Rips_complex
+
+The Rips complex itself is generated in the line::
+
+ rips.generate(2, 50, make_push_back_functor(complex));
+
+which tells it to generate a 2-skeleton of the Rips complex up to
+distance value of 50, and insert the simplices into the previously defined
+vector `complex`.
+
+Subsequent sort is unnecessary since Bron-Kerbosch algorithm that generates the
+complex will actually generate the simplices in lexicographic order; it's there
+for illustration purposes only (the simplices must be sorted
+lexicographically).
+
+The following "paragraph" sets up the filtration with respect to simplex sizes
+(specified by `Generator::Comparison(distances)`), and computes its persistence::
+
+ // Generate filtration with respect to distance and compute its persistence
+ Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
+ Persistence p(f);
+ p.pair_simplices();
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/triangle-zigzag.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,38 @@
+.. _triangle-zigzag-example:
+
+Triangle zigzag example
+=======================
+
+Simple example of a filtered triangle where simplices are first inserted in a
+given order, and then removed in the reverse order is in
+:sfile:`examples/triangle/triangle-zigzag.cpp`. Its Python equivalent
+(:sfile:`examples/triangle/triangle-zigzag.py`) is described next.
+
+.. literalinclude:: ../../examples/triangle/triangle-zigzag.py
+ :language: python
+
+Unlike the :ref:`triangle-example`, here we use :class:`ZigzagPersistence` to
+compute the pairings, and therefore need to store the internal representations
+of the simplicies used by the class. These representation are stored in the
+dictionary ``complex``, which maps the simplices to their representations for
+:class:`ZigzagPersistence`.
+
+The first for loop processes the simplices sorted with respect to
+:func:`data_cmp`. :meth:`ZigzagPersistence.add` invoked within the loop accepts
+the boundary of the newly added cell in its internal representation, which is
+computed by looking up each simplex in the dictionary ``complex``:
+``[complex[ss] for ss in s.boundary]``. If there is a birth, the value to be
+associated with the newly created class is ``b`` (which in this case is simply a
+counter). :meth:`~ZigzagPersistence.add` returns a pair ``(i,d)``. The former
+is an internal representation of the newly added cell, which we immediately
+record with ``complex[s] = i``. The latter is an indicator of whether a death
+has occurred, which happens iff ``d is not None``, in which case ``d`` is the
+birth value passed to :meth:`~ZigzagPersistence.add` whenever the class that
+just died was born. If the death occurred, then we outut the interval ``(d,
+b-1)``.
+
+The second for loop removes simplices in the reverse order of their insertion.
+:meth:`~ZigzagPersistence.remove` takes the index of the cells to be removed
+(looked up in the ``complex`` dictionary: ``complex[s]``), and takes a birth
+value in case a class is born. It return only a death indicator (which again is
+``None`` if no death occurred).
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/examples/triangle.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,43 @@
+.. _triangle-example:
+
+Triangle example
+================
+
+Simple example of a filtered triangle is given in
+:sfile:`examples/triangle/triangle.cpp`. Its equivalent in Python appears in
+:sfile:`examples/triangle/triangle.py`, and we describe it next.
+
+.. literalinclude:: ../../examples/triangle/triangle.py
+ :language: python
+
+After the necessary imports, the ``complex`` is setup explicitly as a list of
+simplices. Each :class:`Simplex` constructor takes an iterable sequence of
+vertices, and optionally a data value.
+
+The complex must be sorted lexicographically to compute persistence using
+:class:`StaticPersistence`, and it is accomplished via a helper comparison function :func:`vertex_cmp`, which compares simplices with respect to the lexicographic ordering::
+
+ complex.sort(vertex_cmp)
+
+A filtration ``f`` is initialized using the :class:`Filtration` class, which
+takes a list of lexicographically sorted simplices and a comparison that defines
+in what order the simplices should come in the filtration. In this case we use
+:func:`data_cmp`, which simply compares simplices' :attr:`~Simplex.data`
+attributes.
+
+:class:`StaticPersistence` is initialized with the filtration, and its method
+:meth:`~StaticPersistence.pair_simplices` pairs the simplices of the
+filtration::
+
+ p = StaticPersistence(f)
+ p.pair_simplices()
+
+Subsequently, we iterate over ``p`` to access a representation of each simplex
+in the filtration order. We output each simplex, its sign, and its pair. Note
+``complex[f[p(i)]]``: ``p(i)`` gives the integer index of the iterator ``i`` in
+the filtration ``f``; in turn ``f[p(i)]`` gives the index of the simplex in the
+(lexicographically ordered) ``complex``. So the entire expression returns
+individual simplices. Naturally, one could use this to access the
+:attr:`~Simplex.data` attribute of the simplices to output the actual
+persistence diagram, as is done in the :ref:`alpha-shape-example` and the
+:ref:`rips-example`.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/ext/cppdocs.py Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,202 @@
+from docutils import nodes, utils
+from docutils.parsers.rst import directives
+from sphinx import addnodes
+from sphinx.util.compat import Directive
+
+from sphinx.directives.desc import CDesc
+
+############
+# CppClass #
+############
+class CppClass(nodes.Part, nodes.Element):
+ @staticmethod
+ def html_visit(visitor, node):
+ visitor.body.append('<dl class="class">')
+ visitor.body.append(visitor.starttag(node, 'dt'))
+ if node['tparams']:
+ visitor.body.append('<em class="property">template</em> <')
+ tparam_names = ['class ' + tp[0] for tp in node['tparams']]
+ visitor.body.append(', '.join(tparam_names))
+ visitor.body.append('><br/>')
+
+ visitor.body.append('<em class="property">class</em> ')
+ visitor.body.append('<tt class="descname">' + node['name'] + '</tt>')
+ visitor.body.append(u'<a class="headerlink" href="#%s" title="%s">\u00B6</a>' %
+ ('cppclass-'+node['name'], _('Permalink to this class')))
+ visitor.body.append('</dt>')
+ visitor.body.append('<dd>')
+
+ if node['derives']:
+ visitor.body.append('<p>Derives from ')
+ for d in node['derives'].split(','):
+ dnode ={ 'type': 'cppclass', 'target': d.strip() }
+ Ref.html_visit(visitor, dnode)
+ visitor.body.append(d)
+ Ref.html_depart(visitor, dnode)
+ visitor.body.append('.</p>')
+
+ visitor.body.append('<p>')
+ for name, desc in node['tparams']:
+ visitor.body.append(name + ' — ' + desc + '<br/>')
+ visitor.body.append('</p>')
+
+
+ @staticmethod
+ def html_depart(visitor, node):
+ visitor.body.append('</dd></dl>')
+
+
+class CppClassDirective(Directive):
+ has_content = True
+ required_arguments = 1
+ final_argument_whitespace = True
+ option_spec = { 'derives': directives.unchanged,
+ 'tparam': directives.unchanged }
+
+ def run(self):
+ env = self.state.document.settings.env
+
+ cppclass = CppClass()
+ cppclass['name'] = self.arguments[0]
+ cppclass['derives'] = self.options.get('derives')
+
+ targetname = '%s-%s' % ('cppclass', cppclass['name'])
+ targetnode = nodes.target('', '', ids=[targetname])
+ self.state.document.note_explicit_target(targetnode)
+
+ indextext = _('%s (C++ class)') % cppclass['name']
+ inode = addnodes.index(entries = [('single', indextext,
+ 'cppclass-' + cppclass['name'],
+ cppclass['name'])])
+
+ self.state.nested_parse(self.content, self.content_offset, cppclass)
+
+ return [inode, targetnode, cppclass]
+
+class TParam(nodes.Element):
+ pass
+
+class TParamDirective(Directive):
+ required_arguments = 1
+ optional_arguments = 1
+ final_argument_whitespace = True
+
+ def run(self):
+ tparam = TParam()
+
+ tparam['name'] = self.arguments[0]
+ if len(self.arguments) > 1:
+ tparam['description'] = self.arguments[1]
+
+ return [tparam]
+
+
+#############
+# CppMethod #
+#############
+class CppMethod(nodes.Part, nodes.Element):
+ @staticmethod
+ def html_visit(visitor, node):
+ visitor.body.append(visitor.starttag(node, 'dt'))
+ visitor.body.append(node['name'])
+ visitor.body.append(u'<a class="headerlink" href="#%s" title="%s">\u00B6</a>' %
+ ('cppmethod-' + node['classname'] + '::' + node['name'], _('Permalink to this class')))
+ visitor.body.append('</dt>')
+ visitor.body.append('<dd>')
+
+ @staticmethod
+ def html_depart(visitor, node):
+ visitor.body.append('</dd></dl>')
+
+
+class CppMethodDirective(Directive):
+ has_content = True
+ required_arguments = 1
+ final_argument_whitespace = True
+
+ def run(self):
+ env = self.state.document.settings.env
+
+ cppmethod = CppMethod()
+ cppmethod['name'] = self.arguments[0] # TODO: parse name
+
+ targetname = '%s-%s' % ('cppmethod', cppmethod['name'])
+ targetnode = nodes.target('', '', ids=[targetname])
+ self.state.document.note_explicit_target(targetnode)
+
+ indextext = _('%s (C++ method)') % cppmethod['name']
+ inode = addnodes.index(entries = [('single', indextext,
+ 'cppmethod-' + cppmethod['name'],
+ cppmethod['name'])])
+
+ self.state.nested_parse(self.content, self.content_offset, cppmethod)
+
+ return [inode, targetnode, cppmethod]
+
+
+class Ref(nodes.Inline, nodes.TextElement):
+ @staticmethod
+ def html_visit(visitor, node):
+ if node['type'] == 'cppclass':
+ visitor.body.append('<a href="#%s-%s">' % (node['type'], node['target']))
+ elif node['type'] == 'cppmethod':
+ # TODO: check if the name is not fully qualified, and has a parent CppClass node,
+ # in that case, prepend the name, otherwise
+ visitor.body.append('<a href="#%s-%s::%s">' % (node['type'], node['classname'], node['target']))
+
+
+ @staticmethod
+ def html_depart(visitor, node):
+ visitor.body.append('</a>')
+
+def cppclass_role(role, rawtext ,text, lineno, inliner, options={}, content=[]):
+ text = utils.unescape(text)
+ node = Ref(text, text, target=text, type='cppclass')
+ node['docname'] = inliner.document.settings.env.docname
+ return [node], []
+
+def cppmethod_role(role, rawtext ,text, lineno, inliner, options={}, content=[]):
+ text = utils.unescape(text)
+ node = Ref(text, text, target=text, type='cppmethod')
+ node['docname'] = inliner.document.settings.env.docname
+ return [node], []
+
+
+def process_classnames(app, doctree, fromdocname):
+ for node in doctree.traverse(CppClass):
+ for method in node.traverse(CppMethod):
+ method['classname'] = node['name']
+ for ref in node.traverse(Ref):
+ if ref['type'] == 'cppmethod':
+ ref['classname'] = node['name']
+
+def process_tparams(app, doctree, fromdocname):
+ for node in doctree.traverse(CppClass):
+ node['tparams'] = []
+ for tparam in node.traverse(TParam):
+ node['tparams'].append((tparam['name'], tparam['description']))
+
+
+def process_cfunction_scope(app, doctree, fromdocname):
+ for node in doctree.traverse():
+ if 'ctype' in node: print node
+ if 'cfunction' in node: print node
+
+def setup(app):
+ app.add_node(Ref, html=(Ref.html_visit, Ref.html_depart))
+
+ app.add_directive('cppclass', CppClassDirective)
+ app.add_node(CppClass, html=(CppClass.html_visit, CppClass.html_depart))
+ app.add_role('cppclass', cppclass_role)
+
+ app.add_directive('tparam', TParamDirective)
+ app.add_node(TParam, html=(lambda v,n: '', lambda v,n: ''))
+
+ app.add_directive('cppmethod', CppMethodDirective)
+ app.add_node(CppMethod, html=(CppMethod.html_visit, CppMethod.html_depart))
+ app.add_role('cppmethod', cppmethod_role)
+
+ app.connect('doctree-resolved', process_classnames)
+ app.connect('doctree-resolved', process_tparams)
+
+ app.connect('doctree-resolved', process_cfunction_scope)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/ext/sfile.py Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,14 @@
+from docutils import nodes, utils
+import posixpath
+
+def sfile_role(typ, rawtext, etext, lineno, inliner, options={}, content=[]):
+ env = inliner.document.settings.env
+ baseuri = env.config.sfile_base_uri
+ text = utils.unescape(etext)
+ refnode = nodes.reference('', '', refuri=posixpath.join(baseuri, text))
+ refnode += nodes.literal(text, text)
+ return [refnode], []
+
+def setup(app):
+ app.add_role('sfile', sfile_role)
+ app.add_config_value('sfile_base_uri', 'http://example.com/source', True)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/ext/sprite_jsmath.py Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,2 @@
+def setup(app):
+ app.add_javascript('jsMath/plugins/spriteImageFonts.js')
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/get-build-install.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,123 @@
+Get, Build, Install
+===================
+
+The most up to date code is available from
+`my Mercurial repository`_.
+If you have Mercurial_, the easiest way to obtain the code is by cloning it:
+
+.. parsed-literal::
+
+ hg clone |dionysus-url|
+ cd Dionysus
+ hg up tip
+
+If you don't have time or desire to deal with Mercurial, you can download the
+`tarball of the entire repository`_. The advantage of using Mercurial is that it
+makes it very easy to keep up with the updates that are periodically committed
+to the repository::
+
+ hg pull -u
+
+
+.. |dionysus-url| replace:: http://hg.mrzv.org/Dionysus/
+
+.. _Mercurial: http://www.selenic.com/mercurial/
+
+.. _`tarball of the entire repository`: http://hg.mrzv.org/Dionysus/archive/tip.tar.gz
+.. _`my Mercurial repository`: http://hg.mrzv.org/Dionysus/
+
+
+Dependencies
+------------
+Dionysus requires the following software:
+
+ :CMake_: for building (version :math:`\geq` 2.6)
+ :Boost_: C++ utilities (version :math:`\geq` 1.36; including Boost.Python used to create
+ Python bindings)
+
+There also seems to be a dependence on the version of GCC, although I don't
+entirely understand it. GCC :math:`\geq` 4.3 definitely works, but some versions
+below that seem not to. There were some major changes in 4.3, so it's not
+entirely surprising, but I don't understand all the subtleties.
+
+Optional dependencies:
+
+ :CGAL_: for alpha shapes (version :math:`\geq` 3.4)
+ :rlog_: used for logging only (not needed by default)
+
+.. :dsrpdb_: for reading PDB files
+ :SYNAPS_: for solving polynomials (for kinetic kernel), which in
+ turn requires GMP_
+
+.. _CMake: http://www.cmake.org
+.. _Boost: http://www.boost.org
+.. _CGAL: http://www.cgal.org
+.. _rlog: http://www.arg0.net/rlog
+.. _dsrpdb: http://www.salilab.org/~drussel/pdb/
+.. _SYNAPS: http://www-sop.inria.fr/galaad/synaps/
+.. _GMP: http://gmplib.org/
+
+
+Building
+--------
+To build the examples as well as the :ref:`Python bindings <python-bindings>`,
+create a directory ``build``. Inside that directory run ``cmake`` and ``make``::
+
+ mkdir build
+ cd build
+ cmake ..
+ make
+
+Instead of ``cmake``, one can run ``ccmake`` for a curses interface. The
+following configuration options are available. One can set them either through
+the curses interface or by passing a flag of the form ``-Doptimize:bool=on`` to
+``cmake``.
+
+ :debug: Turns on debugging compilation
+ :optimize: Turns on compiler optimizations (`on` by default)
+ :logging: Turns on logging facilities
+ :counters: Turns on various built-in counters
+
+Depending on the combination of debugging and optimization, a particular
+``CMAKE_CXX_FLAGS*`` is chosen.
+
+.. tip:: The default settings work fine unless you want to dive into the
+ library's internals with logging or study the performance of various
+ algorithms with counters.
+
+.. todo:: Write sections on logging and counters.
+
+Some parts of Dionysus understand the ``DEBUG_CONTAINERS`` definition which can
+be appended to ``CMAKE_CXX_FLAGS``. If set, the library will use GCC STL's
+debugging containers (from the ``std::__debug`` namespace defined in ``debug/*``
+header files). These containers return safe iterators (the kind that check
+whether they are singular when compared, or additionally whether they are out of
+bounds when dereferenced).
+
+.. todo:: ``ZIGZAG_CONSISTENCY`` definition
+
+
+Install
+-------
+
+At the moment there are no installation procedures. To run the Python code you
+need to have ``.../build/bindings/python`` somewhere in your ``PYTHONPATH``.
+I.e. add::
+
+ export PYTHONPATH=.../build/bindings/python
+
+to your ``~/.bashrc`` (assuming you are using Bash_). Alternatively, run the
+python examples from within ``.../build/bindings/python``::
+
+ python .../Dionysus/examples/triangle/triangle.py
+
+The C++ examples can be run from anywhere. The C++ library consists only of
+header files (no library actually needs to be built), so to compile against it,
+it suffices to add ``-I .../Dionysus/include`` to your ``g++`` flags::
+
+ g++ your-code.cpp -o your-code -I .../Dionysus/include
+
+Proper installation procedures (with ``make install``) will be added in the
+future.
+
+.. _Bash: http://www.gnu.org/software/bash/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/index.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,52 @@
+Welcome to Dionysus' documentation!
+===================================
+
+Dionysus is a C++ library for computing persistent homology. It provides
+implementations of the following algorithms:
+
+* Persistent homology computation [ELZ02]_ [ZC05]_
+* Vineyards [CEM06]_ |cpp-only|
+* Persistent cohomology computation (described in [dSVJ09]_) |cpp-only|
+* Zigzag persistent homology [CdSM09]_
+* :ref:`examples` provide useful functionality in and of themselves:
+
+ * :ref:`Alpha shape construction <alpha-shape-example>` in 2D and 3D
+ * :ref:`Rips complex construction <rips-example>`
+ * Cech complex construction |cpp-only|
+
+.. todo::
+ Document more examples.
+
+The C++ API is currently very poorly documented. One's best source for its
+documentation is its usage in various :ref:`examples` (located in
+:sfile:`examples/`).
+
+The :ref:`Python bindings <python-bindings>` provide both a simple interface to
+the low-level C++ functionality as well as high-level auxilliary routines. Their
+"thinness" is meant to provide the efficiency benefits of C++ together with the
+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
+
+
+..
+ Indices and tables
+ ==================
+
+ * :ref:`genindex`
+ * :ref:`modindex`
+ * :ref:`search`
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/alphashapes.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,42 @@
+.. _alphashapes:
+
+Alpha shapes
+============
+
+There are two functions provided to compute alpha shapes. One in 2D and one in
+3D. Both take a list of points (each a list of coordinates) as input, and fill a
+list with the simplices of the `Delaunay triangulation`_. Each such simplex is
+said to be *attached* (or *regular*) if its dual Voronoi cell does not contain a
+critical point of the distance function to the point set. The smallest value of
+the squared distance function on the dual Voronoi cell of the Delaunay simplex
+is the alpha shape value assigned to it. This value is stored in the simplex's
+`data` attribute; whether it is attached is stored in the `attached` attribute.
+
+.. _`Delaunay triangulation`: http://en.wikipedia.org/wiki/Delaunay_triangulation
+
+
+.. function:: fill_alpha2D_complex(points, complex)
+
+ Appends to the list `complex` the simplices of the 2D Delaunay triangulation
+ on the `points`.
+
+.. function:: fill_alpha3D_complex(points, complex)
+
+ Appends to the list `complex` the simplices of the 3D Delaunay triangulation
+ on the `points`.
+
+
+Example
+-------
+
+The following example generates 10 points on a circle, and computes their
+Delaunay triangulation with corresponding alpha shape values::
+
+ from math import sin, cos, pi
+ points = [[cos(2*pi*t/10), sin(2*pi*t/10)] for t in xrange(10)]
+ complex = []
+ fill_alpha2D_complex(points, complex)
+
+One can extract any given alpha shape with the usual Python list notation::
+
+ alphashape = [s for s in complex if s.data <= .5]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/filtration.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,33 @@
+:class:`Filtration` class
+=========================
+
+.. class:: Filtration
+
+ This class serves as a bridge between a complex represented as a
+ lexicographically sorted list of simplices, and the
+ :class:`StaticPersistence` class which needs to know the order in which the
+ simplices appear in the filtration.
+
+ .. method:: __init__(simplices, cmp)
+
+ Initializes :class:`Filtration` by internally arranging the indices of
+ the elements in the list `simplices` in the order sorted with respect to
+ `cmp`.
+
+ .. method:: __getitem__()
+
+ Random access to the elements of the filtration.
+
+ .. method:: __iter__()
+
+ Iterator over the elements of the filtration, which are simply the
+ indices of the simplices in the original list `lst` sorted with respect
+ to the comparison `cmp`. E.g.::
+
+ simplices = [Simplex([0], 2), ..., Simplex([3,4,5], 3.5)]
+ f = Filtration(simplices, data_dim_cmp)
+ for i in f: print simplices[i]
+
+ .. method:: __len__()
+
+ Size of the filtration.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/overview.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,22 @@
+.. _python-bindings:
+
+Python bindings: module :mod:`dionysus`
+=======================================
+
+.. module:: dionysus
+.. moduleauthor:: Dmitriy Morozov <dmitriy@mrzv.org>
+
+The :ref:`tutorial` describes how to use the bindings. The pages in this section
+document the API of various classes and functions.
+
+The following classes are available in the module:
+
+.. toctree::
+ :maxdepth: 1
+
+ simplex.rst
+ filtration.rst
+ static-persistence.rst
+ alphashapes.rst
+ rips.rst
+ zigzag-persistence.rst
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/rips.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,109 @@
+:class:`Rips` class
+======================
+
+.. class:: Rips
+
+ .. method:: __init__(distances)
+
+ Initializes :class:`Rips` with the given `distances` whose main purpose
+ is to return the distance of two points given their indices. See
+ Distances_ below.
+
+ .. method:: generate(k, max, functor[, seq])
+
+ Calls `functor` with every simplex in the `k`-skeleton of the Rips
+ complex :math:`VR` (`max`). If `seq` is provided, then the complex is
+ restricted to the vertex indices in the sequence.
+
+ .. method:: vertex_coface(v, k, max, functor[, seq])
+
+ Calls `functor` with every coface of the vertex `v` in the `k`-skeleton
+ of the Rips complex :math:`VR` (`max`). If `seq` is provided, then the
+ complex is restricted to the vertex indices in the sequence.
+
+ .. method:: edge_cofaces(u, v, k, max, functor[, seq])
+
+ Calls `functor` with every coface of the edge (`u`, `v`) in the
+ `k`-skeleton of the Rips complex :math:`VR` (`max`). If `seq` is
+ provided, then the complex is restricted to the vertex indices in the
+ sequence.
+
+ .. method:: cmp(s1, s2)
+
+ Compares simplices `s1` and `s2` with respect to their ordering in the
+ Rips complex. Note that like Python's built in `cmp` this is a three
+ possible outsome comparison (-1,0,1) for (:math:`\leq, =, \geq`,
+ respectively).
+
+ .. method:: eval(s)
+
+ Returns the size of simplex `s`, i.e. the length of its longest edge.
+
+
+.. _distances:
+
+Distances
+---------
+
+An instance of `distances` passed to the constructor of :class:`Rips` should
+know its length and the distances between the points. The length should be
+retrievable via ``len(distance)`` and it determines how many points the complex
+is built on. The distances between the points are inferred by the class
+:class:`Rips` by calling `distances` with a pair of vertices as arguments.
+
+For example, the following class represents 10 points on an integer lattice::
+
+ class Distances:
+ def __len__(self):
+ return 10
+
+ def __call__(self, x, y):
+ return math.fabs(y-x)
+
+The bindings provide a pure Python class :class:`PairwiseDistances` to deal with
+explicit points in a Euclidean space. It is defined in
+:sfile:`bindings/python/dionysus/distances.py`::
+
+ class PairwiseDistances:
+ def __init__(self, points, norm = l2):
+ self.points = points
+ self.norm = norm
+
+ def __len__(self):
+ return len(self.points)
+
+ def __call__(self, x, y):
+ return self.norm([x - y for (x,y) in zip(self.points[p1], self.points[p2])])
+
+Another distances class is available that speeds up the computation of the Rips
+complex at the expense of the memory usage: :class:`ExplicitDistances`. It is
+initialized with an instance of any class that behaves like a distances class,
+and it stores all of its distances explicitly to not have to recompute them in
+the future::
+
+ distances = PairwiseDistances(points)
+ distances = ExplicitDistances(distances)
+
+
+Example
+-------
+
+The following example reads in points from a file, and fills the list
+`simplices` with the simplices of the 2-skeleton of the Rips complex built on
+those vertices with distance cutoff parameter 50. Subsequently it computes the
+persistence of the resulting filtration (defined by ``rips.cmp``)::
+
+ points = [for p in points_file('...')]
+ distances = PairwiseDistances(points)
+ rips = Rips(distances)
+ simplices = []
+ rips.generate(2, 50, simplices.append)
+
+ f = Filtration(simplices, rips.cmp)
+ p = StaticPersistence(f)
+ p.pair_simplices()
+
+Essentially the same example is implemented in
+:sfile:`examples/rips/rips-pairwise.py`, although it reads the `k` and `max`
+parameters for the Rips complex on the command line, and uses a trick to speed
+up the computation.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/simplex.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,84 @@
+:class:`Simplex` class
+======================
+
+.. class:: Simplex
+
+ .. method:: __init__(seq[, data])
+
+ Initializes :class:`Simplex` from the given sequence `seq` and
+ optionally real value `data`, e.g.::
+
+ s = Simplex([1,2,3], 7.8)
+
+ .. method:: add(v)
+
+ Adds vertex `v` to the simplex, increasing its dimension by 1.
+
+ .. seealso:: :meth:`~Simplex.join`
+
+ .. attribute:: boundary
+
+ Iterator over the boundary :math:`\partial \sigma` of the simplex,
+ e.g.::
+
+ for sb in s.boundary: print sb
+
+ .. method:: contains(v)
+
+ :returns: `True` iff the simplex contains vertex `v`.
+
+ .. method:: dimension()
+
+ :returns: the dimension of the simplex (one less than its number of
+ vertices).
+
+ .. method:: join(other)
+
+ Joins the current simplex with the `other` simplex. The method copies over
+ the vertices from the `other` simplex.
+
+ .. attribute:: data
+
+ Real value stored in the simplex.
+
+ .. attribute:: vertices
+
+ (Sorted) vertices of the simplex accessible as a sequence, e.g.::
+
+ for v in s.vertices: print v,
+
+ .. method:: __hash__()
+ .. method:: __eq__(other)
+
+ Simplices are hashable, and equality comparable, and therefore can be
+ stored in a dictionary. Simplices are equal if their
+ :attr:`~Simplex.vertices` are the same.
+
+
+Utility functions for manipulating simplices
+--------------------------------------------
+
+The first function :func:`vertex_cmp` is a Python interface to a C++ function.
+The rest are pure Python functions defined in
+:sfile:`bindings/python/dionysis/__init__.py`.
+
+.. function:: vertex_cmp(s1, s2)
+
+ Compares the two simplices with respect to the lexicographic order of their vertices.
+
+.. function:: vertex_dim_cmp(s1, s2)
+
+ Compares the two simplices with respect to their dimension, and lexicographically
+ within the same dimension.
+
+.. function:: data_cmp(s1, s2)
+
+ Compares the two simplices with respect to the data (real values) they
+ store.
+
+.. function:: data_dim_cmp(s1, s2)
+
+ Compares the two simplices with respect to their dimension and within the same
+ dimension with respect to their data.
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/static-persistence.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,67 @@
+:class:`StaticPersistence` class
+================================
+
+.. class:: StaticPersistence
+
+ .. method:: __init__(filtration)
+
+ Initializes :class:`StaticPersistence` with the given
+ :class:`Filtration`. This operation effectively computes the boundary
+ matrix of the complex captured by the filtration with rows and columns
+ sorted with respect to the filtration ordering.
+
+ .. method:: pair_simplices()
+
+ Pairs simplices using the [ELZ02]_ algorithm.
+
+ .. method:: __call__(i)
+
+ Given an SPNode in the internal representation, the method returns its
+ integer offset from the beginning of the filtration. This is useful to
+ lookup the actual name of the simplex in the complex. For example, the
+ following snippet prints out all the unpaired simplices::
+
+ for i in persistence:
+ if i == i.pair: print complex[filtration[persistence(i)]]
+
+ .. method:: __iter__()
+
+ Iterator over the nodes (representing individual simplices). See
+ :class:`SPNode`.
+
+ .. method:: __len__()
+
+ Returns the number of nodes (i.e. the number of simplices).
+
+
+.. class:: SPNode
+
+ The class represents nodes stored in :class:`StaticPersistence`. These nodes
+ are aware of their :meth:`sign` and :attr:`pair` (and :meth:`cycle` if
+ negative after :meth:`StaticPersistence.pair_simplices` has run).
+
+ .. method:: sign()
+
+ Returns the sign of the simplex: `True` for positive, `False` for
+ negative.
+
+ .. attribute:: pair
+
+ Simplex's pair. The pair is set to self if the siplex is unpaired.
+
+ .. method:: cycle()
+
+ If the simplex is negative, its cycle (that it kills) is non-empty, and
+ can be accessed using this method. The cycle itself is an iterable
+ container of :class:`SPNode`. For example, one can print the basis for
+ the (bounding) cycles::
+
+ for i in persistence:
+ for ii in i.cycle(): print complex[filtration[persistence(ii)]]
+
+ .. method:: __eq__(other)
+
+ Returns true if the two nodes are the same. Useful for determining if
+ the node is unpaired (iff ``i == i.pair``), e.g::
+
+ print len([i in persistence if i == i.pair]) # prints the number of unpaired simplices
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/zigzag-persistence.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,93 @@
+:class:`ZigzagPersistence` class
+================================
+
+The class deals with the setting :math:`K_1 \rightarrow K_2 \leftarrow K_3 \rightarrow \dots`
+
+.. class:: ZigzagPersistence
+
+ .. method:: add(boundary, birth)
+
+ Adds a simplex with the given `boundary` to the complex, i.e.
+ :math:`K_{i+1} = K_i \cup \sigma` and `boundary` = :math:`\partial \sigma`.
+ If a new class is born as a result of the addition, `birth` is stored with
+ it for future reference.
+
+ :returns: a pair (`i`, `d`). The first element is the index `i`.
+ It is the internal representation of the newly added simplex,
+ and should be used later for removal or when constructing the
+ boundaries of its cofaces. In other words, `boundary` must
+ consist of these indices. The second element `d` is the death
+ element. It is `None` if a birth occurred, otherwise it
+ contains the value passed as `birth` to
+ :meth:`~ZigzagPersistence.add` or
+ :meth:`~ZigzagPersistence.remove` when the class that just
+ died was born.
+
+ .. method:: remove(index, birth)
+
+ Removes the simplex identified by the given `index` from the complex. If
+ a new class is born as a result of the removal, `birth` is stored with
+ it for future reference.
+
+ :returns: `None` if a birth has occurred, otherwise it contains the value
+ passed as `birth` to :meth:`~ZigzagPersistence.add` or
+ :meth:`~ZigzagPersistence.remove` when the class that just
+ died was born.
+
+The :ref:`triangle-zigzag-example` illustrates the use of :class:`ZigzagPersistence`.
+
+
+Auxilliary functions
+--------------------
+
+A pair of auxilliary functions is provided to help add and remove entire
+collections of simplices. Both are pure Python functions defined in
+:sfile:`bindings/python/dionysus/zigzag.py`.
+
+ .. function:: add_simplices(zigzag, simplices, complex, birth, report_local = False)
+
+ Adds each simplex in `simplices` to the `zigzag`. `complex` is a
+ dictionary mapping simplices to their indices (in `zigzag`'s internal
+ representation). All the newly born classes are given the value of
+ `birth`.
+
+ :returns: list of deaths that occur as a result of `simplices`' removal.
+ Each death is a pair of the dimension of the class and the
+ `birth` value passed when the class was born. By default the
+ deaths equal to `birth` are not reported unless `report_local`
+ is set to `True`.
+
+ .. function:: remove_simplices(zigzag, simplices, complex, birth, report_local = False)
+
+ Same parameters and return as in :func:`add_simplices` except that
+ `simplices` are removed from the `zigzag` and the `complex`.
+
+
+
+:class:`ImageZigzagPersistence` class
+=====================================
+
+The class deals with the setting
+
+.. math::
+ \begin{array}{ccccccc}
+ K_1 & \rightarrow & K_2 & \leftarrow & K_3 & \rightarrow & \dots \\
+ \uparrow & & \uparrow & & \uparrow & \\
+ L_1 & \rightarrow & L_2 & \leftarrow & L_3 & \rightarrow & \dots
+ \end{array}
+
+where the vertical maps are inclusions, i.e. :math:`L_i \subseteq K_i`.
+
+.. class:: ImageZigzagPersistence
+
+ .. method:: add(boundary, subcomplex, birth)
+
+ Interface is the same as in :meth:`ZigzagPersistence.add`. The
+ additional parameter `subcomplex` controls whether the simplex is added
+ to :math:`L` or not. We always have :math:`K_{i+1} = K_i \cup \sigma`.
+ If `subcomplex` is true, then :math:`L_{i+1} = L_i \cup \sigma`,
+ otherwise :math:`L_{i+1} = L_i`.
+
+ .. method:: remove(index, birth)
+
+ Interface is exactly the same as in :meth:`ZigzagPersistence.remove`.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/substitutions.aux Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,2 @@
+.. |cpp-only| replace:: :sup:`(C++ only)`
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/tutorial.rst Mon May 11 12:45:49 2009 -0700
@@ -0,0 +1,197 @@
+.. _tutorial:
+
+Brief Tutorial
+==============
+
+It suffices to import only the necessary commands from module :mod:`dionysus`,
+but for the purposes of the tutorial, we import everything::
+
+ from dionysus import *
+
+Read in a points file (for :func:`points_file` the points should appear one per
+line with coordinates separated with spaces; of course, one can read in the
+points however one wants)::
+
+ points = [p for p in points_file('points_filename')]
+ print "Number of points:", len(points)
+
+
+Complexes
+---------
+
+If the points are in :math:`\mathbb{R}^2` or :math:`\mathbb{R}^3`, one can
+construct an alphashape filtration::
+
+ simplices = []
+ fill_alpha2D_complex(points, simplices) # for 2D, or
+ fill_alpha3D_complex(points, simplices) # for 3D
+
+
+Functions :ref:`fill_alpha*_complex <alphashapes>` fill the ``simplices`` list
+with all the :class:`simplices <Simplex>` of the Delaunay triangulation. Each one has attributes
+``data`` and ``attached`` set to its respective value in the alpha shape (the
+smallest value of the squared distance function on the dual Voronoi cell) and
+whether the simplex is attached or not (i.e. whether its dual cell does not or
+does contain a critical point of the distance function). See :ref:`alphashapes`
+for more details.
+
+As a result, if one wanted only those simplices whose alpha shape value did not
+exceed 10, one could obtain them as follows::
+
+ simplices10 = [s for s in simplices if s.data <= 10]
+
+If the point set lies in higher dimensions, one may construct a Rips complex on
+it. This complex requires only pairwise distances, which makes it very
+versatile. One must first construct an instance of a class providing such
+distances (e.g. :class:`PairwiseDistances` for explicit points, see
+:ref:`distances` for more details), and then pass it to the :class:`Rips`
+complex class::
+
+ distances = PairwiseDistances(points)
+ rips = Rips(distances)
+
+Usually, because of space restrictions, generation of a Rips complex has to be
+restricted to a :math:`k`-skeleton of the complex and some maximal parameter
+:math:`max`. In the following example :math:`k = 3` and :math:`max = 50`::
+
+ simplices = []
+ rips.generate(3, 50, simplices.append)
+
+:meth:`Rips.generate` takes a skeleton and a maximum distance cutoffs, and a
+function which is called with each of the simplices in the complex (in this
+case, Python list `simplices`' ``append`` method).
+
+
+Persistence
+-----------
+
+There are two ways of computing persistence in Dionysus. The first *offline*
+way, required by :class:`StaticPersistence` (and its derivatives), is to set up
+the entire filtration at once, compute its persistence in one operation, and
+then examine the pairings. The second way is to feed simplices one by one in an
+*online* manner and manually keep track of the pairs that are created.
+:class:`ZigzagPersistence` and `CohomologyPersistence` |cpp-only| accept their
+input this way,
+
+
+Offline
+^^^^^^^
+
+For the first approach, i.e. to use :class:`StaticPersistence`, one must put the
+simplices into lexicographic ordering (sort them with respect to
+:func:`vertex_cmp`), and construct a filtration with respect to some ordering
+(for example, :func:`data_dim_cmp` for alpha shapes or :meth:`Rips.cmp` for the
+Rips complex)::
+
+ simplices.sort(vertex_cmp)
+ f = Filtration(simplices, data_dim_cmp) # for the alpha shapes
+ f = Filtration(simplices, rips.cmp) # for the rips complex
+
+Creating an instance of :class:`StaticPersistence` initialized with the
+filtration really initializes a boundary matrix. The subsequent call to
+:meth:`~StaticPersistence.pair_simplices` reduces the matrix to compute the
+persistence of the filtration::
+
+ p = StaticPersistence(f)
+ p.pair_simplices()
+
+Once the simplices are paired, one may examine the pairings by iterating over
+the instance of :class:`StaticPersistence`::
+
+ for i in p:
+ if i.sign():
+ birth = simplices[f[p(i)]]
+ if i == i.pair:
+ print birth.dimension(), birth.data, "inf"
+ continue
+
+ death = simplices[f[p(i.pair)]]
+ print birth.dimension(), birth.data, death.data
+
+The iterator ``i`` over the instance ``p`` of :class:`StaticPersistence` is of type
+:class:`SPNode`, and represents individual simplices taken in the filtration
+order. It knows about its own :meth:`~SPNode.sign` and :attr:`~SPNode.pair`
+(``i.pair == i`` if the simplex is unpaired). Calling ``p`` with the iterator
+``i`` (``p(i)``) gives its integer index in the filtration, which can then be
+used with the instance ``f`` of :class:`Filtration` to get the integer index in
+the lexicographically sorted list `simplices`: ``simplices[f[p(i)]]``. The
+previous code snippet prints out the persistence diagrams of the given
+filtration.
+
+
+Online
+^^^^^^
+
+Class :class:`ZigzagPersistence` accepts additions and removals of the simplices
+one by one, and returns an internal representation of the simplices. When one
+adds a simplex via :meth:`ZigzagPersistence.add`, one must provide its boundary
+in this internal representation together with a *birth value* which
+:class:`ZigzagPersistence` will store in case a class is born as a result of the
+addition. :meth:`~ZigzagPersistence.add` returns a pair ``(i,d)``. ``i`` is the
+internal representation of the newly added simplex, which one must record for
+future use in the boundaries of its cofaces (below it is recorded in the
+dictionary ``complex``). ``d`` is `None` in case of the birth, otherwise, it is the
+previously recorded birth value of the class that dies as a result of the
+addition. The following code adds all the ``simplices`` to a zigzag::
+
+ simplices.sort(data_dim_cmp)
+ complex = {}
+ zz = ZigzagPersistence()
+ for s in simplices:
+ i,d = zz.add([complex[sb] for sb in s.boundary, (s.dimension(), s.data))
+ complex[s] = i
+ if d is not None: # we have a death
+ dimension, birth = d # we previously stored the (dimension, data) pair
+ print dimension, birth, s.data
+
+Similarly, one can remove simplices by calling :meth:`ZigzagPersistence.remove`
+with the internal index previously returned by :meth:`~ZigzagPersistence.add`::
+
+ for s in reversed(simplices):
+ d = zz.remove(complex[s], (s.dimension() - 1, s.data))
+ del complex[s]
+ if d is not None:
+ dimension, birth = d
+ print dimension, birth, s.data
+
+Naturally, :meth:`~ZigzagPersistence.remove` returns only `d`.
+
+If one wants to add or remove an entire set of simplices "at once" (i.e. with
+all births being assigned the same value), there are two helper functions:
+:func:`add_simplices` and :func:`remove_simplices`.
+
+See the :ref:`bindings reference <python-bindings>` for more details, and
+:ref:`triangle-zigzag-example` for an example.
+
+
+.. _speed-up-suggestions:
+
+Speed-up suggestions
+--------------------
+
+Currently, when the choice comes between efficiency and flexibility, the Python
+bindings err on the side of flexibility. There is hope that in the future the
+choice won't really be necessary. Meanwhile, one can use a few techniques that
+speed up computation at the expense of memory:
+
+* To avoid (possibly expensive) computation of distances during Rips complex
+ generation, store :class:`ExplicitDistances` (see :ref:`distances`)::
+
+ distances = PairwiseDistances(points)
+ distances = ExplicitDistances(distances)
+
+* To avoid the computation of simplex sizes in the Rips complex during the
+ initialization of a :class:`Filtration`, store them explicitly in
+ :attr:`Simplex.data` attribute (this is not done by default to save memory);
+ then use :func:`data_dim_cmp` in the initialization of the
+ :class:`Filtration`::
+
+ rips = Rips(distances)
+ simplices = []
+ rips.generate(..., simplices.append)
+ for s in simplices: s.data = rips.eval(s)
+ f = Filtration(simplices, data_dim_cmp)
+
+
+
+.. include:: substitutions.aux