Huge patch on Trac 7477: Matroid theory

[Stefan van Zwam]

Dear Sage developers,

Over the past two years I've been working with several people on bringing Matroid Theory to Sage. We wanted this to be research-level code that will remain relevant for a long time, so we decided to have the design mature for a while before integrating it with Sage. By now we have test-driven the code (it produces output consistent with the theory and with computations carried out by other software for matroid theory), and it is already in use by people in the matroid theory community who install it separately on top of Sage.

The code has now been submitted, see

http://trac.sagemath.org/sage_trac/ticket/7477

There are about 21,000 lines of code, and 500 functions with docstrings. While it is doable to check that it plays nice with the rest of Sage and follows all coding conventions (I'm  sure it does), a full code review sounds like a hard task. So I write you to solicit opinions on how to proceed. Four suggestions have been made on the ticket:

1) Treat it as if it were an spkg, and accept with only a cursory review

2) Try to assemble a team of reviewers, each responsible for part of the code

3) Chop it up into separate, smaller pieces and have those submitted as separate tickets

4) Get the authors of the package to review each others' code

Options (3) and (4) are not very feasible. The core part of the code has close links all over the place, and I think I can't get it much below 9,000 lines of code without a lot of extra work. It would also mean rewriting all docstrings since we can't use the catalog of matroids for our examples. As for (4), Rudi Pendavingh and I are responsible for the bulk of this code, and each of us has had a hand in pretty much every part of it. 

That leaves (1) and (2). I hope to get the other users/developers to review parts of it, but so far I have not heard of any commitments, and I fear it would still take a long time to get it accepted. I'd rather see that time spent building a user base and ironing out any bugs uncovered in actual use.

What do you all think?

Regards,

Stefan.

P.S. The impact on the rest of Sage is the following:

* add a directory sage.matroids

* lazy_import two items, a function "Matroid" and a package "matroids", on startup

* add an entry "Matroid theory" to the reference manual.

[David Joyner]

On Thu, May 2, 2013 at 2:40 PM, Stefan van Zwam <stefan...@gmail.com> wrote:

Dear Sage developers,

Over the past two years I've been working with several people on bringing Matroid Theory to Sage. We wanted this to be research-level code that will remain relevant for a long time, so we decided to have the design mature for a while before integrating it with Sage. By now we have test-driven the code (it produces output consistent with the theory and with computations carried out by other software for matroid theory), and it is already in use by people in the matroid theory community who install it separately on top of Sage.

The code has now been submitted, see

http://trac.sagemath.org/sage_trac/ticket/7477

There are about 21,000 lines of code, and 500 functions with docstrings. While it is doable to check that it plays nice with the

Thank you very much for this addition.

 

rest of Sage and follows all coding conventions (I'm  sure it does), a full code review sounds like a hard task. So I write you to solicit opinions on how to proceed. Four suggestions have been made on the ticket:

1) Treat it as if it were an spkg, and accept with only a cursory review

How easy/feasible would it be to create an spkg and propose it to be added as an optional package?

Is that another possibility?

 

2) Try to assemble a team of reviewers, each responsible for part of the code

3) Chop it up into separate, smaller pieces and have those submitted as separate tickets

4) Get the authors of the package to review each others' code

Options (3) and (4) are not very feasible. The core part of the code has close links all over the place, and I think I can't get it much below 9,000 lines of code without a lot of extra work. It would also mean rewriting all docstrings since we can't use the catalog of matroids for our examples. As for (4), Rudi Pendavingh and I are responsible for the bulk of this code, and each of us has had a hand in pretty much every part of it. 

That leaves (1) and (2). I hope to get the other users/developers to review parts of it, but so far I have not heard of any commitments, and I fear it would still take a long time to get it accepted. I'd rather see that time spent building a user base and ironing out any bugs uncovered in actual use.

What do you all think?

Regards,

Stefan.

