More competitive and comprehensive finite dimensional algebras

[SimonKing]

Issues

• Currently, sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element is Python.

• Multiplication is implemented by

  self.__class__(self.parent(), self._vector * other._matrix)
  

  which for some matrix implementations is not the fastest way:

  sage: v = random_vector(GF(2),2000)
  sage: M = random_matrix(GF(2),2000,2000)
  sage: %timeit v*M
  The slowest run took 8.27 times longer than the fastest. This could mean that an intermediate result is being cached.
  100 loops, best of 3: 4.96 ms per loop
  sage: %timeit M*v
  1000 loops, best of 3: 172 µs per loop
  sage: v = random_vector(GF(3),2000)
  sage: M = random_matrix(GF(3),2000,2000)
  sage: %timeit v*M
  The slowest run took 37.59 times longer than the fastest. This could mean that an intermediate result is being cached.
  1000 loops, best of 3: 1.11 ms per loop
  sage: %timeit M*v
  1000 loops, best of 3: 974 µs per loop
  sage: v = random_vector(QQ,2000)
  sage: M = random_matrix(QQ,2000,2000)
  sage: %timeit v*M
  1 loop, best of 3: 1.51 s per loop
  sage: %timeit M*v
  1 loop, best of 3: 3.32 s per loop
  sage: v = random_vector(GF(9,'x'),2000)
  sage: M = random_matrix(GF(9,'x'),2000,2000)
  sage: %timeit v*M
  1 loop, best of 3: 163 ms per loop
  sage: %timeit M*v
  1 loop, best of 3: 166 ms per loop
  

  (Remark: The timings for GF(9,'x') are with MeatAxe)

• vector*matrix respectively matrix*vector may be a lot slower than matrix*matrix:

  sage: M = random_matrix(GF(9,'x'),2000,2000)
  sage: v = random_matrix(GF(9,'x'),2000,1)
  sage: %timeit M*v
  10 loops, best of 3: 52.4 ms per loop
  sage: v1 = v.transpose()
  sage: %timeit v1*M
  The slowest run took 37.36 times longer than the fastest. This could mean that an intermediate result is being cached.
  10000 loops, best of 3: 14.6 µs per loop
  

  The last example is 163 ms vs. 14.6 µs!!!

• In my applications, I need finite dimensional path algebra quotients. The current implementation corresponds to a finite dimensional path algebra quotient where the underlying quiver has a single vertex.

• Creation of an element also involves construction of its multiplication matrix. From my experience, it can be faster (in some applications) to not construct the matrix of it right away, but be lazy. E.g., when one computes Gröbner bases, only multiplication by power products of the algebra generators are needed, and then it is faster to map the defining vector of the element repeatedly by multiplying with the matrices of the generators, instead of doing matrix by matrix multiplications.

Suggestion

1. Re-implement it in Cython
2. Replace vector*matrix by either row*matrix or matrix*column, being flexible enough that both implementations are supported (because it depends on the field whether row*matrix or matrix*column is fastest).
3. Allow using a quiver with more than one vertices.
4. Be more lazy, i.e. do costly computations such as the construction of the multiplication matrix only when needed.

[jdemeyer]

FYI: PARI/GP has recently gained functionality to deal with general finite dimensional algebras. But it seems that this ticket is mostly about basic arithmetic, so I'm not sure that it's relevant.

[SimonKing]

Replying to jdemeyer:

FYI: PARI/GP has recently gained functionality to deal with general finite dimensional algebras. But it seems that this ticket is mostly about basic arithmetic, so I'm not sure that it's relevant.

This ticket is not necessarily about basic arithmetic. Eventually, I finally want to implement my F5 algorithm for modules over path algebra quotients (after being distracted by other things for several years), starting with finite dimensional quotients.

So, having an implementation of finite dimensional algebras over any field (but for me most importantly over finite not necessarily prime fields!) that is *really* competitive would be nice.

Can you give me a pointer for the PARI/GP implementation? Would it be difficult to wrap in cython?

[jdemeyer]

Replying to SimonKing:

So, having an implementation of finite dimensional algebras over any field (but for me most importantly over finite not necessarily prime fields!) that is *really* competitive would be nice.

An algebra over F_q is just a special case of an algebra over F_p. PARI/GP supports algebras over Q and over F_p.

Can you give me a pointer for the PARI/GP implementation?

See ​https://pari.math.u-bordeaux.fr/dochtml/html-stable/ section "Associative and central simple algebras".

Would it be difficult to wrap in cython?

I don't think so. Depending on your needs, you could also use cypari2.

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

So, having an implementation of finite dimensional algebras over any field (but for me most importantly over finite not necessarily prime fields!) that is *really* competitive would be nice.

An algebra over F_q is just a special case of an algebra over F_p.

I don't buy that. A matrix over F_p should be the same as an algebra over F_q, but the current default in sage differs totally in performance.

Anyway, if you wanted to say that PARI/GP is good in both cases then it's fine.

See ​https://pari.math.u-bordeaux.fr/dochtml/html-stable/ section "Associative and central simple algebras".

Thank you!

Would it be difficult to wrap in cython?

I don't think so. Depending on your needs, you could also use cypari2.

What's that?

[jdemeyer]

cypari2 is a new package with what used to be sage.libs.pari. It is a Cython interface to PARI. It wraps PARI objects in Python objects. Using this is slower than directly calling the PARI library. On the other hand, if your computations spend most of the time in the PARI library, it wouldn't matter for performance.

[SimonKing]

Replying to jdemeyer:

cypari2 is a new package with what used to be sage.libs.pari. It is a Cython interface to PARI. It wraps PARI objects in Python objects. Using this is slower than directly calling the PARI library. On the other hand, if your computations spend most of the time in the PARI library, it wouldn't matter for performance.

Thanks. Are algebras wrapped? If found cypari2 by googling, but a quick look seems to indicate that the conversions are about numbers, not algebras.

[SimonKing]

How can I do the following pari code in Sage?

  { m_i = [0,-1,0, 0;
           1, 0,0, 0;
           0, 0,0,-1;
           0, 0,1, 0];
    m_j = [0, 0,-1,0;
           0, 0, 0,1;
           1, 0, 0,0;
           0,-1, 0,0];
    m_k = [0, 0, 0, 0;
           0, 0,-1, 0;
           0, 1, 0, 0;
           1, 0, 0,-1];
    A = alginit(nfinit(y), [matid(4), m_i,m_j,m_k],  0); }

Defining the matrices works:

