Wednesday, June 20, 2012

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.

Monday, June 4, 2012

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/
which will rebuild the changed components of Sage (compiling Cython code if necessary) and then run the tests defined in 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, index.rst, and lattices.rst therein. The latter basically contains

.. automodule:: sage.lattices.lattice
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
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

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:

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 file in the code repository on github.

Preventing FreeModule from coercing to its base_ring

I've been working on the class structure for lattices recently. As described in the project description and already suggested in a previous discussion on lattices in Sage, the generic Lattice class should inherit from FreeModule in some way. In our current approach, lattices will always have a basis, which is why the lattice base class Lattice_with_basis should probably inherit from FreeModule_submodule_with_basis_pid.

However, the issue with this is that, apparently, FreeModule expects the basis to be in the fraction field of the coefficient ring (if not in the coefficient ring itself). (The coefficient ring is called base_ring in FreeModule, which is a little confusing because this is not the ring the basis has to lie in.) For instance, a ZZ-lattice in RR^n would be constructed as (an inherited version of)
FreeModule(ZZ^n, <basis matrix>)
with a real-valued basis matrix, but then the basis matrix would be coerced to a matrix over QQ^n first, discarding the information about the original field that the lattice was embedded in. This is not desired here—we want to keep the original vector space, e.g. keep calculating with floating point numbers in RR.

After playing around with the way FreeModules are constructed, I implemented a small workaround: When calling the inherited constructor from FreeModule_submodule_with_basis_pid, echelonization and checking (whether the basis vectors are actually in the free module) are disabled to keep the basis "as is". Consequently, they have to be converted from their given form (usually lists) to elements of a vector space; this is already done in the factory function Lattice.

Furthermore, the element_class has to be overridden, because it is used in the constructor to convert the basis elements (even when echelonization and checking are disabled). Apparently, even the semi-private member _element_class has to be set because it is used instead of the interface function element_class, e.g., when random_element is called. Nevertheless, the function element_class in free_module has to be used to get the actual element class (the underlying ring being changed though), because it wraps the Element structure around the ring.

Now, the following works with Lattice:
sage: Lattice([[1.0, 0, 0], [0, 1, 0]], ZZ)
Real-embedded integer lattice of degree 3 and rank 2
Basis matrix:
[ 1.00000000000000 0.000000000000000 0.000000000000000]
[0.000000000000000  1.00000000000000 0.000000000000000]
While using FreeModule coerces to rationals:
sage: FreeModule(ZZ, 3).span([[1.0, 0, 0], [0, 1, 0]])
Free module of degree 3 and rank 2 over Integer Ring
Echelon basis matrix:
[1 0 0]
[0 1 0]
I know this workaround is not an optimal solution—maybe FreeModule itself should be modified if this behaviour is desired. Please let me know if you have any ideas.

Update: This still doesn't work. I came across an error (a pretty severe one, actually) when running the TestSuite on lattices, namely:
sage: L = Lattice([[i, 0], [0, 1]])
sage: + L.an_element()

Unhandled SIGSEGV: A segmentation fault occurred in Sage.
This probably occurred because a *compiled* component of Sage has a bug
in it and is not properly wrapped with sig_on(), sig_off(). You might
want to run Sage under gdb with 'sage -gdb' to debug this.
Sage will now terminate.
/Applications/sage/spkg/bin/sage: line 312:  4708 Segmentation fault
So while constructing lattices this way works perfectly, working with their elements does not. It seems elements still rely on L.basis_ring() being their own "parent" to some extent, which is not true e.g. for ZZ-lattices in RR^n. I really don't know whether the whole idea of subclassing from FreeModule this way is what FreeModule was designed for, so I asked on sage-devel.