P.S. The impact on the rest of Sage is the following:

* add a directory sage.matroids

* lazy_import two items, a function "Matroid" and a package "matroids", on startup

* add an entry "Matroid theory" to the reference manual.

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
To post to this group, send email to sage-...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Volker Braun]

The obvious solution would have been to not produce a patch bomb ;-)  You can have circular dependencies between different tickets, so you can make one for the library of matroids and the other for the core implementation and they depend on each other.

As I said on the ticket, you should be able to review the mathematical correctness among yourselves. You just need somebody who is more familiar with Sage to check that it integrates well.

[Stefan]

How easy/feasible would it be to create an spkg and propose it to be added as an optional package?

Is that another possibility?

It is a possibility, and would be pretty feasible, but an inferior solution. From our point of view (and our audience), having the code be an actual part of Sage is just more convenient. It makes sure we have full access to the bug reporting of Trac, and it lowers the threshold for contributions by other developers. Also, it is harder for other parts of the Sage library (graphs, linear codes, ...) to use the matroid routines when they are hidden in an optional package.

--Stefan.

[Nathann Cohen]

> It is a possibility, and would be pretty feasible, but an inferior solution.
> From our point of view (and our audience), having the code be an actual part
> of Sage is just more convenient. It makes sure we have full access to the
> bug reporting of Trac, and it lowers the threshold for contributions by
> other developers. Also, it is harder for other parts of the Sage library
> (graphs, linear codes, ...) to use the matroid routines when they are hidden
> in an optional package.

Annnnnd future additions to the SPKG will not have to be reviewed
either, and the documentation does not appear on Sage's website.

And they cannot really say to everybody that they can just download
the SPKG and use it in Sage, for if they attempt to use new features
of Sage with their code somebody who runs an old version of Sage can
install their SPKG, which uses functions that are not included in his
version.

I would vote for 1 among all alternatives... and trust the authors. As
simple as that. I also do believe that considering the code they
produced, they would be the first to fix bugs when they will be
reported, be it a small problem or a larger one. It seems to be a work
of care and good will :-)

Nathann

[Rob Beezer]

Dear sage-devel,

Late to the party, as usual.  I've been following the  sage-matroid  list since its inception, helping out where I can.  Stefan, Rudi and others, have been *very* responsive to bug reports on their work-in-progress.  The current users are a very knowledgeable bunch (both in mathematics and combinatorial computing), so the discussions about design decisions have been very thoughtful.  I have a lot of confidence in the quality of this contribution.

While I had hoped this did not become a "patch bomb", I do agree with much of what Nathann has said.

(1) The code is being used currently for research problems, so subtle bugs are being actively ferreted out.  These are being fixed rapidly.

(2) The (computational) matroid community will be a good addition to Sage's stable.

(3) Development will be much faster and broader when this work is integrated into Sage proper.  Much of the traffic on  sage-matroid  involves pecularities of implementing/running this as an add-on.

(4) Stressing other parts of Sage (eg matrix code) should lead to improvements there.

So I would be in favor of accepting this en masse if it had a knowledgeable review for fitting into Sage, as Nathann and Volcker propose.  I also hope that the matroid community can commit to working with other developers in improving the matrix code (and perhaps other areas) to the point where they do not need to reimplement existing functionality.  We are already seeing some of this at http://trac.sagemath.org/14627, so hopefully that can continue.

Finally, thanks to the  sage-matroid  group for all their hard work on this.  It is going to be another jewel in Sage's crown.

Rob

[kcrisman]

Finally, thanks to the  sage-matroid  group for all their hard work on this.  It is going to be another jewel in Sage's crown.

Hear, hear!  As long as it's not a Silmaril...

[Dima Pasechnik]

On Friday, 24 May 2013 04:31:46 UTC+8, Rob Beezer wrote:

So I would be in favor of accepting this en masse if it had a knowledgeable review for fitting into Sage, as Nathann and Volcker propose. 