sage: Mi = pari("""
....: [0,-1,0, 0;
....:            1, 0,0, 0;
....:            0, 0,0,-1;
....:            0, 0,1, 0]""")
sage: Mj = pari("""
....: [0, 0,-1,0;
....:            0, 0, 0,1;
....:            1, 0, 0,0;
....:            0,-1, 0,0];""")
sage: Mk = pari("""
....: [0, 0, 0, 0;
....:            0, 0,-1, 0;
....:            0, 1, 0, 0;
....:            1, 0, 0,-1];""")
sage: Mj
[0, 0, -1, 0; 0, 0, 0, 1; 1, 0, 0, 0; 0, -1, 0, 0]
sage: Mi*Mk
[0, 0, 1, 0; 0, 0, 0, 0; -1, 0, 0, 1; 0, 1, 0, 0]

However, my attempts to define the algebra failed:

sage: pari("nfinit(y)").alginit([pari("matid(4)"), Mi,Mj,Mk], 0)
...
PariError: incorrect priority in nffactor: variable 0 >= y
sage: pari.alginit(pari("nfinit(y)"), [pari("matid(4)"), Mi,Mj,Mk], 0)
...
AttributeError: 'cypari2.pari_instance.Pari' object has no attribute 'alginit'

[jdemeyer]

Replying to SimonKing:

Are algebras wrapped?

Everything is wrapped. That's the nice thing about cypari2: it automagically wraps all GP functions (assuming that the argument types are understood; PARI functions taking functions as input are not supported for example).

It's probably not so relevant, but it's best to use #23518.

a quick look seems to indicate that the conversions are about numbers, not algebras.

Well yes, it's a Python package so it only supports a limited set of conversions for standard Python types. Sage does support more conversions though. In any case, don't be blinded by the conversions, that's really just a detail. If PARI supports your algorithms, surely the conversion between Sage and PARI should be possible to implement.

[jdemeyer]

Use a different variable name :-)

sage: pari("nfinit(t)").alginit(["matid(4)", Mi, Mj, Mk], 0)

Note that you don't need pari() in the arguments of the method: the inputs are automatically converted to PARI. This means that strings are evaluated.

With #23518, you should be able to write this as

sage: pari.alginit("nfinit(t)", ["matid(4)", Mi, Mj, Mk], 0)

or

sage: pari.alginit(pari.nfinit("t"), [pari.matid(4), Mi, Mj, Mk], 0)

[SimonKing]

Replying to jdemeyer:

Use a different variable name :-)

sage: pari("nfinit(t)").alginit(["matid(4)", Mi, Mj, Mk], 0)

Why is that? pari("nfinit(y)") by itself does work.

[SimonKing]

What's wrong with this?

sage: M1 = random_matrix(GF(2),1024,1024)
sage: M2 = random_matrix(GF(2),1024,1024)
sage: A = pari(GF(2)).alginit(["matid(1024)", M1, M2],0)
...
PariError: incorrect type in alginit (t_POL)

So, does the first argument of alginit really needs to be a number field, not a finite field?

[SimonKing]

Reading that over finite fields one uses algtableinit, I tried

sage: MT = pari(["matid(1024)", M1,M2])
sage: A = MT.algtableinit(2)
...
PariError: incorrect type in algtableinit (t_VEC)

which didn't work either. Why?

[SimonKing]

The following doesn't look promising:

sage: pM1 = pari(M1)
sage: pM2 = pari(M2)
sage: %timeit M1*M2
1000 loops, best of 3: 1.08 ms per loop
sage: %timeit pM1*pM2
1 loop, best of 3: 201 ms per loop

[jdemeyer]

Replying to SimonKing:

Reading that over finite fields one uses algtableinit, I tried

sage: MT = pari(["matid(1024)", M1,M2])
sage: A = MT.algtableinit(2)
...
PariError: incorrect type in algtableinit (t_VEC)

which didn't work either. Why?

Don't you need a list of 1024 matrices in dimension 1024?

[SimonKing]

Replying to jdemeyer:

Don't you need a list of 1024 matrices in dimension 1024?

Ouch, yes! Apparently I thought of constructing the algebra that is generated by these matrices (which I guess has dimension much larger than 1024).

[SimonKing]

Replying to jdemeyer:

Replying to SimonKing:

Reading that over finite fields one uses algtableinit, I tried

sage: MT = pari(["matid(1024)", M1,M2])
sage: A = MT.algtableinit(2)
...
PariError: incorrect type in algtableinit (t_VEC)

which didn't work either. Why?

Don't you need a list of 1024 matrices in dimension 1024?

But in any case, Pari complains about type (t_VEC), not about value. So, the number of matrices isn't the only problem.

[jdemeyer]

Replying to SimonKing:

But in any case, Pari complains about type (t_VEC)

Don't take the PARI error messages too literally...

[SimonKing]

I think I shouldn't be overambitious in this ticket: The finite-dimensional path algebra quotient stuff should be for a later ticket, where I will probably sub-class the now cythonised finite dimensional algebra elements.

Changes:

• The element class is now in Cython.
• The multiplication matrix of an element is now lazy, which means that initialisation of an element will be substantially faster.
• An element is encoded by a matrix with a single row, not by a vector. As I have demonstrated, multiplying a single row matrix by a square matrix is always faster than multiplying a vector by a square matrix.

Caveat: It seems to me that OVER QQ multiplying a square matrix with a single-column matrix is faster than multiplying a single-row matrix with a square matrix. But that's the only case where I found that behaviour. Since I am not interested in the rational case anyway, I don't plan to provide a special implementation over QQ where elements are encoded as single-column matrices rather than single-row matrices.

PS: For me, tests pass...

---

New commits:

​a34ee26

Faster implementation of finite dimensional algebras

[SimonKing]

Back to pari:

sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1, 0], [0, 1]]), Matrix([[0, 1], [0, 0]])])
sage: pari([x.matrix() for x in A.gens()]).algtabaleinit()
Traceback (most recent call last):
...
PariError: incorrect type in algtableinit (t_VEC)

What's wrong?

[SimonKing]

Replying to SimonKing:

What's wrong?

I don't know why, but apparently Pari doesn't accept some associative algebras that are perfectly fine in Sage, while it accept others:

sage: B = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,1]]), Matrix([[0,0,0], [1,0,1], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [1,0,1]])])
sage: pari([x.matrix() for x in B.gens()]).algtableinit()
[0, 0, 0, 0, 0, 0, [1, 0, 0; 0, 1, 0; 0, 0, 1], [1, 0, 0; 0, 1, 0; 0, 0, 1], [[1, 0, 0; 0, 1, 0; 0, 0, 1], [0, 0, 0; 1, 0, 1; 0, 0, 0], [0, 0, 0; 0, 0, 0; 1, 0, 1]], 0, [3, 0, 1]]
sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: pari([x.matrix() for x in A.gens()]).algtableinit()
---------------------------------------------------------------------------
PariError                                 Traceback (most recent call last)
<ipython-input-31-0e66e23e7476> in <module>()
----> 1 pari([x.matrix() for x in A.gens()]).algtableinit()

