Implement integration and Taylor series for lazy series #36233
Conversation
|
This needs a rebase, but it's not urgent for me. |
tscrim
force-pushed
the
lazy_series/integration
branch
from
adc4538 to
98a354e
Compare
|
Here is the rebased version. No problem. This was something that was useful for solving differential equations as FPS. Actually, that makes me wonder if we should add a method to the lazy (Laurent) series parent to return the Taylor expansion of a function as a lazy series. |
Yes, that would certainly be useful. Another method I had in FriCAS which is extremely useful is a very general (and therefore possibly slow) facility to compute the solution of a functional equation, given initial values. For example: It might be even nicer if one could replace My implementation simply created a power series with undetermined coefficients and complained whenever it could not compute the next coefficient by either solving a linear equation or using an initial value. Possibly one could even relax "linear equation" to "equation with unique solution". |
|
keyword: LazyPowerSeries |
tscrim
force-pushed
the
lazy_series/integration
branch
from
98a354e to
419701c
Compare
|
Here is an initial version of the functional equation. It has not been well tested, but I got a version that uses the undefined series as the variable. It might not be the most robust, but I don't have good examples to really test it on. I also added the Taylor series method. |
|
This is extremely cool! I have lots of examples, but I should be doing other work this week. If you want to, here are examples to try:
|
|
Thank you. Those were good tests, which uncovered some bugs. I think the overall effect is that it is now slower since I have to substitute what I know on each coefficient (rather than directly plugging it in beforehand), but it is now more robust. I also added a test for a Laurent series function equation solution. The next thing to consider is what the best API is. Right now this is a bit clunky, but I don't know how to make it better. This also won't be able to do things with general input (such as using generic functions), but it should cover most reasonable cases... |
|
I just started to play with it. It is wonderful, I am impressed! Here are some very vague ideas concerning the user interface:
So maybe, as a method of Other things: should possibly be Here is a bug that I do not understand: There should be a solution, however (this is Jacobi, see Stanley EC2, page 282, exercise 6.63b): |
|
PS: I tested all the other examples I had, without any further bugs. I am impressed. |
|
Another bug, from a current collaboration with M. Ludwig and A. Freyer: However, here is a solution: Moreover, the following might also be a bug: |
|
I looked a little at the code.
Now |
|
It seems to me that the following should be possible: Think of the desired series Hm, I just see that modifying the caches of Then, computing the next term of |
|
Here is a proof of concept. After applying the patch below, do: diff --git a/src/sage/data_structures/stream.py b/src/sage/data_structures/stream.py
index e5de010fbad..0f023aee101 100644
--- a/src/sage/data_structures/stream.py
+++ b/src/sage/data_structures/stream.py
@@ -147,6 +147,12 @@ class Stream():
"""
self._true_order = true_order
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ pass
+
@lazy_attribute
def _approximate_order(self):
"""
@@ -447,6 +453,15 @@ class Stream_inexact(Stream):
yield self.get_coefficient(n)
n += 1
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ if self._is_sparse:
+ self._cache = {i: function(v) for i, v in self._cache.items()}
+ else:
+ self._cache = [function(v) for v in self._cache]
+
def order(self):
r"""
Return the order of ``self``, which is the minimum index ``n`` such
@@ -1277,6 +1292,14 @@ class Stream_uninitialized(Stream_inexact):
self._initializing = False
return result
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ if self._target is not None:
+ self._target.map_cache(function)
+
class Stream_functional_equation(Stream_inexact):
r"""
@@ -1439,6 +1462,13 @@ class Stream_unary(Stream_inexact):
self._series = series
super().__init__(is_sparse, true_order)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._series.map_cache(function)
+
def __hash__(self):
"""
Return the hash of ``self``.
@@ -1554,6 +1584,14 @@ class Stream_binary(Stream_inexact):
self._right = right
super().__init__(is_sparse, False)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._left.map_cache(function)
+ self._right.map_cache(function)
+
def __hash__(self):
"""
Return the hash of ``self``.
@@ -2370,6 +2408,13 @@ class Stream_plethysm(Stream_binary):
f = Stream_map_coefficients(f, lambda x: p(x), is_sparse)
super().__init__(f, g, is_sparse)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._powers = [g.map_cache(function) for g in self._powers]
+
@lazy_attribute
def _approximate_order(self):
"""
@@ -3251,6 +3296,12 @@ class Stream_shift(Stream):
self._shift = shift
super().__init__(series._true_order)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ self._series.map_cache(function)
+
@lazy_attribute
def _approximate_order(self):
""" |
|
TL;DR, I am not sure which would be best. I am not opposed to switching to the cache sub if you definitively think it is better. In an ideal world, we would be able to solve things out of order, but knowing how many coefficient to compute in order to get a unique solution is an undecidable problem. However you are right that the code does not currently check that the next variable is the next in order., which is the only sensible thing to do for a computation to succeed. However, from your examples, I was able to quash a number of bugs. I had to do some hacks for the brittleness of the (infinite) polynomial ring substitution, which took more time than almost anything... I also removed the |
|
I don't know whether substituting elements in the cache is faster - I would think that it is if the dependence on previous coefficients is complicated, and particularly so if after substitution we have just rational numbers. But I guess it would be necessary to code it and try some examples. I would expect that the Jacobi example would be interesting. Currently it takes already seconds to compute the coefficient of 49, the number of monomials before substitution is about 150000, whereas after substitution it is (of course) just 1. In any case, I am quite convinced meanwhile that it would be good to call the method |
|
I shouldn't have done this, but at least I have some success. Computing This is really only a proof of concept. In particular, this only works if it is indeed enough to check only the last element of the cache (which I believe is true). Initially I thought that I would have to go through all of the cache, which would make things extremely slow. In any case, I dislike my implementation, there ought to be a better way. Moreover, there is at least one bug which I do not understand yet - when you execute the cosine example on line 2599 of Finally, there seems to be a performance bug in the diff --git a/src/sage/data_structures/stream.py b/src/sage/data_structures/stream.py
index 88c765569b..87e5042eeb 100644
--- a/src/sage/data_structures/stream.py
+++ b/src/sage/data_structures/stream.py
@@ -147,6 +147,12 @@ class Stream():
"""
self._true_order = true_order
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ pass
+
@lazy_attribute
def _approximate_order(self):
"""
@@ -447,6 +453,17 @@ class Stream_inexact(Stream):
yield self.get_coefficient(n)
n += 1
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ if self._cache:
+ if self._is_sparse:
+ i = max(self._cache)
+ self._cache[i] = function(self._cache[i])
+ else:
+ self._cache[-1] = function(self._cache[-1])
+
def order(self):
r"""
Return the order of ``self``, which is the minimum index ``n`` such
@@ -1299,6 +1316,14 @@ class Stream_uninitialized(Stream_inexact):
ret._target = temp
return ret
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ if self._target is not None:
+ self._target.map_cache(function)
+
class Stream_functional_equation(Stream_inexact):
r"""
@@ -1343,12 +1368,23 @@ class Stream_functional_equation(Stream_inexact):
initial_values = []
super().__init__(False, true_order)
self._F = F
- self._base = R
self._initial_values = initial_values
self._approximate_order = approximate_order
+
+ from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing
+ self._P = InfinitePolynomialRing(R, names=('FESDUMMY',))
+ x = self._P.gen()
+ def get_coeff(n):
+ if n < len(self._initial_values):
+ return self._initial_values[n]
+ return x[n]
+
self._uninitialized = uninitialized
self._uninitialized._approximate_order = approximate_order
- self._uninitialized._target = self
+ self._uninitialized._target = Stream_function(get_coeff,
+ is_sparse=False,
+ approximate_order=approximate_order,
+ true_order=True)
def iterate_coefficients(self):
"""
@@ -1368,56 +1404,40 @@ class Stream_functional_equation(Stream_inexact):
"""
yield from self._initial_values
- from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing
- P = InfinitePolynomialRing(self._base, names=('FESDUMMY',))
- x = P.gen()
- PFF = P.fraction_field()
- offset = self._approximate_order
-
- def get_coeff(n):
- return x[n-offset]
-
- sf = Stream_function(get_coeff, is_sparse=False, approximate_order=offset, true_order=True)
- self._F = self._F.replace(self._uninitialized, sf)
+ x = self._P.gen()
+ PFF = self._P.fraction_field()
n = self._F._approximate_order
- data = list(self._initial_values)
while True:
coeff = self._F[n]
if coeff.parent() is PFF:
coeff = coeff.numerator()
else:
- coeff = P(coeff)
+ coeff = self._P(coeff)
V = coeff.variables()
- # Substitute for known variables
- if V:
- # The infinite polynomial ring is very brittle with substitutions
- # and variable comparisons
- sub = {str(x[i]): val for i, val in enumerate(data)
- if any(str(x[i]) == str(va) for va in V)}
- if sub:
- coeff = coeff.subs(sub)
- V = P(coeff).variables()
-
- if len(V) > 1:
- raise ValueError(f"unable to determine a unique solution in degree {n}")
-
if not V:
if coeff:
raise ValueError(f"no solution in degree {n} as {coeff} != 0")
n += 1
continue
+ if len(V) > 1:
+ raise ValueError(f"unable to determine a unique solution in degree {n}")
+
+ # the following test is wrong
+ if V[0] != x[n + len(self._initial_values)]:
+ print(f"the solutions to the coefficients must be computed in order")
+ # raise ValueError(f"the solutions to the coefficients must be computed in order")
+
# single variable to solve for
hc = coeff.homogeneous_components()
if not set(hc).issubset([0,1]):
raise ValueError(f"unable to determine a unique solution in degree {n}")
- if str(hc[1].lm()) != str(x[len(data)]):
- raise ValueError(f"the solutions to the coefficients must be computed in order")
- val = -hc.get(0, P.zero()).lc() / hc[1].lc()
+
+ val = -hc.get(0, self._P.zero()).lc() / hc[1].lc()
# Update the cache
- data.append(val)
+ self._F.map_cache(lambda c: c.subs({V[0]: val}))
yield val
n += 1
@@ -1461,6 +1481,13 @@ class Stream_unary(Stream_inexact):
self._series = series
super().__init__(is_sparse, true_order)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._series.map_cache(function)
+
def __hash__(self):
"""
Return the hash of ``self``.
@@ -1576,6 +1603,14 @@ class Stream_binary(Stream_inexact):
self._right = right
super().__init__(is_sparse, False)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._left.map_cache(function)
+ self._right.map_cache(function)
+
def __hash__(self):
"""
Return the hash of ``self``.
@@ -2392,6 +2427,13 @@ class Stream_plethysm(Stream_binary):
f = Stream_map_coefficients(f, lambda x: p(x), is_sparse)
super().__init__(f, g, is_sparse)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ super().map_cache(function)
+ self._powers = [g.map_cache(function) for g in self._powers]
+
@lazy_attribute
def _approximate_order(self):
"""
@@ -3273,6 +3315,12 @@ class Stream_shift(Stream):
self._shift = shift
super().__init__(series._true_order)
+ def map_cache(self, function):
+ r"""
+ Apply ``function`` to each element in the cache.
+ """
+ self._series.map_cache(function)
+
@lazy_attribute
def _approximate_order(self):
""" |
I was thinking about this but could not find a solution other than using you Suppose that we do It seems to me that the fundamental problem is as follows:
Do you see another way to achieve this? Or do you think that the I am quite convinced that we want to specialize the placeholders in the cache whenever a coefficient is determined, because of when keeping them undetermined, the size of the equations may grow extremely fast, as for example in the Jacobi functional equation. |
|
Sorry for being slow to respond. I think the issue you are seeing is because you cannot assume that only the largest degree in a series's cache will be the only one you need to substitute at any given time. At the end, that is true, but potentially not for the intermediaries. I cannot think of another way around it other than going through the entire case of everything in the stream-tree. For what's optimal, I would say it likely depends on the situation. Yet, I suspect doing the replace will generally be faster to large coefficients. It also is easier to prove the code is correct IMO. The replace might also be easier to generalize to the multivariate case. I was able to get a significant speedup: By simply substituting in the known values into the We might be able to get some other speedups by controlling ourselves the polynomial rings. It might make some of the code more robust too... However, I think that is best for a followup. I also moved the How fast is the Jacobi example with your fricas implementation? Perhaps this might be a good benchmark. |
|
I finally managed to adapt the idea of substituting the variables in all caches recursively. The result is quite extreme. Computing Moreover, I did not even do any obvious optimizations yet. One thing I do not really know is, whether it could happen that not only the last element of the cache contains the new variable. Thinking of van der Hoeven's algorithm, #34616, with respect to this question makes my head spin. Since this is such a huge improvement, I'd rather not bury it in a separate pull request for which I won't really have the time for. One idea to make what's misnamed |
|
That is certainly quite a big improvement! Thank you for doing it. I have added the patch. For future reference, you could always just make a PR on my fork based off my branch (so it will be included in this PR). Just to check this is to be taken in conjunction with my current approach of a substitution rather than a replacement (ahem...substitution...) for it? If it is meant to be a strict a replacement, then I am not sure how to prove correctness (and I suspect I could find an example where this breaks through some cancellation of higher order terms than the ones computed). |
Thank you! Yes, it was just a proof of concept so far. Since you say you like it, I will improve on it (and make a pull request).
Yes, exactly. As mentioned above (#36233 (comment)) the replace mechanism seems to be necessary.
It would be wonderful if you could cook up an example where, when computing the next term of |
|
I am in the process of adding doc and tests to your patch. Please wait a bit before making any changes. I will let you know when I leave for the night, at which time I won't work more on this for the day. |
|
@mantepse Okay, I finished putting on the polish and bringing stuff to full test coverage. I probably won't do any more work on this for today. I included a remaining of the method, but I am not sure it is much better. Anyways, feel free to decide on a name you are happy with. For the example, I was thinking something about taking a derivative of something then evaluating the result at |
tscrim
force-pushed
the
lazy_series/integration
branch
from
c62ee41 to
d91b508
Compare
|
is this ready? |
|
No, everything related to implicit definitions of power series should be removed, because this will be done in #37033, and it will be done differently. |
tscrim
force-pushed
the
lazy_series/integration
branch
from
2688cb5 to
3ab190b
Compare
|
@mantepse Sorry for taking a while to split the two. Here is the branch without the implicit definition stuff and all commits squashed down to avoid accidentally deleting stuff with your branch. If you remove all of the integration/Taylor series stuff from your branch and squash everything down, then we shouldn't get any conflicts. |
tscrim
changed the title
| interprets it as the number of times to integrate the function. | ||
| If ``variable`` is not the variable of the Laurent series, then | ||
| the coefficients are integrated with respect to ``variable``. | ||
|
|
There was a problem hiding this comment.
I'm afraid that this conflicts with the conventions in https://doc.sagemath.org/html/en/reference/power_series/sage/rings/multi_power_series_ring_element.html#sage.rings.multi_power_series_ring_element.MPowerSeries.integral
I think it is good to have the possibility of providing integration constants, but it is more important to be consistent with multivariate power series. So maybe we could have a keyword only parameter constants?
desolve has other conventions, but from my experience (and I use it quite a bit because of my teaching) they are a nightmare.
There was a problem hiding this comment.
I just realised that this is the integral for univariate series. I am not sure whether my comment still applies.
There was a problem hiding this comment.
The multivariate power series is a clear outlier, and it also doesn't support the same functionality as other cases:
sage: R.<x,y> = QQ['t'][[]]
sage: t = R.base_ring().gen()
sage: x.integral(t) # errors out instead of t*x
From this, I would say the multivariate power series needs to be brought into line.
I will address the rest of the comments later when I have a bit more time.
src/sage/rings/lazy_series_ring.py
Outdated
| sage: L.taylor(f) | ||
| 1 + z + z^2 + z^3 + z^4 + z^5 + z^6 + O(z^7) | ||
|
|
||
| For inputs as symbolic functions, this uses the generic implementation, |
There was a problem hiding this comment.
maybe replace For inputs as symbolic functions with If f is in the symbolic ring
There was a problem hiding this comment.
That actually changes things in a very subtle but meaningful way. In particular, polynomials (e.g., QQ['x,y']) are also in SR but are not symbolic functions. So I think it is better to not change this.
|
I have addressed all of the things I agree with and commented on why I did not make those changes. I also implemented the multivariate Taylor series by using |
|
Documentation preview for this PR (built with commit fe6bac8; changes) is ready! 🎉 |
|
Thank you! (The one doctest failure is almost certainly unrelated.) |
|
I agree and I cannot reproduce locally. Thank you. |
…ries
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes sagemath#1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->
We implement an `integral()` method for lazy Laurent and power series.
If we raise an error if we perform $\int t^{-1} dt$.
We also prove a way to construct a lazy (power) series as the Taylor
series of a function.
<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes sagemath#12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->
### 📝 Checklist
<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->
- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.
### ⌛ Dependencies
<!-- List all open PRs that this PR logically depends on
- sagemath#12345: short description why this is a dependency
- sagemath#34567: ...
-->
- sagemath#35485: Uses the `is_uninitialized()` checking.
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
URL: sagemath#36233
Reported by: Travis Scrimshaw
Reviewer(s): Martin Rubey, Travis Scrimshaw
…ries
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes sagemath#1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->
We implement an `integral()` method for lazy Laurent and power series.
If we raise an error if we perform $\int t^{-1} dt$.
We also prove a way to construct a lazy (power) series as the Taylor
series of a function.
<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes sagemath#12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->
### 📝 Checklist
<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->
- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.
### ⌛ Dependencies
<!-- List all open PRs that this PR logically depends on
- sagemath#12345: short description why this is a dependency
- sagemath#34567: ...
-->
- sagemath#35485: Uses the `is_uninitialized()` checking.
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
URL: sagemath#36233
Reported by: Travis Scrimshaw
Reviewer(s): Martin Rubey, Travis Scrimshaw
…ries
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes sagemath#1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->
We implement an `integral()` method for lazy Laurent and power series.
If we raise an error if we perform $\int t^{-1} dt$.
We also prove a way to construct a lazy (power) series as the Taylor
series of a function.
<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes sagemath#12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->
### 📝 Checklist
<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->
- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.
### ⌛ Dependencies
<!-- List all open PRs that this PR logically depends on
- sagemath#12345: short description why this is a dependency
- sagemath#34567: ...
-->
- sagemath#35485: Uses the `is_uninitialized()` checking.
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
URL: sagemath#36233
Reported by: Travis Scrimshaw
Reviewer(s): Martin Rubey, Travis Scrimshaw
We implement an$\int t^{-1} dt$ .
integral()method for lazy Laurent and power series. If we raise an error if we performWe also prove a way to construct a lazy (power) series as the Taylor series of a function.
📝 Checklist
⌛ Dependencies
is_uninitialized()checking.