I second this, too.

Dima

 

[Simon King]

I third this, three :)

In other words, +1.

[Volker Braun]

On Thursday, May 23, 2013 9:31:46 PM UTC+1, Rob Beezer wrote:

So I would be in favor of accepting this en masse if it had a knowledgeable review for fitting into Sage, as Nathann and Volcker propose.

I agree with everything except the spelling of my name ;-)

 

[Nathann Cohen]

> I agree with everything except the spelling of my name ;-)

Sooooooooooo how do we do ? We check that the ticket is up to date,
does not break doctests, and we set it to "positive review" with as
reviewer "The Sage Team" ? :-P

Otherwise they may have to rebase their patch every second day.

Nathann

[Volker Braun]

I think the following needs to be addressed:

  * segmentation faults due to recent changes in Sage

  * either replace the private reimplementation of matrices or give a good reason for why it is necessary

Then it is imho ready for inclusion.

[David Joyner]

On Thu, May 23, 2013 at 4:31 PM, Rob Beezer <goo...@beezer.cotse.net> wrote:
> Dear sage-devel,
>

...

>
> So I would be in favor of accepting this en masse if it had a knowledgeable
> review for fitting into Sage, as Nathann and Volcker propose. I also hope

+1

..

>
> Rob
>
>

[William Stein]

And it has 100% doctest coverage.  (It does, right?)

> ..
>
> >
> > Rob

[Michael Welsh]

On 25/05/2013, at 1:09 AM, William Stein <wst...@gmail.com> wrote:
>
> And it has 100% doctest coverage. (It does, right?)

Yes, it does. We made sure of that before submitting.

[Travis Scrimshaw]

Hey,

  I'm for accepting en masse since it is (mostly) self-contained, provided the doc builds cleanly too. From a skim-through of the patch, it seems like there is some documentation formatting I'd like to see changed (and some things that might be misformatted). However we can fix this up over time since this is minor compared to the functionality.

On Friday, May 24, 2013 2:22:35 AM UTC-7, Volker Braun wrote:

I think the following needs to be addressed:

  * segmentation faults due to recent changes in Sage

  * either replace the private reimplementation of matrices or give a good reason for why it is necessary

Then it is imho ready for inclusion.

Wasn't the second point addressed in http://trac.sagemath.org/14627, so shouldn't the matrices be replaced?

Best,

Travis

PS - the extra c in Volcker makes your name sound ckooler :P

[Volker Braun]

On Saturday, May 25, 2013 2:39:36 AM UTC+1, Travis Scrimshaw wrote:

  * either replace the private reimplementation of matrices or give a good reason for why it is necessary 

Wasn't the second point addressed in http://trac.sagemath.org/14627, so shouldn't the matrices be replaced?

Well I hope so but there might be another performance bottleneck that the matroid code hits. Maybe one of the authors can tell us more about it?

PS - the extra c in Volcker makes your name sound ckooler :P

Oh boy ;-)

 

[Stefan]

On Saturday, May 25, 2013 6:38:00 AM UTC-4, Volker Braun wrote:

On Saturday, May 25, 2013 2:39:36 AM UTC+1, Travis Scrimshaw wrote:

  * either replace the private reimplementation of matrices or give a good reason for why it is necessary 

Wasn't the second point addressed in http://trac.sagemath.org/14627, so shouldn't the matrices be replaced?

Well I hope so but there might be another performance bottleneck that the matroid code hits. Maybe one of the authors can tell us more about it?

I'd have to check carefully. There's likely to be some performance loss even when all obvious bottlenecks are accounted for, because for our special classes BinaryMatrix, TernaryMatrix, QuaternaryMatrix we use inline get() and set() methods that bypass the Sage finite field elements. We want this code to push the limits of matroid computation, so we'd hate to incur even a 10% speed loss.

And as Travis said, issues with the documentation etc. can be fixed over time. We intend to use, improve, and expand this code for a long, long time.