cypari2/auto_gen.pxi in cypari2.gen.Gen_auto.algtableinit()

cypari2/handle_error.pyx in cypari2.handle_error._pari_err_handle()

PariError: incorrect type in algtableinit (t_VEC)

Is there a clear reason or are associative algebras in Pari just broken?

[SimonKing]

Strangely, the patchbot doesn't test it (although it needs review). I'm trying to ?kick it (I don't know if that's still needed nowadays, I learned about ?kick a couple of years ago).

[chapoton]

you need to fill the Author field, otherwise it is deemed unsafe

[SimonKing]

Replying to chapoton:

you need to fill the Author field, otherwise it is deemed unsafe

Oops, sorry! And thank you for filling it.

[SimonKing]

The patchbot reports a failure in sage.geometry.polyhedron.backend_normaliz.Polyhedron_normaliz.integral_hull. Is that related? In any case, it seems that it just changes because of a different sorting of the result (and if it does involve finite dimensional algebras, the sorting of the output would change, because sorting is now defined and is compatible with sorting of matrices).

[chapoton]

#23586

[SimonKing]

I notice a problem: Old pickles won't unpickle.

When I save a finite dimensional algebra element in the old version, and load it in the new version, I get

sage: x = load('backup.sobj')
sage: x
<repr(<sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebra_with_category.element_class at 0x7f836819c0b8>) failed: AttributeError: 'NoneType' object has no attribute 'list'>
sage: x._vector is None
True
sage: x.__dict__

{'_default_filename': '/home/king/Sage/git/sage/backup.sobj',
 '_matrix': [0 1 0]
 [0 0 0]
 [0 0 0],
 '_vector': (0, 1, 0)}
sage: type(x)
<class 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebra_with_category.element_class'>
sage: type(x) is A.element_class
False
sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x = A.1
sage: type(x)
<type 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebraElement'>

So, the old pickle gives the a python class that is derived from the new cython class and is not the same as the new element_class.

How can this be cured? There is unpickle_override, but I don't know if this helps in the current situation.

[SimonKing]

Would it help to implement a __setstate__ method that handles assignment of the old data to the new cdef public slots?

[tscrim]

IIRC, what I've done in the past when changing a Python class to a cdef class is to use register_unpickle_override. I think something else you can do is just create an essentially dummy Python class that does nothing other than be a subclass of the Cython class.

Also, something else to be aware of is checking you are inheriting things from the category. IIRC, when Cythonizing CombinatorialFreeModuleElement, Nicolas did some special workup to get category inheritance to work properly (my solution was the dummy Python class).

[tscrim]

Also, I get similar behavior as mentioned in the ticket description for Z and Q:

sage: v = random_vector(R, 2000)
sage: M = random_matrix(R, 2000)
sage: %timeit M * v
sage: %timeit v * M
sage: %timeit matrix(v) * M
sage: %timeit M * matrix(v).transpose()
sage: vp = matrix(v)
sage: %timeit vp * M
sage: vp = vp.transpose()
sage: %timeit M * vp

yields for R = QQ:

1 loop, best of 3: 3.09 s per loop
1 loop, best of 3: 1.49 s per loop
1 loop, best of 3: 331 ms per loop
1 loop, best of 3: 210 ms per loop
1 loop, best of 3: 330 ms per loop
1 loop, best of 3: 211 ms per loop

and R = ZZ:

1 loop, best of 3: 657 ms per loop
10 loops, best of 3: 116 ms per loop
10 loops, best of 3: 60.1 ms per loop
10 loops, best of 3: 32.8 ms per loop
10 loops, best of 3: 59.6 ms per loop
10 loops, best of 3: 32.3 ms per loop

Interestingly for R = ZZ['x'] (with averaging 10 trials of 700x700 matrices):

1 loop, best of 3: 252 ms per loop
1 loop, best of 3: 197 ms per loop
1 loop, best of 3: 257 ms per loop
1 loop, best of 3: 210 ms per loop
1 loop, best of 3: 255 ms per loop
1 loop, best of 3: 212 ms per loop

Basically the multiplication of (dense) matrix * vector (and vice versa) is using the naïve algorithm and does not take advantage of specialized implementations:

        cdef Py_ssize_t i
        return sum([self.column(i, from_list=True) * v[i]
                    for i in xrange(self._ncols)], M(0))

While _matrix_times_vector_ clearly deserves specialized implementations (same for _vector_times_matrix_), it will take some more specialized code to improve it. The amazing thing is Q has a specialized implementation for _vector_times_matrix_, but it is not as fast as going back and forth with a matrix (flint is doing something really good)! O_o That would be a project for a separate ticket.

[SimonKing]

__setstate__ will work.

Question: When unpickling, the dictionary contains an item _default_filename, which provides the file name of the pickle. Should I simply ignore it? I think so...

I agree that the vector vs. matrix problem should be treated separately. I think currently it is best to do "single row matrix" times "square matrix".

[SimonKing]

PS: Apart from _default_filename, should other unknown items simply be ignored? Or better raise a TypeError? when an unexpected property occurs at unpickling?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3c39dac

Unpickling old f.d. algebra elements

[SimonKing]

The test for the new __setstate__ method is indirect, but I tested that indeed a pickle created with the old version correctly unpickles in the new version (but: It results in a python class derived from FiniteDimensionalAlgebraElement, rather than in the cython class FiniteDimensionalAlgebraElement itself).

While I was at it, I introduced an unpickle function. Reasons: It is slightly faster than the previous way, and if the multiplication matrix has previously been computed then the matrix will be taken into account in the pickle (thus, avoiding a re-computation).

[SimonKing]

PS: __setstate__ should in real applications only be invoked when unpickling old pickles. In that case, the resulting element will be a python class, thus, will have a __dict__. Then, the special cdef slots are filled and all remaining attributes will be put into __dict__.

[SimonKing]

According to the title, this ticket is about "more competitive" implementation. So, here are some timings. Old implementation:

sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 172.77 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 81.9 µs per loop
sage: %timeit y+z
The slowest run took 8.10 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 66.4 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 15.88 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 348 µs per loop

New implementation:

sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 380.33 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 9.52 µs per loop
sage: %timeit y+z
The slowest run took 17.05 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 5.64 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 30.65 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 261 µs per loop

So, I kept my promise...

[SimonKing]

Let's consider different fields (note that I am using meataxe, for GF(9,'x')).

Old:

sage: A = FiniteDimensionalAlgebra(GF(2), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 95.24 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 96.3 µs per loop
sage: %timeit y+z
10000 loops, best of 3: 84.7 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 10.92 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 637 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 153.07 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 77 µs per loop
sage: %timeit y+z
The slowest run took 5.54 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 64.2 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 26.46 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 375 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(9,'x'), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 201.31 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 82.8 µs per loop
sage: %timeit y+z
The slowest run took 6.70 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 65.2 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 13.58 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 513 µs per loop

New:

sage: A = FiniteDimensionalAlgebra(GF(2), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 162.81 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 17.9 µs per loop
sage: %timeit y+z
The slowest run took 19.10 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.13 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 21.34 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 553 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 183.38 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 13 µs per loop
sage: %timeit y+z
The slowest run took 17.76 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 7.99 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 7.77 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 270 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(9,'x'), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 150.55 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 19.3 µs per loop
sage: %timeit y+z
The slowest run took 16.32 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.92 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 25.80 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 434 µs per loop

So, it is consistently faster, apparently for all fields.

[tscrim]

Great, thank you for the timings.

Some things that will need to be done before positive review:

• You should have a pxd file with the appropriate declarations.
• I would get rid of the is_vector and just replace it with the appropriate isinstance test for clarity (and the compiler does not always inline inline functions).
• I would not mark unpickle_FiniteDimensionalAlgebraElement as inline as it probably never will be inlined.
• In this part of the initialization

  if A == elt.parent():
      self._vector = elt._vector.base_extend(k)
      self.__matrix = elt._matrix.base_extend(k)
  

  I think you should not compute the matrix for elt, but instead keep it lazy:

  if A == elt.parent():
     mat = (<FiniteDimensionalAlgebraElement> elt).__matrix
     if mat is None:
         self.__matrix = None
     else:
         self.__matrix = mat.base_extend(k)
  

• You should generally avoid the matrix() constructor as it is slower than creating the corresponding MatrixSpace and directly constructing the element from that (at least last time I checked). This way avoids a lot of the logic for parsing arguments that is there in matrix(). So you could do changes like:

  -self._vector = matrix(self._parent.base_ring(), 1,len(v), list(v))
  +self._vector = MatrixSpace(self._parent.base_ring(), len(v)).one()
  

• You should add a typecast here for speed:

  -cdef FiniteDimensionalAlgebraElement a = ~self
  +cdef FiniteDimensionalAlgebraElement a = <FiniteDimensionalAlgebraElement> (~self)
  

  Also, on the line below you know that a.__matrix is set (if invertible), so you can avoid the extra Python call.
• Add a typecase for other in _mul_.
• In _matrix, you should add a cdef int i (Cython makes better code with this). Actually, better to do

  cdef int i
  cdef tuple table = <tuple> A.table()
  ret = sum(self._vector[0,i] * table[i] for i in xrange(A.degree()))
  self.__matrix = matrix(A.base_ring(), ret)
  

• As Jeroen said on sage-devel, by having a __dict__, the class is slower than without it. I suspect there are probably better ways to deal with the unpickling than to have a __dict__. I can look into this if you want.

[SimonKing]

Replying to tscrim:

• As Jeroen said on sage-devel, by having a __dict__, the class is slower than without it. I suspect there are probably better ways to deal with the unpickling than to have a __dict__. I can look into this if you want.

You misunderstand what I was doing.

We are talking here about the __setstate__ method, that will only be called when someone unpickles an old pickle of finite-dimensional algebra elements. Pickling of the new implementation relies on __reduce__ and an unpickle function that has nothing to do with __setstate__.

In the old implementation, an element was an instance of an element_class derived from a Python class. In other words, it is a dynamic class whose base is the former Python class FiniteDimensionalAlgebraElement. Hence, when unpickling it, the class will still be a dynamic (Python) class whose base is the (now!) Cython class FiniteDimensionalAlgebraElement.

Therefore, when unpickling an old pickle, we can be sure that the class has a __dict__, simply because it is a Python class! Nonetheless, just to be on the safe side, I check whether it really has a __dict__, and only if it has, I use it to store additional attributes that the user might have assigned to the element before pickling it.

In other words, FiniteDimensionalAlgebraElement has no __dict__, and we are merely talking about an attempt to recover as gracefully as possible in the unlikely event that the user has an old pickle containing additional attributes. But if that situation occurs, the class will be a derived Python class that has a __dict__ -- I am not adding it.

[SimonKing]

Replying to tscrim:

• You should generally avoid the matrix() constructor as it is slower than creating the corresponding MatrixSpace and directly constructing the element from that (at least last time I checked).

I can confirm, although the difference isn't huge. Non-square matrix is faster with MatrixSpace:

sage: %timeit M = matrix(GF(7), 1,1000, range(1000))
The slowest run took 5.83 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 66.3 µs per loop
sage: %timeit M = MatrixSpace(GF(7), 1,1000)(range(1000))
The slowest run took 5.47 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 56.4 µs per loop

Square unit matrix (mind that we can not take .one(), as we want a mutable copy):

sage: %timeit M = matrix(GF(7), 1000,1000, 1)
1000 loops, best of 3: 738 µs per loop
sage: %timeit M = MatrixSpace(GF(7), 1000)(1)
The slowest run took 5.07 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 712 µs per loop

[tscrim]

Replying to SimonKing:

I see. I am not completely convinced that you want to rebuild the old class rather than the new class. The reason you may unpickle to a different parent is because of the UniqueFactory of GF(p) having the Sage version as part of the key. If you do it over, e.g., Z, then it ends up being two "different" element classes of the same parent. Now as soon as you do an operation, it goes to the new class. However, I am still willing to set a positive review since this is close enough to a pathological example and everything seems to work.

I am encountering a difference because of extension classes versus regular classes:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: A.an_element()[0]  # TypeError - defined in the category and not inherited

So no matter what, you currently need a Python class for the element_class.

One way you can add a test is take an old pickle string and then load that (which is how I am testing things).

[tscrim]

Replying to SimonKing:

Replying to tscrim:

• You should generally avoid the matrix() constructor as it is slower than creating the corresponding MatrixSpace and directly constructing the element from that (at least last time I checked).

Square unit matrix (mind that we can not take .one(), as we want a mutable copy):

I don't see where you need a mutable copy. Can you show me where that is? Also, one() does some extra checks, but the fastest way is identity_matrix() (all MS(1) does is call MS.identity_matrix().__copy__()).

[SimonKing]

Replying to tscrim:

I don't see where you need a mutable copy.

I thought at some point of the code matrix entries are assigned to after creation, but it seems that it is in fact not the case. In that case, one should .zero() and .one(), as it is cached.

[SimonKing]

Replying to tscrim:

Replying to SimonKing:

I see. I am not completely convinced that you want to rebuild the old class rather than the new class.

And I don't see how I could change that. I'd be happy to unpickle the old dynamic class that is derived from the old Python class FiniteDimensionalAlgebraElement just as plain FiniteDimensionalAlgebraElement. There would be no problem if the old pickle would state that the class is the element class of the algebra (because that in fact is FiniteDimensionalAlgebraELement, not just a derived class).

However, the old pickle apparently states that the class is a dynamic class obtained from a certain class, and does not state that it is the element class of a specific parent (as it should!). Indeed, in the old version:

sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x = A.1
sage: type(A.1)._reduction

(<function dynamic_class at 0x7f6a647bbc08>,
 ('FiniteDimensionalAlgebra_with_category.element_class',
  (<class 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebraElement'>,
   <class 'sage.categories.category.JoinCategory.element_class'>),
  None,
  None,
  None))

So, how could I change that without changing the old pickle?

The reason you may unpickle to a different parent...

Am I?

I am encountering a difference because of extension classes versus regular classes:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: A.an_element()[0]  # TypeError - defined in the category and not inherited

That's totally odd. By how the category framework works, all methods defined in the category should still be available for an element class that is a cython class.

Is it not the case that the element class of a parent P coincides with P.Element if P.Element is a cython class (which is no the case for finite dimensional algebras, but used not to be the case in the old implementation).

So no matter what, you currently need a Python class for the element_class.

No. It is perfectly normal for an element class to be Cython.

[SimonKing]

Replying to tscrim:

• You should generally avoid the matrix() constructor as it is slower than creating the corresponding MatrixSpace and directly constructing the element from that (at least last time I checked).

For reference:

sage: %timeit M = MatrixSpace(QQ,1,5)()
The slowest run took 15.03 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 5.06 µs per loop
sage: %timeit M = matrix.zero(QQ,1,5)
The slowest run took 13.26 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 7.55 µs per loop
sage: %timeit M = MatrixSpace(QQ,1,5).zero()
The slowest run took 10.71 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 12.3 µs per loop
sage: %timeit M = matrix(QQ,1,5,0)
The slowest run took 15.02 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 13.2 µs per loop
sage: %timeit M = MatrixSpace(QQ,5)(1)
The slowest run took 191.83 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.88 µs per loop
sage: %timeit M = matrix.identity(QQ,5)
The slowest run took 15.86 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 7.5 µs per loop
sage: %timeit M = MatrixSpace(QQ,5).one()
The slowest run took 9.16 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 18.8 µs per loop
sage: %timeit M = matrix(QQ,5,5,1)
The slowest run took 9.29 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 14.8 µs per loop

That's odd. .zero() and .one() are cached methods. Why are they SLOWER than constructing a matrix from scratch?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b5d4357

Further performance tweaks in f.d. algebras

[SimonKing]

In the new (final?) commit, I am addressing your objections as follows.

Replying to tscrim:

Great, thank you for the timings.

New timings:

sage: A = FiniteDimensionalAlgebra(GF(2), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 138.54 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 16.6 µs per loop
sage: %timeit y+z
The slowest run took 16.78 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 5.97 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 9.31 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 509 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 2667.19 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 12 µs per loop
sage: %timeit y+z
The slowest run took 18.26 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.51 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 6.92 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 252 µs per loop
sage: A = FiniteDimensionalAlgebra(GF(9,'x'), [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 165.25 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 19 µs per loop
sage: %timeit y+z
The slowest run took 8.29 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.15 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 8.84 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 402 µs per loop
sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x,y,z = A.gens()
sage: %timeit x*y
The slowest run took 165.79 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 8.2 µs per loop
sage: %timeit y+z
The slowest run took 15.54 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6 µs per loop
sage: %timeit loads(dumps(x))
The slowest run took 26.51 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 251 µs per loop                                                                                                                                                          

=> Thank you for your hints on performance! It resulted in a noticeable improvement.

• You should have a pxd file with the appropriate declarations.

Done.

• I would get rid of the is_vector and just replace it with the appropriate isinstance test for clarity (and the compiler does not always inline inline functions).

Done.

• I would not mark unpickle_FiniteDimensionalAlgebraElement as inline as it probably never will be inlined.

Done.

• In this part of the initialization

  if A == elt.parent():
      self._vector = elt._vector.base_extend(k)
      self.__matrix = elt._matrix.base_extend(k)
  

  I think you should not compute the matrix for elt, but instead keep it lazy:

Done.

• You should generally avoid the matrix() constructor as it is slower than creating the corresponding MatrixSpace and directly constructing the element from that (at least last time I checked).

Done.

Also, on the line below you know that a.__matrix is set (if invertible), so you can avoid the extra Python call.

Done.

• Add a typecase for other in _mul_.

Done.

• In _matrix, you should add a cdef int i (Cython makes better code with this). Actually, better to do

  cdef int i
  cdef tuple table = <tuple> A.table()
  ret = sum(self._vector[0,i] * table[i] for i in xrange(A.degree()))
  self.__matrix = matrix(A.base_ring(), ret)
  

Done.

• As Jeroen said on sage-devel, by having a __dict__, the class is slower than without it. I suspect there are probably better ways to deal with the unpickling than to have a __dict__.

I didn't, see discussion above.

As pointed out on sage-devel, P.element_class should be pickled by construction and not by constituents. That's a bug in the category framework that will generally be a problem when a former python element class becomes a cython element class. So, that's for a different ticket.

Something else: All tests pass :-)

[tscrim]

Replying to SimonKing:

Replying to tscrim:

Replying to SimonKing:

I see. I am not completely convinced that you want to rebuild the old class rather than the new class.

And I don't see how I could change that. I'd be happy to unpickle the old dynamic class that is derived from the old Python class FiniteDimensionalAlgebraElement just as plain FiniteDimensionalAlgebraElement. There would be no problem if the old pickle would state that the class is the element class of the algebra (because that in fact is FiniteDimensionalAlgebraELement, not just a derived class).

Ugg...the pickle is much more fugly than I was thinking. So it will be what it will be. At least, I don't see a way we can get around it either.

So, how could I change that without changing the old pickle?

The reason you may unpickle to a different parent...

Am I?

You specifically, no. IMO, it is a shortcoming of UniqueFactory (which is used for the GF(p) constructor).

I am encountering a difference because of extension classes versus regular classes:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: A.an_element()[0]  # TypeError - defined in the category and not inherited

That's totally odd. By how the category framework works, all methods defined in the category should still be available for an element class that is a cython class.

However, these are special Python methods and behave different than normal methods. I ran into this same problem on #22632 (I say _mul_, but it really comes from the lack of the __mul__ IIRC), and this was the only way I could work around it. I'm not sure if #23435 would actually be a proper fix for this or not.

Is it not the case that the element class of a parent P coincides with P.Element if P.Element is a cython class (which is no the case for finite dimensional algebras, but used not to be the case in the old implementation).

So no matter what, you currently need a Python class for the element_class.

No. It is perfectly normal for an element class to be Cython.

However, it is the most straightforward way to get special methods implemented in the category to work while still retaining the Cythonization.

[tscrim]

It is a backwards incompatible change in output of vector() returning a matrix instead of a vector. You should change vector() to return self._vector[0] (and maybe add a quick comment there to say self._vector is now a 1xn matrix).

[tscrim]

Replying to SimonKing:

That's odd. .zero() and .one() are cached methods. Why are they SLOWER than constructing a matrix from scratch?

You are not really constructing a matrix from scratch but instead doing a copy. However, what is really strange is this:

sage: MS = MatrixSpace(QQ, 2)
sage: %timeit MS.one()
The slowest run took 9.92 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 5.36 µs per loop
sage: %timeit MS.identity_matrix()
The slowest run took 47.13 times longer than the fastest. This could mean that an intermediate result is being cached.
10000000 loops, best of 3: 60.7 ns per loop
sage: MS.one == MS.identity_matrix
True
sage: MS.one is MS.identity_matrix
True

I guess that is coming from the alias and @cached_method not as well-behaved with aliases. However, that is a huge difference (and something I will need to remember).

[jdemeyer]

Replying to tscrim:

I am encountering a difference because of extension classes versus regular classes:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: A.an_element()[0]  # TypeError - defined in the category and not inherited

General request: when posting this kind of stuff on Trac, please do include the tracebacks. Right now, in order to follow this discussion, I would have to compile this branch just to see what kind of error you are talking about. Not fun!

[jdemeyer]

Replying to tscrim:

I'm not sure if #23435 would actually be a proper fix for this or not.

I don't think that #23435 would really fix anything regarding this.

[tscrim]

Replying to jdemeyer:

Replying to tscrim:

I am encountering a difference because of extension classes versus regular classes:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: A.an_element()[0]  # TypeError - defined in the category and not inherited

General request: when posting this kind of stuff on Trac, please do include the tracebacks. Right now, in order to follow this discussion, I would have to compile this branch just to see what kind of error you are talking about. Not fun!

There is no (functional) traceback:

sage: A.an_element()[0]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-26af788a2600> in <module>()
----> 1 A.an_element()[Integer(0)]

TypeError: 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebraElement' object has no attribute '__getitem__'

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0a88fb4

Avoid incompatible change in FiniteDimensionalAlgebraElement.vector()

[jdemeyer]

Replying to tscrim:

There is no (functional) traceback:

Well, the error message makes sense. Now I know that it's about __getitem__.

[SimonKing]

The latest commit reverts the backwards incompatible change in .vector().

[jdemeyer]

Replying to tscrim:

You are not really constructing a matrix from scratch but instead doing a copy. However, what is really strange is this:

The difference is this:

sage: MS = MatrixSpace(QQ, 2)
sage: MS.__dict__

{'_MatrixSpace__is_sparse': False,
 '_MatrixSpace__matrix_class': <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'>,
 '_MatrixSpace__ncols': 2,
 '_MatrixSpace__nrows': 2,
 '_implementation': 'flint',
 '_reduction': (<class 'sage.matrix.matrix_space.MatrixSpace'>,
  (Rational Field, 2, 2, False, 'flint'),
  {}),
 'identity_matrix': Cached version of <function identity_matrix at 0x7f1014913b90>}

identity_matrix is an attribute on the object MS, while one is an attribute on type(MS).

I recently created getattr_debug which is really useful in situations like this:

sage: MS = MatrixSpace(QQ, 2)
sage: getattr_debug(MS, "one")
getattr_debug(obj=Full MatrixSpace of 2 by 2 dense matrices over ~~~, name='one'):
  type(obj) = <class 'sage.matrix.matrix_space.MatrixSpace_with_category'>
  object has __dict__ slot (<type 'dict'>)
  found 'one' in dict of <class 'sage.matrix.matrix_space.MatrixSpace'>
  got <sage.misc.cachefunc.CachedMethod object at 0x7~~~ (<type 'sage.misc.cachefunc.CachedMethod'>)
  attribute is ordinary descriptor (has __get__)
  object __dict__ does not have 'one'
  calling __get__()
  returning Cached version of <function identity_matrix at ~~~ (<type 'sage.misc.cachefunc.CachedMethodCallerNoArgs'>)
Cached version of <function identity_matrix at 0x7f1014913b90>
sage: getattr_debug(MS, "identity_matrix")
getattr_debug(obj=Full MatrixSpace of 2 by 2 dense matrices over ~~~, name='identity_matrix'):
  type(obj) = <class 'sage.matrix.matrix_space.MatrixSpace_with_category'>
  object has __dict__ slot (<type 'dict'>)
  found 'identity_matrix' in dict of <class 'sage.matrix.matrix_space.MatrixSpace'>
  got <sage.misc.cachefunc.CachedMethod object at 0x7~~~ (<type 'sage.misc.cachefunc.CachedMethod'>)
  attribute is ordinary descriptor (has __get__)
  found 'identity_matrix' in object __dict__
  returning Cached version of <function identity_matrix at ~~~ (<type 'sage.misc.cachefunc.CachedMethodCallerNoArgs'>)
Cached version of <function identity_matrix at 0x7f1014913b90>

[tscrim]

Replying to SimonKing:

The latest commit reverts the backwards incompatible change in .vector().

Thanks. I think because it is an implementation detail, it should be a proper code comment rather than a docstring note.

The only issue now is the __getitem__.

[SimonKing]

Replying to tscrim:

Replying to SimonKing:

The latest commit reverts the backwards incompatible change in .vector().

Thanks. I think because it is an implementation detail, it should be a proper code comment rather than a docstring note.

Shall I change it?

The only issue now is the __getitem__.

OK, I'll provide a __getitem__ that mimics what used to be inherited from the category.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8afaa05

FiniteDimensionalAlgebraElement.__getitem__ replacing the category inherited method

[SimonKing]

Replying to SimonKing:

Replying to tscrim:

Thanks. I think because it is an implementation detail, it should be a proper code comment rather than a docstring note.

Shall I change it?

Done anyway...

The only issue now is the __getitem__.

OK, I'll provide a __getitem__ that mimics what used to be inherited from the category.

Done.

[tscrim]

Replying to SimonKing:

Replying to SimonKing:

Replying to tscrim:

Thanks. I think because it is an implementation detail, it should be a proper code comment rather than a docstring note.

Shall I change it?

Done anyway...

Thanks. I think it is better to hide such details from the front-end user, but it might be useful for someone reading the code.

The only issue now is the __getitem__.

OK, I'll provide a __getitem__ that mimics what used to be inherited from the category.

Done.

It is a more systemic problem:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: len(A.an_element())
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-3f0d37d36d4e> in <module>()
----> 1 len(A.an_element())

TypeError: object of type 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebraElement' has no len()

Granted, I think this might be the only other special method defined in the category.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​2fd7595

Add FiniteDimensionalAlgebraElement.__len__

[SimonKing]

Replying to tscrim:

It is a more systemic problem:

sage: A = FiniteDimensionalAlgebra(ZZ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
sage: len(A.an_element())
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-3f0d37d36d4e> in <module>()
----> 1 len(A.an_element())

TypeError: object of type 'sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element.FiniteDimensionalAlgebraElement' has no len()

Granted, I think this might be the only other special method defined in the category.

Yes, I forgot. These are in fact the only two magical methods provided by ModulesWithBasis.ElementMethods.

---

New commits:

​2fd7595

Add FiniteDimensionalAlgebraElement.__len__

---

New commits:

​2fd7595

Add FiniteDimensionalAlgebraElement.__len__

[tscrim]

We definitely should override __getitem__ as that will be a lot faster than the generic implementation. However, there is a change in behavior:

sage: A.an_element()[5]
IndexError: matrix index out of range

Before this would return 0. It is arguable that it is the wrong behavior, and you don't have to change it, but I am just pointing it out. However, the different behavior does make the iterator work, but that is a matter of a lack of implementation. Really, there should be an implementation at the category level, but this would create an inconsistency with free module elements (i.e., created by vector()) and subclasses of IndexedFreeModuleElement.

I think it would be good to decide on and implement the behavior for __iter__ here:

1 - Iterate over the vector and yield the coefficients. 2 - Iterate over pairs (k, v), where k is the index/key and v is a non-zero coefficient.

My vote would be for 2 to be consistent with many (all?) the other algebras in Sage.

Subsequently, your __len__ is a change in behavior, but it needs to be consistent with the __iter__. So we will need to choose.

[tscrim]

Oh, one other thing that extension classes without a __dict__ cannot do is add arbitrary attributes. So this no longer works:

sage: A.an_element().rename("foo")
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-6-996df49490f5> in <module>()
----> 1 _.rename("foo")

/home/travis/sage/src/sage/structure/sage_object.pyx in sage.structure.sage_object.SageObject.rename (/home/travis/sage/src/build/cythonized/sage/structure/sage_object.c:2354)()
    171                 self.__custom_name = str(x)
    172             except AttributeError:
--> 173                 raise NotImplementedError("object does not support renaming: %s" % self)
    174 
    175     def reset_name(self):

NotImplementedError: object does not support renaming: e0

However, I am okay with that not working as that is merely a convenience feature and most likely never used on these elements.

[SimonKing]

Replying to tscrim:

Oh, one other thing that extension classes without a __dict__ cannot do is add arbitrary attributes. So this no longer works:

sage: A.an_element().rename("foo")

Seriously? Strange. I thought that renaming is supported on the level of SageObject. Aha, I just read:

           To support them for a specific class, add a
           ``cdef public __custom_name`` attribute.

I find it strange that it has not been added on the level of SageObject.

[SimonKing]

Replying to tscrim:

We definitely should override __getitem__ as that will be a lot faster than the generic implementation. However, there is a change in behavior:

sage: A.an_element()[5]
IndexError: matrix index out of range

Before this would return 0. It is arguable that it is the wrong behavior,

I do believe not raising an IndexError is wrong in that case.

However, the different behavior does make the iterator work, but that is a matter of a lack of implementation.

Here I am unsure what you are referring to by "different behaviour". As a matter of fact, iterator works now:

sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x = A.1
sage: list(x)
[0, 1, 0]

I think it would be good to decide on and implement the behavior for __iter__ here:

1 - Iterate over the vector and yield the coefficients. 2 - Iterate over pairs (k, v), where k is the index/key and v is a non-zero coefficient.

My vote would be for 2 to be consistent with many (all?) the other algebras in Sage.

I don't think it is possible to be consistent, as even multi- and univariate polynomials aren't:

sage: R.<x,y> = QQ[]
sage: list(x^2+2*y)
[(1, x^2), (2, y)]
sage: R.<x> = QQ[]
sage: list(x^8)
[0, 0, 0, 0, 0, 0, 0, 0, 1]

Subsequently, your __len__ is a change in behavior, but it needs to be consistent with the __iter__. So we will need to choose.

A finite dimensional algebra is a finite dimensional vector space, and thus I believe one should be consistent with vectors, not with multivariate polynomials. And that is the case with the current implementation.

[SimonKing]

Replying to SimonKing:

Replying to tscrim:

My vote would be for 2 to be consistent with many (all?) the other algebras in Sage.

I don't think it is possible to be consistent, as even multi- and univariate polynomials aren't:

sage: R.<x,y> = QQ[]
sage: list(x^2+2*y)
[(1, x^2), (2, y)]
sage: R.<x> = QQ[]
sage: list(x^8)
[0, 0, 0, 0, 0, 0, 0, 0, 1]

Also note that it is (v, k) not (k, v).

[tscrim]

I wasn't thinking of polynomial rings (which are not [currently] in the category of ModulesWithBasis), but of subclasses of CombinatorialFreeModule:

sage: E = algebras.Exterior(QQ, 2)
sage: E.an_element()
2*e0 + 4*e1 + 1
sage: list(E.an_element())
[((), 1), ((1,), 4), ((0,), 2)]

sage: F = algebras.Free(QQ, 3, 'x')
sage: list(F.an_element())
[(x1, 3), (x0, 2), (1, 2)]
sage: F.an_element()
2 + 2*x0 + 3*x1

[tscrim]

Replying to SimonKing:

Replying to tscrim:

We definitely should override __getitem__ as that will be a lot faster than the generic implementation. However, there is a change in behavior:

sage: A.an_element()[5]
IndexError: matrix index out of range

Before this would return 0. It is arguable that it is the wrong behavior,

I do believe not raising an IndexError is wrong in that case.

It probably is a side-effect of a generic implementation of _coefficient_fast, which is called from __getitem__, that doesn't do any safety checks (because it needs to support infinite bases). So, yes, I agree that raising an error is the right thing.

However, the different behavior does make the iterator work, but that is a matter of a lack of implementation.

Here I am unsure what you are referring to by "different behaviour". As a matter of fact, iterator works now:

sage: A = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0,0], [0,1,0], [0,0,0]]), Matrix([[0,1,0], [0,0,0], [0,0,0]]), Matrix([[0,0,0], [0,0,0], [0,0,1]])])
sage: x = A.1
sage: list(x)
[0, 1, 0]

