sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.33k stars 453 forks source link

revert for LazyTaylorSeries and LazySymmetricFunction is missing #34383

Closed mantepse closed 1 year ago

mantepse commented 2 years ago
sage: L.<z> = LazyTaylorSeriesRing(ZZ)
sage: (z-z^2).revert()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
...
AttributeError: 'LazyTaylorSeriesRing_with_category.element_class' object has no attribute 'revert'

Same for LazySymmetricFunctions. compositional_inverse might be a good alias.

Depends on #32324 Depends on #34453 Depends on #34494

CC: @tscrim

Component: combinatorics

Keywords: LazyPowerSeries

Author: Martin Rubey

Branch/Commit: f3f011f

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/34383

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 0d9f02f to ae99175

DaveWitteMorris commented 2 years ago
comment:38

I created #34441 for the issue about ZZ.variable_names() that was mentioned in comment:28 and comment:30. IMHO it's a bug.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from ae99175 to baff215

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

baff215plethysm with tensor products, part 2 of 2
mantepse commented 2 years ago
comment:40

The only missing method is reversion for Taylor series now.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from baff215 to 9f200b4

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

dff4fc4Adding a completion of graded algebras by formal series.
4dd1217Merge branch 'public/rings/lazy_graded_algebras-34407' of trac.sagemath.org:sage into t/34383/revert_for_lazytaylorseries_and_lazysymmetricfunction_is_missing
9f200b4fix oversights
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

43adffffix bug in LazyTaylorSeries.polynomial, add revert
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 9f200b4 to 43adfff

mantepse commented 2 years ago
comment:43

There remains a bug in LazyTaylorSeries.__call__:

At the very beginning, we deal with the case where f is a polynomial:

        # f has finite length
        if isinstance(self._coeff_stream, Stream_exact) and not self._coeff_stream._constant:
            # constant polynomial
            poly = self.polynomial()
            if poly.is_constant():
                return P(poly)
            return P(poly(g))

However, this is too eager. In the example below, we would like to compute g.define(z / ((1-z)(g))), which fails, whereas g.define(z / (1-g)) works. Also, if we compute in the LazyLaurentSeriesRing instead, there is no problem.

It would be great if LazyLaurentSeries.__call__ and LazyTaylorSeries.__call__ and perhaps even LazySymmetricFunction.__call__ could share some code, to avoid these inconsistencies.

Example:

sage: L.<z> = LazyTaylorSeriesRing(QQ)
sage: f = 1 - z
sage: g = L(None, valuation=1)
sage: g.define(z / (1 - g))
sage: g
z + z^2 + 2*z^3 + 5*z^4 + 14*z^5 + 42*z^6 + 132*z^7 + O(z^8)

sage: g = L(None, valuation=1)
sage: g.define(z / f(g))
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/sage-develop/src/sage/data_structures/stream.py:329, in Stream_inexact.__getitem__(self, n)
    328 try:
--> 329     c = self._cache[n]
    330 except KeyError:

KeyError: 1

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 g.define(z / f(g))

File ~/sage-develop/src/sage/rings/lazy_series.py:3508, in LazyTaylorSeries.__call__(self, check, *g)
   3506     if poly.is_constant():
   3507         return P(poly)
-> 3508     return P(poly(g))
   3510 # f now has (potentially) infinitely many terms
   3511 # Lift the resulting parent to a lazy series (if possible)
   3512 # Also make sure each element of g is a LazyModuleElement
   3513 from sage.rings.polynomial.polynomial_ring import PolynomialRing_general

File ~/sage-develop/src/sage/rings/polynomial/polynomial_rational_flint.pyx:552, in sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint.__call__()
    550             return acb_z
    551 
--> 552     return Polynomial.__call__(self, *x, **kwds)
    553 
    554 cpdef Polynomial truncate(self, long n):

File ~/sage-develop/src/sage/rings/polynomial/polynomial_element.pyx:881, in sage.rings.polynomial.polynomial_element.Polynomial.__call__()
    879 d = pol.degree()
    880 
--> 881 if d <= 0 or (isinstance(a, Element) and R.is_exact() and a.is_zero()):
    882     return cst # with the right parent thanks to the above coercion
    883 elif pol._parent is R and a.is_gen():