--Stefan.

[Stefan]

I'd have to check carefully. There's likely to be some performance loss even when all obvious bottlenecks are accounted for, because for our special classes BinaryMatrix, TernaryMatrix, QuaternaryMatrix we use inline get() and set() methods that bypass the Sage finite field elements. We want this code to push the limits of matroid computation, so we'd hate to incur even a 10% speed loss.

Oh, we also use bitpacking for GF(3) and GF(4), which in Sage is done only for binary matrices currently, if I'm not mistaken.

[Volker Braun]

On Saturday, May 25, 2013 3:51:06 PM UTC+1, Stefan wrote:

There's likely to be some performance loss even when all obvious bottlenecks are accounted for, because for our special classes BinaryMatrix, TernaryMatrix, QuaternaryMatrix we use inline get() and set() methods that bypass the Sage finite field elements.

Only if get/setting matrix elements is actually the bottleneck. But once you are actually doing something with the matrices then e.g. Strassen-Winograd will run circles around the naive O(N^3) multiplication that you have in there.

[Stefan van Zwam]

Multiplication occurs only a handful of times (when computing invariants). The real bottleneck is this method (focus on self._A which is the matrix):

    cdef bint __exchange(self, long x, long y):

        """

        Put element indexed by ``x`` into basis, taking out element ``y``. Assumptions are that this is a valid basis exchange.

        .. NOTE::

            Safe for noncommutative rings.

        """

        cdef long px, py, r

        px = self._prow[x]

        py = self._prow[y]

        piv = self._A.get_unsafe(px, py)

        pivi = piv ** (-1)

        self._A.row_scale(px, pivi)

        self._A.set_unsafe(px, py, pivi + self._one)       # pivoting without column scaling. Add extra so column does not need adjusting

        for r in xrange(self._A.nrows()):            # if A and A' are the matrices before and after pivoting, then

            a = self._A.get_unsafe(r, py)       # ker[I A] equals ker[I A'] except for the labelling of the columns

            if a and r != px:

                self._A.row_add(r, px, -a)

        self._A.set_unsafe(px, py, pivi)

        self._prow[y] = px

        self._prow[x] = py

        BasisExchangeMatroid.__exchange(self, x, y)

--Stefan.

[Volker Braun]

If that is actually the bottleneck then the relevant matrices just need to override add_multiple_of_row_c with a specialized version instead of using the generic one from matrix0.

It seems that renaming Sage's add_multiple_of_row() into row_add() etc. in your code is just making it unnecessarily difficult to switch to Sage matrices.

[Stefan van Zwam]

Hi Volker, and everyone else on sage-devel,

I did some speed tests with an old version of the code that still used Sage matrices (commit 6d9e689 on bitbucket, with a few tweaks to make it run on Sage 5.10.beta4). Each of these tests performs a matroid computation through extensive pivoting.

I ran three versions of the test:

* Sage matrices
* Same algorithm, using a LeanMatrix instance
* Optimized algorithm for this class of matroids

for four different classes of matroids:

* Binary
* Ternary
* Quaternary
* Regular

Without patch 14627 (on Sage 5.9.beta3):

Comparing the speed of various classes.
==========================================================
Binary matroids:
Sage matrix:
5 loops, best of 5: 101 ms per loop
lean_matrix.BinaryMatrix:
5 loops, best of 5: 21.2 ms per loop
BinaryMatroid and BinaryMatrix:
5 loops, best of 5: 3.08 ms per loop
==========================================================
Ternary matroids:
Sage matrix:
5 loops, best of 5: 5.7 s per loop
lean_matrix.TernaryMatrix:
5 loops, best of 5: 1.29 s per loop
TernaryMatroid:
5 loops, best of 5: 7.98 ms per loop
==========================================================
Quaternary matroids:
Sage matrix:
5 loops, best of 5: 186 ms per loop
lean_matrix.QuaternaryMatrix:
5 loops, best of 5: 205 ms per loop
QuaternaryMatroid and QuaternaryMatrix:
5 loops, best of 5: 13.4 ms per loop
==========================================================
Regular matroids:
Sage matrix:
5 loops, best of 5: 753 ms per loop
lean_matrix.IntegerMatrix:
5 loops, best of 5: 49.3 ms per loop
RegularMatroid and IntegerMatrix:
5 loops, best of 5: 19.8 ms per loop