The old behavior for __getitem__ returned 0, where there seems to be a default __iter__ by Python that uses __getitem__ with ints until it raises an IndexError.

I think it would be good to decide on and implement the behavior for __iter__ here:

1 - Iterate over the vector and yield the coefficients. 2 - Iterate over pairs (k, v), where k is the index/key and v is a non-zero coefficient.

My vote would be for 2 to be consistent with many (all?) the other algebras in Sage.

I don't think it is possible to be consistent, as even multi- and univariate polynomials aren't:

sage: R.<x,y> = QQ[]
sage: list(x^2+2*y)
[(1, x^2), (2, y)]
sage: R.<x> = QQ[]
sage: list(x^8)
[0, 0, 0, 0, 0, 0, 0, 0, 1]

Good point. Always fun the differences between uni- and multivariate polynomials.

Subsequently, your __len__ is a change in behavior, but it needs to be consistent with the __iter__. So we will need to choose.

A finite dimensional algebra is a finite dimensional vector space, and thus I believe one should be consistent with vectors, not with multivariate polynomials. And that is the case with the current implementation.

See comment:72 for other algebras. Polynomials error out when trying to do, e.g., len(x^2+x) and sidestep the issue. It comes down to len meaning length as a column vector or length of the support.

---

I am okay with the current behavior, even if I would prefer it to agree with ModulesWithBasis objects/CombinatorialFreeModule subclasses, but I just wanted to bring up that it is a change. Since you prefer the current implementation and there is a green patchbot, I am setting a positive review.