File ~/sage-develop/src/sage/structure/element.pyx:1064, in sage.structure.element.Element.is_zero()
   1062         implement ``__bool__`` instead.
   1063     """
-> 1064     return not self
   1065 
   1066 def _cache_key(self):

File ~/sage-develop/src/sage/rings/lazy_series.py:618, in LazyModuleElement.__bool__(self)
    616     if any(a for a in self._coeff_stream._cache):
    617         return True
--> 618 if self[self._coeff_stream._approximate_order]:
    619     return True
    620 raise ValueError("undecidable as lazy Laurent series")

File ~/sage-develop/src/sage/rings/lazy_series.py:254, in LazyModuleElement.__getitem__(self, n)
    252     step = n.step if n.step is not None else 1
    253     return [R(self._coeff_stream[k]) for k in range(start, n.stop, step)]
--> 254 return R(self._coeff_stream[n])

File ~/sage-develop/src/sage/data_structures/stream.py:331, in Stream_inexact.__getitem__(self, n)
    329         c = self._cache[n]
    330     except KeyError:
--> 331         c = self.get_coefficient(n)
    332         self._cache[n] = c
    333 else:

File ~/sage-develop/src/sage/data_structures/stream.py:876, in Stream_uninitialized.get_coefficient(self, n)
    857 def get_coefficient(self, n):
    858     """
    859     Return the ``n``-th coefficient of ``self``.
    860 
   (...)
    874         1
    875     """
--> 876     return self._target[n]

TypeError: 'NoneType' object is not subscriptable
tscrim commented 2 years ago
comment:44

Actually, it is not a bug in the lazy Taylor series code. In fact, it i working “correctly” because it is trying a very basic test to see if it is zero or not, the test of which is out of our hands, which at the time it does not know. I would actually expect the same problem in the LazyLaurentSeriesRing.

To fix this, the quick-n-dirty fix would be to have an uninitialized series have __bool__() simply return True as it represents a variable that has yet to be specialized . This would make it consistent with variables of a polynomial ring or symbolic variables. )Plus, in nearly all use cases, this will end up defining a non-zero series anyways.) Actually, because of that, I would even say this is the correct thing to do, and not simply a hack solution.

A more broader change to fix the problem would be to have some kind of lazily-evaluated version of composition. Not only would it increase the code complexity when there is arguably a valid solution, it also would not be compatible with polynomials that can be evaluated at elements of the lazy series ring R with the result being coerced into R. Note that the polynomial itself does ‘’not’’ even need to convert (much less coerce) into R. So this would just sweep this problem under the rug (this includes even the extreme change of removing the exact implementation or feeding off to the polynomial code).

Having written all of this down, I am fully convinced we should just make it work by saying an uninitialized series is considered to be non-zero by having __bool__() return True in that case.

mantepse commented 2 years ago
comment:45

I'm afraid I simply do not understand your argument. You say that you "would actually expect the same problem in the LazyLaurentSeriesRing, but I don't see the problem there:

sage: L.<z> = LazyLaurentSeriesRing(QQ)
sage: f = 1-z
sage: g = L(None, valuation=1)
sage: g.define(z / f(g))
sage: g
z + z^2 + 2*z^3 + 5*z^4 + 14*z^5 + 42*z^6 + 132*z^7 + O(z^8)

In fact, this is (essentially) a doctest.

In LazyLaurentSeries.__call__, if f is a polynomial, we use the following test before composing the polynomial with g:

I don't see where g.__bool__ would come into play here.

mantepse commented 2 years ago
comment:46

I tried to do a consistency check. I'm sure it could be improved.

    sage: def check(L, z, verbose=False):
    ....:     # division
    ....:     lf = [0, L(0), 1, L(1), z, 1 + z, 2 + z + z^2]
    ....:     lg = [3, L(3), 1 + z, 2 + z + z^2]
    ....:     for f in lf:
    ....:         for g in lg:
    ....:             try:
    ....:                 h = f / g
    ....:                 if verbose: print("(%s) / (%s) = %s" % (f, g, h))
    ....:             except Exception as e:
    ....:                 print("%s in (%s) / (%s)" % (e, f, g))
    ....:     # composition
    ....:     f = 2 + z + z^2
    ....:     l = [(f, 0), (f, L(0)), (f, 2 + z + z^2), (f, L(3)/(1 - 2*z))]
    ....:     f = L(3)/(2 - 3*z)
    ....:     l.extend([(f, 0), (f, L(0)), (f, 3*z/(1 - 2*z))])
    ....:     for f, g in l:
    ....:         try:
    ....:             h = f(g)
    ....:             if verbose: print("(%s)(%s) = %s" % (f, g, h))
    ....:         except Exception as e:
    ....:             print("%s in (%s)(%s)" % (e, f, g))
    ....:     # reversion
    ....:     l = [2 + 3*z, 3*z + 2*z^2, 3*z/(1 - 2*z - 3*z^2)]
    ....:     for f in l:
    ....:         try:
    ....:             h = f.revert()
    ....:             if verbose: print("(%s)^{(-1)} = %s" % (f, h))
    ....:         except Exception as e:
    ....:             print("%s in (%s).revert()" % (e, f))

The current results are as follows:

sage: L.<z> = LazyLaurentSeriesRing(QQ)
sage: check(L, z)

sage: L.<z> = LazyTaylorSeriesRing(QQ)
sage: check(L, z)
unable to evaluate the series at (0,) in (3/2 + 9/4*z + 27/8*z^2 + 81/16*z^3 + 243/32*z^4 + 729/64*z^5 + 2187/128*z^6 + O(z^7))(0)
'NoneType' object is not subscriptable in (3*z + 2*z^2).revert()

sage: p = SymmetricFunctions(QQ).p()
sage: L = LazySymmetricFunctions(p)
sage: check(L, L(p[1]))
can only compose with a positive valuation series in (2*p[] + p[1] + (p[1,1]))(3*p[] + 6*p[1] + (12*p[1,1]) + (24*p[1,1,1]) + (48*p[1,1,1,1]) + (96*p[1,1,1,1,1]) + (192*p[1,1,1,1,1,1]) + O^7)
no conversion of this rational to integer in (3/2*p[] + 9/4*p[1] + (27/8*p[1,1]) + (81/16*p[1,1,1]) + (243/32*p[1,1,1,1]) + (729/64*p[1,1,1,1,1]) + (2187/128*p[1,1,1,1,1,1]) + O^7)(0)
tscrim commented 2 years ago
comment:47

Replying to @mantepse:

I'm afraid I simply do not understand your argument. You say that you "would actually expect the same problem in the LazyLaurentSeriesRing, but I don't see the problem there:

sage: L.<z> = LazyLaurentSeriesRing(QQ)
sage: f = 1-z
sage: g = L(None, valuation=1)
sage: g.define(z / f(g))
sage: g
z + z^2 + 2*z^3 + 5*z^4 + 14*z^5 + 42*z^6 + 132*z^7 + O(z^8)

In fact, this is (essentially) a doctest.

In LazyLaurentSeries.__call__, if f is a polynomial, we use the following test before composing the polynomial with g:

  • if g is not a LazyModuleElement, we compose
  • otherwise, if g is a LazyCauchyProductSeries and g._coeff_stream is Stream_exact and not g._coeff_stream._constant, we also do a special case
  • finally, if none of the above applies, we use the generic algorithm.

In some respects, we get lucky here. I would expect the following to work:

sage: L.<z> = LazyLaurentSeriesRing(QQ)
sage: f = 1 - z
sage: g = L(None, valuation=1)
sage: f(g) == f.polynomial()(g)

In particular, the f.polynomial()(g) fails in exactly the same way as above.

The only reason we need to do anything special is for the negative exponents when g is exact. Everything else is just an optimization (as the comment there indicates, it should be replaceable by return poly(g)) to take advantage of caching (it could also benefit by reducing the number of multiplications too).

I don't see where g.__bool__ would come into play here.

It is there in the error message as there is a call to is_zero().

tscrim commented 2 years ago
comment:48

I am a bit happy as it seems on a simple example that my optimization is faster:

sage: L.<z> = LazyLaurentSeriesRing(QQ)
sage: f = z^-1 + 2 + z + z^2 + 2*z^3
sage: g = z^-2 / (1 - z)

sage: %time (f(g))[300]
CPU times: user 25.6 ms, sys: 10 µs, total: 25.6 ms
Wall time: 25.2 ms
94862
sage: %time (f.polynomial()(g))[300]
CPU times: user 92.1 ms, sys: 0 ns, total: 92.1 ms
Wall time: 91.8 ms
94862

sage: %time (f(g))[3000]
CPU times: user 2.22 s, sys: 0 ns, total: 2.22 s
Wall time: 2.22 s
9048062
sage: %time (f.polynomial()(g))[3000]
CPU times: user 8.2 s, sys: 0 ns, total: 8.2 s
Wall time: 8.2 s
9048062

All come behold the power of the cache and rejoice. :P

tscrim commented 2 years ago
comment:49

Here is the fix with doctests added.


New commits:

cfec1e0Making uninitialized series g have bool(g) == True. Better define() semantic with exact and non lazy series elements.
tscrim commented 2 years ago

Changed branch from u/mantepse/revert_for_lazytaylorseries_and_lazysymmetricfunction_is_missing to public/rings/lazy_series_revert-34383

tscrim commented 2 years ago

Changed commit from 43adfff to cfec1e0

tscrim commented 2 years ago

Reviewer: Travis Scrimshaw

mantepse commented 2 years ago
comment:50

Cool, and thank you for the explanation! I missed the fact that composing a polynomial with an uninitialised series actually makes sense even when using the generic polynomial code.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

bdc4537make the linter happier
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from cfec1e0 to bdc4537

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

07ca029add a consistency check
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from bdc4537 to 07ca029

mantepse commented 2 years ago
comment:53

There are still some tests which (trivially) fail. I am going to work on them for the next 2 hours or so.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 07ca029 to 1d559cc

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

1d559ccimprove composing Taylor series, start making plethysm more general
mantepse commented 2 years ago
comment:55

Taylor series should be OK now.

There are still slight problems with plethysm. Some of them might be best fixed in the implementation of ordinary symmetric functions:

sage: p = SymmetricFunctions(QQ).p()
sage: p[2,2,1](3)
...
TypeError: only know how to compute plethysms between symmetric functions or tensors of symmetric functions

However, I am not completely sure who should handle plethysms of the form f o g where f is an ordinary symmetric function and g is a lazy symmetric function. The situation is similar to the case of polynomials, but not the same. Essentially, we need to be able to express g as a (possibly constant, possibly lazy) symmetric function in the power sum basis.

One possibility would be to handle this in Stream_plethysm. However, this goes against the general philosophy that the "exact" case is handled in lazy_series.py. On the other hand, implementing yet another version of plethysm in LazySymmetricFunction.__call__ seems like a lot of code duplication to me.

tscrim commented 2 years ago
comment:56

Replying to @mantepse:

Taylor series should be OK now.

Great!

There are still slight problems with plethysm. Some of them might be best fixed in the implementation of ordinary symmetric functions:

sage: p = SymmetricFunctions(QQ).p()
sage: p[2,2,1](3)
...
TypeError: only know how to compute plethysms between symmetric functions or tensors of symmetric functions

Yea, this one annoys me from time to time that it doesn't at least try to coerce the element into a symmetric functions. This is for a separate ticket, but +1 for trying the coercion (or at least handling elements from the base ring)

However, I am not completely sure who should handle plethysms of the form f o g where f is an ordinary symmetric function and g is a lazy symmetric function. The situation is similar to the case of polynomials, but not the same. Essentially, we need to be able to express g as a (possibly constant, possibly lazy) symmetric function in the power sum basis.

The symmetric function f definitely should. There is no mechanism within Python (or within the coercion framework) to automatically pass control over to g.

One possibility would be to handle this in Stream_plethysm. However, this goes against the general philosophy that the "exact" case is handled in lazy_series.py. On the other hand, implementing yet another version of plethysm in LazySymmetricFunction.__call__ seems like a lot of code duplication to me.

I don't know how Stream_plethysm could even get control of this without having f actually handle the object. At which point, there is an easy solution: convert f to a LazySymmetricFunction and run __call__ on that. I will push a fix in just a minute.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

c47b060Doing some changes to normal symmetric functions plethysm to support coefficients and lazy symmetric functions.
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 1d559cc to c47b060

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

5776c61Adding the category of commutative rings to the category of tensor products of commutative algebras.
5f21f3dMerge branch 'public/categories/tensor_commutative_rings-34453' into public/rings/lazy_series_revert-34383
a2e9ed7Making the test for coercion of plethysms wrt tesnor products more robust.
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from c47b060 to a2e9ed7

tscrim commented 2 years ago
comment:59

Replying to @tscrim:

Replying to @mantepse:

There are still slight problems with plethysm. Some of them might be best fixed in the implementation of ordinary symmetric functions:

sage: p = SymmetricFunctions(QQ).p()
sage: p[2,2,1](3)
...
TypeError: only know how to compute plethysms between symmetric functions or tensors of symmetric functions

Yea, this one annoys me from time to time that it doesn't at least try to coerce the element into a symmetric functions. This is for a separate ticket, but +1 for trying the coercion (or at least handling elements from the base ring)

I ended up fixing this here as it was basically part of the changes I needed to make for the lazy series. We will not be able to return an element in the original parent for "constants" (broadly interpreted) without doing some cumbersome (and well-definedness is another issue too).

The check for tensor products is a bit brittle, so I've tried to make it a bit more robust (say, if the commutative ring is also a tensor product).

In trying to find a doctest for this, I found out that tensor products of commutative algebras do not know they are commutative rings:

sage: X = algebras.Shuffle(QQ, 'ab')
sage: Y = algebras.Shuffle(QQ, 'bc')
sage: X in CommutativeRings()
True
sage: Y in CommutativeRings()
True
sage: T = tensor([X,Y])
sage: T in CommutativeRings()
False

This is now #34453. While not a strict dependency for this ticket, it means I get to write the doctest I wanted.

One possibility would be to handle this in Stream_plethysm. However, this goes against the general philosophy that the "exact" case is handled in lazy_series.py. On the other hand, implementing yet another version of plethysm in LazySymmetricFunction.__call__ seems like a lot of code duplication to me.

I don't know how Stream_plethysm could even get control of this without having f actually handle the object. At which point, there is an easy solution: convert f to a LazySymmetricFunction and run __call__ on that. I will push a fix in just a minute.

I have now done this.

tscrim commented 2 years ago

Changed dependencies from #32324 to #32324, #34453

mantepse commented 2 years ago
comment:60

So, the idea is now to make Stream_plethysm handle the plethysm of a symmetric function with finite support and lazy symmetric functions of valuation 0?

(That's fine with me, I just want confirmation.)

(sorry for being short, I just woke up :-)

tscrim commented 2 years ago
comment:61

Replying to @mantepse:

So, the idea is now to make Stream_plethysm handle the plethysm of a symmetric function with finite support and lazy symmetric functions of valuation 0?

Essentially yes because it already must do so.

(sorry for being short, I just woke up :-)

No problem.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from a2e9ed7 to 94219d4

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

94219d4fix plethysm with f having finite support
mantepse commented 2 years ago
comment:63

It seems that we are done for now!

I am sure that code and especially tests and documentation could and should be polished, but maybe it would be good to move forward now. At least concerning functionality, the code is a huge improvement over what we had before this ticket!

Perhaps it would make sense to make more tests systematic, and also check which tests are now redundant.

I am very grateful for all your help and expertise, Travis!

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 94219d4 to 6745cf1

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

6745cf1Some last details and optimizations to Stream_plethysm.
tscrim commented 2 years ago
comment:65

Thank you for the changes and all your hard work.

Some comments about my changes:

Isn't the := syntax great? I had been missing that so much.

If my changes are good, then I we can set a positive review and move on to the next thing. We can do more polishing of the documentation and tests later.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

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

3d6d39ffix typo in comment
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 6745cf1 to 3d6d39f

mantepse commented 2 years ago
comment:67

Perfect, thanks!

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 3d6d39f to dafdcb4

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

dafdcb4fix doctests
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

731072aremove unnecessary code leftover from merge, and remove unneccesary import detected by pyflakes
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from dafdcb4 to 731072a