With patch 14627 (on Sage 5.10.beta4):

Comparing the speed of various classes.
==========================================================
Binary matroids:
Sage matrix:
5 loops, best of 5: 91.2 ms per loop
lean_matrix.BinaryMatrix:
5 loops, best of 5: 21.4 ms per loop
BinaryMatroid and BinaryMatrix:
5 loops, best of 5: 3.24 ms per loop
==========================================================
Ternary matroids:
Sage matrix:
5 loops, best of 5: 240 ms per loop
lean_matrix.TernaryMatrix:
5 loops, best of 5: 1.2 s per loop
TernaryMatroid:
5 loops, best of 5: 8.17 ms per loop
==========================================================
Quaternary matroids:
Sage matrix:
5 loops, best of 5: 434 ms per loop
lean_matrix.QuaternaryMatrix:
5 loops, best of 5: 162 ms per loop
QuaternaryMatroid and QuaternaryMatrix:
5 loops, best of 5: 11.3 ms per loop
==========================================================
Regular matroids:
Sage matrix:
5 loops, best of 5: 687 ms per loop
lean_matrix.IntegerMatrix:
5 loops, best of 5: 45.7 ms per loop
RegularMatroid and IntegerMatrix:
5 loops, best of 5: 18.6 ms per loop


Conclusions:

The only speedup obtained from 14627 is for prime fields (and it is HUGE: a factor 20). But our custom class still outperforms Sage's matrices by at least a factor 3 much of the time. With the more specialized routines that exploit the bit packed representation of the matrix, the speedup is even more significant.

It will require significant effort to bring Sage's matrices up to the same speed. For ternary matroids, it won't happen until Martin Albrecht has generalized M4Ri(e) to fields of other characteristic. So I propose to keep the LeanMatrix code for now, opening a ticket titled something like "Bring Sage matrices into the matroid code without performance loss" for a later date.

--Stefan.

P.S. It looks like GF(4)-matrices had a speed regression between 5.9 and 5.10?!

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "sage-matroid" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-matroid...@googlegroups.com.

[Martin Albrecht]

Hi all,

I didn't follow the thread but for me it would be helpful to get a self
contained description of what the performance issue is with current Sage
matrices:

What operations at what dimensions are slow?

> P.S. It looks like GF(4)-matrices had a speed regression between 5.9 and
> 5.10?!

Can you provide an example so I can test?

Cheers,
Martin

--
name: Martin Albrecht
_pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
_otr: 47F43D1A 5D68C36F 468BAEBA 640E8856 D7951CCF
_www: http://martinralbrecht.wordpress.com/
_jab: martinr...@jabber.ccc.de

[Volker Braun]

Overall timings aren't really that useful, without profiling there is no way to tell what part is slow. Can you at least give your own matrix class the same interface as Sage matrices? It should only be a matter of switching from ... import ... as MyMatrix. If they are not interchangeable then nobody will be able to profile it and no progress will be made.

[Volker Braun]

I'm under the following impression, somebody correct me if I'm wrong: Their only real application for matrices are as a container of row vectors. And Sage's add_multiple_of_row is just the generic implementation for most (all?) backends.

[Travis Scrimshaw]

Hey Stefan,

   From just looking at the function, there is only one case in the for loop that is non-trivial, so I'm thinking the function could be restructured as:

    cdef bint __exchange(self, long x, long y):
        """
        Put element indexed by ``x`` into basis, taking out element ``y``. Assumptions are that this is a valid basis exchange.

        .. NOTE::

            Safe for noncommutative rings.
        """
        cdef long px, py,