[SimonKing]

Thank you!

[tscrim]

Thank you for your work on this. I look forward to having the quotients as well.

[jdemeyer]

You are using cdef int in a few places where you might want cdef Py_ssize_t. Otherwise you are restricting to a size of 2³¹ for no good reason.

[SimonKing]

Replying to jdemeyer:

You are using cdef int in a few places where you might want cdef Py_ssize_t. Otherwise you are restricting to a size of 2³¹ for no good reason.

The restriction has already been there before:

sage: M = MatrixSpace(GF(2),1,2^32)()
sage: M[0,2^32-1]
---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
<ipython-input-2-0b1808df77ca> in <module>()
----> 1 M[Integer(0),Integer(2)**Integer(32)-Integer(1)]

/home/king/Sage/git/sage/src/sage/matrix/matrix0.pyx in sage.matrix.matrix0.Matrix.__getitem__ (build/cythonized/sage/matrix/matrix0.c:7018)()
    946                 if not PyIndex_Check(col_index):
    947                     raise TypeError("index must be an integer")
--> 948                 col = col_index
    949                 if col < 0:
    950                     col += ncols

OverflowError: value too large to convert to int

[SimonKing]

Replying to jdemeyer:

You are using cdef int in a few places where you might want cdef Py_ssize_t. Otherwise you are restricting to a size of 2³¹ for no good reason.

