Update on class structure, lattice constructors, plotting

Sorry that I haven't posted to this blog for a while. That doesn't mean that nothing happened in the meanwhile—on the contrary, I got quite a lot done in the last weeks of this Summer of Code.

First of all, I completely redesigned the class structure for lattices in Sage. Accounting for the discussion on sage-devel, a lattice in Sage is now a FreeQuadraticModule_ambient_pid with a given inner product matrix ("Gram matrix") and the "standard" basis [(1, 0, 0, ...), (0, 1, 0, ...), ...]. Additionally, a basis for the embedding of the lattice is stored. When constructing a lattice, either a Gram matrix or an embedded basis can be given, and the respective other property will be calculated (either using Cholesky decomposition or via the standard real inner product).

A lattice can also be constructed using a quadratic form or an order in an algebraic number field. The inner product matrix will be calculated accordingly.

There are only two lattice classes now: Lattice_with_basis is meant to be the base class for all lattices in Sage, while Lattice_ZZ represents lattices with base ring ZZ and integral Gram matrix. Storing the Gram matrix as well and not only the basis has the advantage that some lattices with floating-point basis can still be represented exactly, e.g., the embedding of ZZ[sqrt(2)] as a lattice in the sense of Minkowski's geometry of numbers. LLL reduction is primarily performed on the Gram matrix, and the embedded basis is transformed accordingly.

I also added new ways of creating lattices: random_lattice uses random_matrix to create a random basis of a given dimension and turns it into a ZZ-lattice. special_lattice can be used to create particular "special" named lattices such as the Leech lattice. The list of integrated named lattices is probably not complete, but easily extensible. Again, lattices can be specified using their basis or their Gram matrix.

Finally, 2-dimensional and 3-dimensional lattices can be plotted now. All combinations of basis vectors are plotted as points, and the extensive graphics capabilities of Sage allow to modify the plot further. For instance, the Voronoi cell can easily be added, too:

This graphics was constructed using the following code:

L = special_lattice('SimpleCubic')
graphics = plot_lattice(L)
V = L.voronoi_cell()
graphics += V.plot()
graphics.show(viewer='tachyon')

How do you like it?

Rebasing to Sage 5.2

As there were some changes to Sage since version 5.0 that also affect my code, I rebased my code from Sage 5.0 to 5.2. That meant installing the latest version, cloning a new branch "lattices", and manually copying over the lattices files and re-doing respective changes to some global files. I retained the 5.0-compatible version as a separate branch, whereas the main branch on github is 5.2-targeted. Of course, it would have been much easier if all of Sage was managed through git, so I'm really looking forward to that happening.

Closest vectors

I haven't written about the latest development yet: closest vectors in lattices. This is done using the algorithm in the paper "A Deterministic Single Exponential Time Algorithm for Most Lattice Problems based on Voronoi Cell Computations" by Daniele Micciancio and Panagiotis Voulgaris. As the title suggests, it uses the Voronoi cell V. The idea is

1. to deal with the special case where the target vector is in 2V, and
2. to solve the general problem with polynomially many calls to this special case.

See the paper (esp. Algorithm 1) for the details.

Currently, I am working on polishing up documentation, and maybe we still need to change the class structure a little.

Voronoi cells

Recently, I have implemented the diamond-cutting algorithm to calculate the Voronoi cell of a lattice as described in the paper Computing the Voronoi Cell of a Lattice: The Diamond-Cutting Algorithm by Emanuele Viterbo and Ezio Biglieri.

Basically, it starts with a parallelotope defined by hyperplanes through (and normal to) the halves of the basis vectors (and their negated counterparts). Then it gradually combines basis vectors and cuts off  further half spaces, much like a diamond is cut.

The algorithm considers vectors within a certain distance from the origin. This radius can be specified explicitly to speed things up significantly; if not, a "conservative" safe choice is made.

In the current implementation, I am using Sage's existing Polyhedron class to perform the cutting. It seems reasonably fast.

Emanuele Viterbo was so friendly as to provide his implementation of the algorithm from 1995. Comparing running times, it does not seem like it is worth switching from the Polyhedron class to a custom-tailored version of cutting (i.e. intersecting polyhedrons with half spaces). For example, calculating the Voronoi cell of the following lattice:

2.45 0.00 0.00 0.00 0.00 0.00 
-0.41 2.42 0.00 0.00 0.00 0.00 
-0.41 -0.48 2.37 0.00 0.00 0.00 
-0.41 -0.48 -0.59 2.29 0.00 0.00 
-0.41 -0.48 -0.59 -0.76 2.16 0.00 
-0.41 -0.48 -0.59 -0.76 -1.08 1.87

with ball radius 3 takes 51.4 sec using Emanuele's program but only 13.5 sec with my implementation using Polyhedron on my 2.4 GHz Intel Core 2 Duo Mac. (104 cutting operations are performed in both cases.)

A nice consequence of using the Polyhedron class is that it already provides methods to deal with the resulting Voronoi cell. For example, a representation either by (in)equalities or by vertices/edges/faces can be computed easily.

I had to take special care to deal with non-quadratic basis matrices. To deal with "over-specified" lattices (e.g. redundant basis vectors as in (1,0), (2,0), (0,1)), an LLL-reduction is performed prior to diamond-cutting. (This helps to speed up the algorithm anyway.)

"Non-full" lattices (less basis vectors than dimension of embedded space) are dealt with in a special way, too:

1. To get a quadratic basis matrix, "artificial" basis vectors that represent infinity are added first (they are simply chosen to be longer than any other basis vector);
2. LLL-reduction is performed to get a quadratic basis matrix.
3. After diamond-cutting, all inequalities induced by artificial basis vectors are removed again.

This allows the following:

sage: L = Lattice([[2, 0, 0], [0, 2, 0]])
sage: V = L.voronoi_cell()
sage: V.Hrepresentation()
(An inequality (-1, 0, 0) x + 1 >= 0, An inequality (0, -1, 0) x + 1 >= 0, An inequality (1, 0, 0) x + 1 >= 0, An inequality (0, 1, 0) x + 1 >= 0)

Of course, this is still sort of a workaround—maybe there is a better solution to the problem of non-quadratic basis matrices.

Update on basis coercion, LLL and Fincke-Pohst

As written in a previous post, I tried to prevent FreeModule from coercing its basis vectors to the fraction field of its ambient space, mainly motivated by the fact that I wanted floating-point numbers to remain "inexact" and not being converted to rational numbers. However, the resulting discussion on sage-devel made clear that this is actually desirable and we do want this coercion to happen. This is a pity somehow because it renders some of my code useless, but admittedly, it was a "hack" anyway.

So now Lattice_with_basis inherits from FreeModule_submodule_with_basis_pid without changing much of its behaviour. Only echelonization of the basis matrix is turned off, and the representation (_repr_) is different, of course.

Most of the action happens in Lattice_ZZ_in_RR now: I added code to reduce the given basis using the LLL algorithm. Because of our restriction to QQ (no RR vectors), the patch for Ticket #12051 is already enough. (Nevertheless, I might still add an LLL implementation for RR matrices.) Basis reduction can be disabled by handing reduce_basis=False to the Lattice constructor/factory function.

Furthermore, I added a shortest_vectors method to Lattice_ZZ_in_RR that uses Pari's Fincke-Pohst implementation qfminim to find (a certain number of) short vectors (of a given maximum length). The builtin method gram_matrix in FreeModule is used to get the basis inner product matrix. I might add a generator version as well at a later point.

The lattice discriminant is already implemented by FreeModule using the appropriate Gram matrix as well.

Now I'm starting to work on Voronoi cells based on the diamond-cutting algorithm.

My Eclipse/PyDev/git workflow for Sage

I want to briefly describe how I code for Sage so that others can comment on how to improve the workflow or maybe learn something from it.

For Python (and actually for all programming languages, meanwhile), I've always been using Eclipse as an IDE. With PyDev you get pretty nice Python syntax highlighting, code outlines, and features such as code completion, although I have never used that a lot. PyDev has even added Cython support for .pyx files recently—yeah!

I have added the complete Sage tree (which is located in /Applications/sage/ on my Mac) as an Eclipse project. As recommended and described earlier, I created a new branch ("sandbox") using sage -clone lattices; consequently, I'm working on the files in sage/devel/sage-lattices.

Test-driven development makes much sense in general, and especially when working on Sage:

• Due to the fact that you have to rebuild Sage every time you change something, working completely "interactive" isn't possible.
• You don't want to type the code to "bootstrap" your tests (e.g. creating lattices) every single time you change something.
• Eventually, lots of examples are demanded in the documentation anyway.

This is why I usually add some examples that I want to be working in Python docstrings, then code something, and then run the tests. Running the tests is done by

./sage -bt devel/sage-lattices/sage/lattices/lattice.py

which will rebuild the changed components of Sage (compiling Cython code if necessary) and then run the tests defined in lattice.py. Usually it says "All tests passed!" in the end; otherwise it gives an overview of the test cases (examples) that failed.

To see how the documentation that is generated from docstrings looks like, I sometimes rebuild it. For some strange reasons, the documentation generator (Sphinx) does not always recognize that I changed something in the source files. To avoid having to rebuild the complete reference manual (which takes quite a while) every time, I created a small documentation part that only includes the lattices chapter. This is easily done by adding a folder to doc/en and creating files conf.py, index.rst, and lattices.rst therein. The latter basically contains

Lattices
========

.. automodule:: sage.lattices.lattice
   :members:
   :undoc-members:
   :show-inheritance:

Then, the relevant documentation can be rebuilt using

sage --docbuild lattices html -S -a,-E

(-S passes arguments on to Sphinx, -a,-E erase and recreate everything.)

On a side note: When building the HTML documentation on my Mac, I got a bunch of warnings saying "Incompatible library version: libfontconfig.1.dylib requires version 13.0.0 or later, but libfreetype.6.dylib provides version 10.0.0". It seems that the system-wide version of dvipng that I installed through Fink referenced a library "libfreetype" that was loaded through Sage, but with an insufficiently recent version. I simply resolved this by replacing the libfreetype* files in sage/local/lib with the Fink files from /sw/lib. Now the HTML documentation displays beautiful PNG graphics for LaTeX formulas.

I tried to use the integrated Eclipse/PyDev Python debugger with Sage, but ran into some problems. First, you have to launch Eclipse through the Sage shell, i.e. using

sage -sh
/path/to/eclipse

But then, after selecting sage-python as interpreter, when I tried to run the doctests in debug mode, I still got some weird ImportErrors. I could reproduce them using sage-python alone (leaving Eclipse out) and reported the issue on sage-trac.

Despite test-driven development, sometimes I still want to play around with things in Sage interactively.  To reload my lattice module after having done

from sage.lattices import lattice

at one point, I still have to rebuild everything (using the mentioned sage -bt ... command which tests things as well), and then reload it using

lattice = reload(lattice)

Then I can go on creating lattices, e.g. by

L = lattice.Lattice([[1, 0], [0, 1]])

and experiment with them.

When I am done working on something, I do a git commit to my sage repository using

git commit -am "my commit message"

Using the Eclipse git interface is terribly slow when working on the entire Sage tree (even though only a few files are checked in), so I'm using the terminal for this. Usually I also push commits to the github repository using

git push

immediately.

Usually I have at least three terminal windows open: one entirely for running the doctests, one for handling git, and one for "experiments" in Sage. To finally include a picture on this blog as well, here's how this might look like:

Please let me know if you have any suggestions.

Finally, some useful documentation links:

• Sage Developer's Guide
• Editing and building documentation
• Unit testing for Sage objects

Class structure

The class structure of the upcoming lattice module is now as follows: Lattice_with_basis inherits from FreeModule_submodule_with_basis_pid (with some tweaks as described in my previous post) and is the general class for all lattices with bases. This is what Robert Miller called SubLatticeModule in a rather old discussion about lattices, although it is more general as it does not restrict to ZZ-lattices. Would it make sense to implement a general LatticeModule (in Robert Miller's notion) that does not have a concrete basis? How would such a LatticeModule be constructed and what could it offer?

Inheriting from Lattice_with_basis there is Lattice_ZZ that specializes on lattices with integer coefficients (ZZ-Lattices), probably the most usual form of lattices.

Furthermore, there is the specialization on lattices embeddable in the real vector space, Lattice_ZZ_in_RR (again, a very common form of lattices). The actual type of the vectors in the lattice could still be ZZ^n or QQ^n (not necessarily RR^n). Most of the algorithms implemented in this project will probably lie in this class.

I have already added a few lines of documentation and some examples/doctests.

You can see the current state of the code in the main lattice.py file in the code repository on github.