r
        cdef LeanMatrix A = self._A
        px = self._prow[x]
        py = self._prow[y]
        piv = A.get_unsafe(px, py)

        pivi = piv ** (-1)

        # This would be a simple modification of the previous, removing the for loop
        #a = pivi + self._one
        #A.row_scale(px, pivi)
        #A.set_unsafe(px, py, a)       # pivoting without column scaling. Add extra so column does not need adjusting

        # if A and A' are the matrices before and after pivoting, then

        # ker[I A] equals ker[I A'] except for the labelling of the columns  

        #if a:
        #    A.row_add(px, px, -a)
        #A.set_unsafe(px, py, pivi)

        # This is with some other (possible) (micro) optimizations
        if pivi + self._one != 0:

            A.row_scale(px, -pivi * pivi)

        else:
            A.row_scale(px, pivi)

        A.set_unsafe(px, py, pivi)

        self._prow[y] = px
        self._prow[x] = py
        BasisExchangeMatroid.__exchange(self, x, y)

For the second one, I use the fact that

pivi * A_ij - (pivi + 1) * pivi * A_ij = -pivi^2 * A_ij

All I've done is unraveled the outcome of the code by following the logic (interpret that as I haven't tested it), but this should net (again, provided I didn't make any mistakes) a nice speed increase.

Also, if you're not really using matrices, such as doing a significant number of multiplications/using category model/etc., but instead as a nice data structure, I'd be more inclined to say to strip out the non speed-essential methods, making it as lightweight and specific as possible and converting to proper matrices when needed, and include it into sage (after perhaps changing the name). How does this sound to everyone?

Best,

Travis

[Dima Pasechnik]

On 2013-05-28, Travis Scrimshaw <tsc...@ucdavis.edu> wrote:
> ------=_Part_1805_33045155.1369720622849
> Content-Type: text/plain; charset=ISO-8859-1

it looks like something that should be well-known how to implement
well, as in the simplex methods for linear programming one does
that kind of operations a lot.

Also, did you try implementing this using numpy arrays?
This would give quite a speedup, IMHO.
I imagine there could be a BLAS function (something like daxpy.f) called from numpy,
to do the most expensive part (and BLAS can be much faster than C code on
good CPUs, due to parallelization).

> I'd be more inclined to say to strip out the non
> speed-essential methods, making it as lightweight and specific as possible
> and converting to proper matrices when needed, and include it into sage
> (after perhaps changing the name). How does this sound to everyone?
>
> Best,
> Travis

Dima>

[Volker Braun]

On Tuesday, May 28, 2013 4:52:29 PM UTC+1, Dima Pasechnik wrote:

Also, did you try implementing this using numpy arrays? 

Numpy doesn't support bitpacked Z/n Z so thats of no use.

[Stefan]

Hi Travis,

Thanks for your suggestion! I'll experiment with it later today.

Also, if you're not really using matrices, such as doing a significant number of multiplications/using category model/etc., but instead as a nice data structure, I'd be more inclined to say to strip out the non speed-essential methods, making it as lightweight and specific as possible and converting to proper matrices when needed, and include it into sage (after perhaps changing the name). How does this sound to everyone?

That is, in fact, what we've done. The class LeanMatrix and its subclasses only implement the essential (to us) operations, using cdef methods for speed. They are only used internally: the user inputs Sage matrices, and whenever the code returns a matrix, it's a Sage matrix (some advanced, underscored, methods are an exception, but those are mainly for internal use). 

--Stefan.

[Stefan]

That method is the generic version, which works for arbitrary rings (including non-commutative ones). The BinaryMatroid, TernaryMatroid, QuaternaryMatroid classes have optimized versions already.

--Stefan.