And doesn't Py_ssize_t denote some sort of signed integer? I guess I need an unsigned one. Py_size_t?

[SimonKing]

One could even argue that using cdef int prevents one from crashing:

sage: M = MatrixSpace(GF(2),2^32)
sage: M(1)
------------------------------------------------------------------------
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/cysignals/signals.so(+0x5e75)[0x7fb34704ee75]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/cysignals/signals.so(+0x5ec5)[0x7fb34704eec5]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/cysignals/signals.so(+0x8e97)[0x7fb347051e97]
/lib/x86_64-linux-gnu/libpthread.so.0(+0x11390)[0x7fb35378d390]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/matrix/matrix_mod2_dense.so(+0xb102)[0x7fb1b8d71102]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/matrix/matrix0.so(+0x407ee)[0x7fb1bc5a67ee]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x48bd)[0x7fb353aa18ad]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/misc/cachefunc.so(+0x16433)[0x7fb34572a433]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/misc/cachefunc.so(+0x166c7)[0x7fb34572a6c7]
/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/misc/cachefunc.so(+0x1698c)[0x7fb34572a98c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7fb3539ed483]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x56da)[0x7fb353aa26ca]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(+0x8567c)[0x7fb353a1e67c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7fb3539ed483]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(+0x657ec)[0x7fb3539fe7ec]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7fb3539ed483]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(+0xc1985)[0x7fb353a5a985]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyObject_Call+0x43)[0x7fb3539ed483]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x56da)[0x7fb353aa26ca]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCode+0x19)[0x7fb353aa69a9]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x7428)[0x7fb353aa4418]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalFrameEx+0x8020)[0x7fb353aa5010]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCodeEx+0x81c)[0x7fb353aa688c]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyEval_EvalCode+0x19)[0x7fb353aa69a9]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyRun_FileExFlags+0x8a)[0x7fb353acaa4a]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(PyRun_SimpleFileExFlags+0xd7)[0x7fb353acbdf7]
/home/king/Sage/git/sage/local/lib/libpython2.7.so.1.0(Py_Main+0xc3e)[0x7fb353ae24be]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf0)[0x7fb352cc2830]
python(_start+0x29)[0x400719]
------------------------------------------------------------------------
Attaching gdb to process id 6413.

Failed to run gdb.
Failed to run gdb.
Install gdb for enhanced tracebacks.
------------------------------------------------------------------------
Unhandled SIGSEGV: A segmentation fault occurred.
This probably occurred because a *compiled* module has a bug
in it and is not properly wrapped with sig_on(), sig_off().
Python will now terminate.
------------------------------------------------------------------------
Segmentation fault (core dumped)

[jdemeyer]

Replying to SimonKing:

The restriction has already been there before:

The fact that a certain bug exists already is not a good reason to add more bugs of a similar nature.

[jdemeyer]

Replying to SimonKing:

And doesn't Py_ssize_t denote some sort of signed integer?

Yes, it does. However, Python generally uses Py_ssize_t for sizes (say, of lists or tuples). Also, in Sage, the number of rows/columns of a matrix is already stored as Py_ssize_t:

cdef class Matrix(ModuleElement):
    # All matrix classes must be written in Cython
    cdef Py_ssize_t _nrows
    cdef Py_ssize_t _ncols

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3cc3a27

Use Py_ssize_t to iterate over matrix indices

[SimonKing]

Replying to jdemeyer: Also, in Sage, the number of rows/columns of a matrix is already stored as Py_ssize_t:

cdef class Matrix(ModuleElement):
    # All matrix classes must be written in Cython
    cdef Py_ssize_t _nrows
    cdef Py_ssize_t _ncols

OK, that's a good reason. I changed to Py_ssize_t.

[tscrim]

Jeroen, any more comments (or things not addressed)?

[jdemeyer]

Not really. Anyway, I have not looked deeply into this ticket, I was just pointing a few random things.

[tscrim]

I let this drop off. Whoops. Back to positive review.