#4102 closed enhancement (fixed)
make bessel_J symbolic
Reported by: | jwmerrill | Owned by: | burcin |
---|---|---|---|
Priority: | major | Milestone: | sage-5.11 |
Component: | calculus | Keywords: | sd48 |
Cc: | jason, kcrisman, benjaminfjones, dsm, burcin, eviatarbach | Merged in: | sage-5.11.rc0 |
Authors: | Benjamin Jones, Eviatar Bach, Volker Braun | Reviewers: | Karl-Dieter Crisman, Burcin Erocal |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #9880 | Stopgaps: |
Description (last modified by )
The motivation for this is
sage: plot(bessel_J(1, x), (x, 1, 10)) Traceback (most recent call last): ... TypeError: Unable to convert x (='1-1/8*x^2+1/192*x^4-1/9216*x^6+1/737280*x^8-1/88473600*x^10+1/14863564800*x^12-1/3329438515200*x^14+1/958878292377600*x^16+O(x^17)') to real number.
The problem is that special functions, or at least bessel_J
, can't currently be partially evaluated--that is, called with a SymbolicExpression
as an argument. The model of good behavior is polylog
, for which the above method produces a perfectly nice plot
sage: plot(polylog(1,x),(x,.1,.9)) #makes a fine plot
See discussion at http://groups.google.com/group/sage-support/browse_thread/thread/1b985b080ba2369e/7dee9eed953857f5#7dee9eed953857f5
Release manager:
Apply:
Attachments (16)
Change History (104)
comment:1 Changed 11 years ago by
comment:2 Changed 11 years ago by
- Keywords jason added
comment:3 Changed 11 years ago by
- Cc jason added
- Keywords jason removed
comment:4 Changed 8 years ago by
- Cc kcrisman added
- Report Upstream set to N/A
comment:5 Changed 8 years ago by
- Summary changed from Make special functions behave like PrimitiveFunction to make bessel_J symbolic
This ticket description was too broad. We have lots of tickets on fixing symbolic functions, see symbolics/functions for an overview.
See #3426 and #4230 for other tickets related to bessel functions. It would make sense to fix all these together, with a base class that handles the common properties of bessel functions (differentiation), and subclasses that implement numerical evaluation, etc. for each type.
comment:6 Changed 7 years ago by
- Cc benjaminfjones added
- Description modified (diff)
comment:7 Changed 7 years ago by
See also #10636, which I somehow never saw before. This sage-support thread yields an interesting related error:
sage: var('k') k sage: Z = sum( ((-1)^k*(x^(2*k+1))/factorial(2*k+1)),k,0,oo) sage: Z 1/2*sqrt(pi)*sqrt(2)*sqrt(x)*bessel_j(1/2, x) sage: Z(x=3) 1/2*sqrt(pi)*sqrt(2)*sqrt(3)*bessel_j(1/2, 3) sage: Z(x=3).n() --------------------------------------------------------------------------- TypeError Traceback (most recent call last) sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:20822)() TypeError: cannot evaluate symbolic expression numerically sage: besse bessel_I bessel_J bessel_K bessel_Y
Note that we apparently aren't yet converting Maxima's Bessel properly to 'our' uppercase version because of this.
comment:8 follow-ups: ↓ 9 ↓ 11 Changed 7 years ago by
This wouldn't be hard to implement using #11143 as a template, but what to do about names? I guess new names for the symbolic Bessel functions should be chosen and deprecation notices added to the existing Bessel_J
etc.
What new names do people like in the interim?
bessel_J_symbolic
bessel_J_sym
bessel_Js
???
comment:9 in reply to: ↑ 8 ; follow-up: ↓ 10 Changed 7 years ago by
Replying to benjaminfjones:
What new names do people like in the interim?
bessel_J_symbolic
bessel_J_sym
bessel_Js
???
What is wrong with taking over bessel_{I,J,K,Y}
? Why do we need interim names?
Thanks for looking into this.
comment:10 in reply to: ↑ 9 Changed 7 years ago by
What is wrong with taking over
bessel_{I,J,K,Y}
? Why do we need interim names?
I concur.
bessel_J_symbolic
One could use that as the "symbolic" class one and then do the usual foo = Foo_Class()
?
Thanks for looking into this.
Agreed!
comment:11 in reply to: ↑ 8 Changed 7 years ago by
Replying to benjaminfjones:
What new names do people like in the interim?
I agree with Burcin: just use bessel_J
(etc) if possible.
But regardless of the name, I'd love to see this get in -- it will make teaching a PDE class a bit easier. (No more lambda
expressions in the plots...) Thanks for working on this!
comment:12 Changed 7 years ago by
- Cc dsm added
I thought I'd post some work in progress for all the Bessel function fans in the crowd. There is still work to be done, e.g.
- getting symbolic integration via Maxima to work
- adding maxima, scipy, PARI, etc. to the implemented backend systems
- more documentation and doctests
but at this point the patch applies cleanly to 5.4.1 and all the doctests in sage/functions/special.py
pass.
If you are interested, have a look at the doctests included with the Bessel
function for examples and let me know if you see any serious problems (other than the deficits I've listed above).
Thanks!
comment:13 follow-up: ↓ 15 Changed 7 years ago by
Nice! A few comments of the type you solicited:
- Why
typ
and nottype
? Some Python reserved thing, maybe? But it looks like a typ-o to the (quickly reading) end user. - I'd like to be able to plot
f(x) = Bessel(0)
but maybe that doesn't make sense? I guess a variable is necessary... anyway, just throwing it out there. f = maxima(Bessel(typ='K')(x,y))
turns out great, but does it convert back? Likef.derivative('y')
is-(bessel_k(x+1,y)+bessel_k(x-1,y))/2
, but does it then (when put in the_sage_
method) go back to "our" uppercase Bessel functions?- Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.
- At least some of the error messages should be in doctests, maybe the ones with the wrong type and a non-implemented system.
class_attrs['_conversions'] = {}
--- what do we do with this in Maxima, then? Maybe it's better to raise an error; Maxima tends to otherwise just take things as new variables, which could be dangerous.- How many of the currently-deleted doctests do you think would be worth preserving in the long run? Any deprecation needed here?
Anyway, clearly a lot of planning and looks very promising!
comment:14 follow-up: ↓ 16 Changed 7 years ago by
Thanks for the patch. Bessel functions are symbolic with so few lines of code. Amazing. :)
I like the metaclass idea. That never occured to me before as an option for functions with parameters, like the order here. I have a few questions:
- what is the advantage of creating a new symbolic function for each (type, order) pair, instead of having a generic function for each type that takes the order as a parameter? The latter would be similar to the design of the
psi()
function in GiNaC/Pynac.
- wouldn't this approach make life harder if in the future we want to add exact evaluation method for Bessel_J only?
I am also concerned about blowing up the symbolic function registry with these instances. Note that we keep an entry in a C++ list with a pointer to all the possible custom methods, and a Python dictionary with a function -> Pynac id mapping for each subclass of BuiltinFunction
that is instantiated.
BTW, I suggest moving this code to a new file sage/functions/bessel.py
.
comment:15 in reply to: ↑ 13 ; follow-up: ↓ 17 Changed 6 years ago by
Replying to kcrisman:
Thanks for the comments. I finally found some time to get back to this ticket :)
Nice! A few comments of the type you solicited:
- Why
typ
and nottype
? Some Python reserved thing, maybe? But it looks like a typ-o to the (quickly reading) end user.
That was intentional, I was trying not to shadow the builtin name.
- I'd like to be able to plot
f(x) = Bessel(0)
but maybe that doesn't make sense? I guess a variable is necessary... anyway, just throwing it out there.
I think f(x) = Bessel(0,x)
is better, I don't like x
being implicit on the right hand side.
f = maxima(Bessel(typ='K')(x,y))
turns out great, but does it convert back? Likef.derivative('y')
is-(bessel_k(x+1,y)+bessel_k(x-1,y))/2
, but does it then (when put in the_sage_
method) go back to "our" uppercase Bessel functions?
Good question, I haven't tested it out. I think I'm going to rewrite most of the code anyway so I'll keep this in mind.
- Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.
???
- At least some of the error messages should be in doctests, maybe the ones with the wrong type and a non-implemented system.
Good point, I'll make sure to doctest the exceptions.
class_attrs['_conversions'] = {}
--- what do we do with this in Maxima, then? Maybe it's better to raise an error; Maxima tends to otherwise just take things as new variables, which could be dangerous.
I don't know about this. I had this in the case of the single parameter functions like Bessel(1,'K')
that Maxima doesn't have equivalents for (it just has the two parameter functions). I think I'll scrap the idea of having all infinitely many one-parameter Bessel functions registered as their own symbolic functions (see Burcin's comments below).
- How many of the currently-deleted doctests do you think would be worth preserving in the long run? Any deprecation needed here?
Ideally all of them. Some will need to be modified because of the change in numerical back-end.
Anyway, clearly a lot of planning and looks very promising!
comment:16 in reply to: ↑ 14 Changed 6 years ago by
Replying to burcin:
Thanks for looking at it!
Thanks for the patch. Bessel functions are symbolic with so few lines of code. Amazing. :)
I like the metaclass idea. That never occured to me before as an option for functions with parameters, like the order here. I have a few questions:
- what is the advantage of creating a new symbolic function for each (type, order) pair, instead of having a generic function for each type that takes the order as a parameter? The latter would be similar to the design of the
psi()
function in GiNaC/Pynac.
The advantage I had in mind was just code reuse. In hindsight though, this approach makes the code more complicated and less maintainable that it needs to be. I'm going to refactor it into four generic functions in sage/functions/bessel.py
as you suggest.
- wouldn't this approach make life harder if in the future we want to add exact evaluation method for Bessel_J only?
Good point..
I am also concerned about blowing up the symbolic function registry with these instances. Note that we keep an entry in a C++ list with a pointer to all the possible custom methods, and a Python dictionary with a function -> Pynac id mapping for each subclass of
BuiltinFunction
that is instantiated.
Also a good point. This occurred to me, but I didn't think through the consequences very much. I can see it being a problem if a user can inadvertently create a very large number of symbolic functions at runtime by doing something innocuous like:
sage: point([ (k, Bessel(k, 'J')(pi)) for k in range(1000) ]) ... (1000 new symbolic function classes created!)
BTW, I suggest moving this code to a new file
sage/functions/bessel.py
.
Good idea. New patch coming soon...
comment:17 in reply to: ↑ 15 ; follow-up: ↓ 18 Changed 6 years ago by
- Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.
???
Such as this and this. Just for forward-compatibility (instead of the percent business). Problem is that when I looked for ways to get around braces "naturally" occurring in LaTeX, there weren't necessarily a lot of them. (Ways, that is.)
comment:18 in reply to: ↑ 17 Changed 6 years ago by
Replying to kcrisman:
- Maybe Python 3 string formatting? Though I am not sure how to mix that with LaTeX braces.
???
Such as this and this. Just for forward-compatibility (instead of the percent business). Problem is that when I looked for ways to get around braces "naturally" occurring in LaTeX, there weren't necessarily a lot of them. (Ways, that is.)
Oh, right. I just didn't see what in the code you were referring to. I see now.
Anyway, to one of your earlier questions, with the new code I'm about to post the following conversions to and from Maxima work great:
sage: mb = maxima(bessel_I(1,x)) sage: mb.sage() bessel_I(1, x) sage: x,y = var('x,y') sage: f = maxima(Bessel(typ='K')(x,y)) sage: f.derivative('x') %pi*csc(%pi*x)*('diff(bessel_i(-x,y),x,1)-'diff(bessel_i(x,y),x,1))/2-%pi*bessel_k(x,y)*cot(%pi*x) sage: m = f.derivative('x') sage: m.sage() -1/2*(x*bessel_I(-x, y)/y - x*bessel_I(x, y)/y + bessel_I(-x - 1, y) + bessel_I(x - 1, y))*pi*csc(pi*x) - pi*cot(pi*x)*bessel_K(x, y)
comment:19 Changed 6 years ago by
Just as an FYI, this ask.sagemath question wants lots and lots of precision on Bessel Y. Which I think will be provided here via mpmath - just pointing out we should doctest it.
comment:20 Changed 6 years ago by
comment:21 Changed 6 years ago by
Yet more work in progress, added lots of doctests, in particular to show that problems in #4230 and #3426 are fixed by this patch. One feature this patch adds is the ability to solve Bessel's diffeq (using Maxima) and get a usable symbolic solution back (see the doctests in the module level docstring and in the Bessel() function).
Added brief mathematical exposition on the module docstring as well as some usage EXAMPLES.
I think the patch is more or less ready for initial review. All tests withing the new file sage/functions/bessel.py
pass on my machine and sage-5.6.beta3, but there are several doctests failing in other places that reference the Bessel functions. I'll work on patching those up for the next iteration.
comment:22 follow-up: ↓ 25 Changed 6 years ago by
Dumb comments since I don't have time for proper review - and more importantly, there are no obvious horrible things (though I haven't gone in depth with branches yet). If all this really works and there are no typos, I think this will be a VERY nice addition.
unqiue
typoVerfify
typoincreasin
typo- I don't know why, but the math following "It follows from Bessel's..." doesn't look right in the doc (the
:math:
directive is not parsed) - Is \operatorname{bessel\_I}(\alpha, x) standard or should we just just the I sub whatever?
`test_relation()`
needs to be in double back ticks or have a link to the appropriate place in the ref manual- Trac tickets should use
:trac:
syntax - Does
f(x) = Bessel(0): plot(f, (x, 1, 10))
work, or must one useBessel(0)(x)
? The example after that makes it look like maybe not... - Bessel Y and Bessel K need a little filling out - and one of the
:math:
directives doesn't show up at all, the other isn't parsed right for some reason again - I'd personally get rid of the whitespace changes in sage/functions/special.py - unlikely to have an effect, but not really necessary.
- How should we deal with the removal of the
"maxima"
and"pari"
arguments? I'm not sure if it's really feasible to have a deprecation period for that, but we should discuss it. - I suppose we don't need to keep all old tests, but some of them are rather different and might be worth putting in a TESTS section somewhere, just so that we don't have some unexpected regression only they catch.
- The switch to alpha from nu - I would rather not deprecate this either, but in principle someone could have used it as a keyword argument in the past...
comment:23 Changed 6 years ago by
- Reviewers set to Karl-Dieter Crisman
- Status changed from new to needs_review
Here's the only failures I got with 5.7.beta3.
sage -t "devel/sage-main/sage/calculus/desolvers.py" ********************************************************************** File "/Users/.../sage-5.7.beta3/devel/sage-main/sage/calculus/desolvers.py", line 252: sage: desolve(x^2*diff(y,x,x)+x*diff(y,x)+(x^2-4)*y==0,y) Expected: k1*bessel_j(2, x) + k2*bessel_y(2, x) Got: k1*bessel_J(2, x) + k2*bessel_Y(2, x) ********************************************************************** 1 items had failures: 1 of 62 in __main__.example_1 ***Test Failed*** 1 failures. [10.1 s]
Hmm, should we keep the lowercase versions around, or was that actually an error that we never parsed those? Apparently it was pure Maxima output, looking at the old code, so just replace with the actual return value.
sage -t "devel/sage-main/sage/calculus/wester.py" ********************************************************************** File "/Users/.../sage-5.7.beta3/devel/sage-main/sage/calculus/wester.py", line 39: sage: bessel_J (2, 1+I) Expected: 0.0415798869439621 + 0.247397641513306*I Got: bessel_J(2, I + 1) ********************************************************************** 1 items had failures: 1 of 200 in __main__.example_0 ***Test Failed*** 1 failures. [4.4 s]
Easy enough to fix - we could even add the n()
version there. Actually, we probably should, since the test is probably testing for something - we'll have to just read that quick.
sage -t "devel/sage-main/sage/functions/special.py" ********************************************************************** File "/Users/.../sage-5.7.beta3/devel/sage-main/sage/functions/special.py", line 521: sage: n(bessel_J(3,10,"maxima")) Exception raised: Traceback (most recent call last): sage: n(bessel_J(3,10,"maxima")) File "function.pyx", line 354, in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:3976) TypeError: Symbolic function bessel_J takes exactly 2 arguments (3 given) ********************************************************************** File "/Users/.../sage-5.7.beta3/devel/sage-main/sage/functions/special.py", line 536: sage: n(bessel_J(3,10,"maxima")) Exception raised: n(bessel_J(Integer(3),Integer(10),"maxima"))###line 536: sage: n(bessel_J(3,10,"maxima")) File "function.pyx", line 354, in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:3976) TypeError: Symbolic function bessel_J takes exactly 2 arguments (3 given) ********************************************************************** 2 items had failures: 1 of 5 in __main__.example_8 1 of 5 in __main__.example_9 ***Test Failed*** 2 failures. [3.8 s]
Hmm, here is where that potential deprecation I mentioned might come in.
sage -t "devel/sage-main/sage/symbolic/random_tests.py" ********************************************************************** File "/Users/.../sage-5.7.beta3/devel/sage-main/sage/symbolic/random_tests.py", line 236: sage: print "ignore this"; random_expr(50, nvars=3, coeff_generator=CDF.random_element) # random Exception raised: <snip> File "misc.pyx", line 209, in sage.structure.misc.getattr_from_other_class (sage/structure/misc.c:1488) AttributeError: 'sage.rings.complex_interval.ComplexIntervalFieldElement' object has no attribute 'arcsin' **********************************************************************
Apparently the usual craziness with random expression has reached new heights with these extra functions. I suggest we just try a different seed or something. Not every random expression will be meaningful, for instance.
comment:24 Changed 6 years ago by
- Status changed from needs_review to needs_work
- Work issues set to doctest fixes, various minor documentation and other issues
Needs work:
- doctest fixes, various minor documentation and other issues
And needs review on the code itself :-)
But now it's really looking tractable to finish this off.
comment:25 in reply to: ↑ 22 Changed 6 years ago by
- Cc burcin added
Replying to kcrisman:
Thanks for looking over the patch!
Dumb comments since I don't have time for proper review - and more importantly, there are no obvious horrible things (though I haven't gone in depth with branches yet). If all this really works and there are no typos, I think this will be a VERY nice addition.
unqiue
typoVerfify
typoincreasin
typo- I don't know why, but the math following "It follows from Bessel's..." doesn't look right in the doc (the
:math:
directive is not parsed)
Fixed, it was a one space indentation problem.
- Is \operatorname{bessel\_I}(\alpha, x) standard or should we just just the I sub whatever?
I was unsure about this. I think now it makes most sense to just stick with I_\alpha when in a math block.
`test_relation()`
needs to be in double back ticks or have a link to the appropriate place in the ref manual- Trac tickets should use
:trac:
syntax- Does
f(x) = Bessel(0): plot(f, (x, 1, 10))
work, or must one useBessel(0)(x)
? The example after that makes it look like maybe not...
The way the code works you have to specify the variable on the right hand side, e.g.:
sage: f(x) = Bessel(0)(x) sage: f(x) = bessel_I(0, x)
This is because Bessel()
returns lambda functions under the hood. Personally, I prefer having the variable dependence made explicit, but I'm open to suggestions if there are other
ways you can think of to get that functionality.
- Bessel Y and Bessel K need a little filling out - and one of the
:math:
directives doesn't show up at all, the other isn't parsed right for some reason again
I'll work on Bessel Y and K, I realize now that you mention it that I didn't fill in the definitions there :). Anyone with good doctest suggestions for these is welcome to chime in, I'll add the tests if you can think of interesting or relevant ones.
- I'd personally get rid of the whitespace changes in sage/functions/special.py - unlikely to have an effect, but not really necessary.
Yep, you're right.
- How should we deal with the removal of the
"maxima"
and"pari"
arguments? I'm not sure if it's really feasible to have a deprecation period for that, but we should discuss it.
This is a question that's been troubling me. I don't know how to make the deprecation happen and make progress here without renaming the new functions and then swapping them in when the deprecation period is over. The infrastructure we have for adding symbolic functions (inheriting from BuiltinFunction
, etc) doesn't support any kind of system arguments without overriding call()
.
At least existing code that uses the Bessel functions as they are now in Sage won't break unless they explicitly use the system argument.
Maybe @burcin has a suggestion here?
- I suppose we don't need to keep all old tests, but some of them are rather different and might be worth putting in a TESTS section somewhere, just so that we don't have some unexpected regression only they catch.
Good point, I'll see what I can salvage.
- The switch to alpha from nu - I would rather not deprecate this either, but in principle someone could have used it as a keyword argument in the past...
It doesn't make a difference to me, I used alpha because that's what's on the wikipedia page :) Abramowitz & Stegun uses n or \nu and so does the mpmath documentation. I can easily switch \alpha to \nu with some editor-fu.
comment:26 Changed 6 years ago by
- Status changed from needs_work to needs_review
- Work issues doctest fixes, various minor documentation and other issues deleted
I've uploaded two new patches addressing comments by @kcrisman. The two latest patches are finally ready for complete review.
Patchbot:
Apply trac_symbolic_bessel_v5.patch trac_symbolic_bessel_doctests.patch
comment:27 Changed 6 years ago by
Some additional notes regarding your comments, @kcrisman:
- The doctest of
desolve
that was failing was returning unconverted Maxima output, that's what thebessel_j
, etc are. With the v5 patch we now have an honest Sage symbolic function as the return value in that case (also, see the diff eq doctest examples I wrote).
- I changes from using \alpha as the order to \nu for consistency with mpmath and Maxima.
- I changed the random expression seed to 53 (my favorite random integer) and that test no longer raises an exception.
comment:28 Changed 6 years ago by
Nice work. I do feel like we should ask on sage-devel about deprecating the evaluation system keywords versus (as in the current version) simply removing them but leaving room for their return...
assert _system == 'mpmath'
I suppose we could at least doctest this as a todo with the error message it receives with the "wrong" input.
I still want to give it a final go-over, but various testing of your tests and code makes me think this is ready. What a job!
comment:29 Changed 6 years ago by
I spent some time today thinking about the deprecation issue. Here's my summary:
The current patches introduce two API changes. First, the new bessel_I, bessel_J, etc
functions take two positional arguments whereas the old ones take 2 positional arguments and two optional keyword arguments (algorithm
and prec
). The second API change is the same change in arguments, but to the constructor Bessel
.
I can add deprecation warnings to the Bessel
constructor easily and have it call the old bessel_?
functions during the deprecation period. On the other hand, I don't know how to implement deprecation for the old bessel_?
functions. I can imagine trying to override BuiltinFunction
's call method, or turning the new bessel_?
functions into wrappers which call the new ones if two arguments are used, and give a deprecation warning and call the old versions if more than 2 arguments are given.
I can:
- try to implement the deprecation
- ask on sage-devel for a waiver from the deprecation policy in this case
- other?
What say you all?
comment:30 Changed 6 years ago by
After chatting with @burcin, I decided to do 1. The new patch implements deprecation by retaining all the old code with names prefixed by an underscore. I preserved all the old doctests (which reference the underscored names). In the new code, I overrode __call__
if more than 2 arguments are given.
Ready for review!
comment:31 Changed 6 years ago by
Patchbot, apply trac_symbolic_bessel_v5.patch trac_symbolic_bessel_doctests.patch trac_symbolic_bessel_deprecation.patch
comment:32 Changed 6 years ago by
Can you remind me what the plan was for algorithm as a keyword? Wasn't there some plan to reintroduce it once we got the hooks in properly, so that people could easily compare e.g. Pari, mpmath, Mma (if available), etc.? Doesn't mean this patch couldn't go in, just trying to get things straight.
Quick glance at the latest patch looks good. What a lot of underscores; did you write a script?
comment:33 Changed 6 years ago by
The algorithm keyword in _eval_ requires a work in progress branch of pynac that has probably bit rot by now. I have it on my agenda (since last year!) to look back into it..
comment:34 Changed 6 years ago by
Right, I understand - my point is that maybe the deprecation message should be slightly different, something like
deprecation(4102, 'precision argument is deprecated\nalgorithm argument will only be available as a named keyword in the future and is currently in mothballs')
or something somewhat more professional/accurate than that.
comment:35 Changed 6 years ago by
I see, sure, how about just:
"precision argument is deprecated, algorithm argument is currently deprecated, but will be available as a named keyword in the future"
comment:36 Changed 6 years ago by
I see, sure, how about just: "precision argument is deprecated, algorithm argument is currently deprecated, but will be available as a named keyword in the future"
That sounds good, except that I would put a semicolon instead of a comma in the first comma spot.
comment:37 Changed 6 years ago by
Latest attachment makes a tiny doc change to fix a latex mistake.
comment:38 Changed 6 years ago by
Patchbot apply trac_symbolic_bessel_v5.patch trac_symbolic_bessel_doctests.patch trac_symbolic_bessel_deprecation.patch
comment:39 follow-up: ↓ 40 Changed 6 years ago by
Patchbot shows test failures on sage-5.9.beta2
comment:40 in reply to: ↑ 39 Changed 6 years ago by
Patchbot shows test failures on sage-5.9.beta2
In particular, they are relevant:
sage -t sage/functions/bessel.py ********************************************************************** File "sage/functions/bessel.py", line 977, in sage.functions.bessel.Bessel Failed example: f = desolve(diffeq, y, [1, 1, 1]); f Expected: (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) - (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) Got: (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_j(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1)) - (bessel_j(0, 1) + bessel_j(1, 1))*bessel_Y(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1)) ********************************************************************** File "sage/functions/bessel.py", line 983, in sage.functions.bessel.Bessel Failed example: fp.subs(x=1).n() Exception raised: Traceback (most recent call last): File "/mnt/storage2TB/patchbot/Sage/sage-5.9.beta2/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 460, in _run self.execute(example, compiled, test.globs) File "/mnt/storage2TB/patchbot/Sage/sage-5.9.beta2/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 819, in execute exec compiled in globs File "<doctest sage.functions.bessel.Bessel[29]>", line 1, in <module> fp.subs(x=Integer(1)).n() File "expression.pyx", line 4381, in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:21011) TypeError: cannot evaluate symbolic expression numerically ********************************************************************** 1 item had failures: 2 of 46 in sage.functions.bessel.Bessel [258 tests, 2 failures, 16.6 s] ---------------------------------------------------------------------- sage -t sage/functions/bessel.py # 2 doctests failed ----------------------------------------------------------------------
I don't like that some are lowercase and some uppercase. I think that the second error is just a result of that - the derivative won't function properly.
comment:41 Changed 6 years ago by
There is something strange going on. Here is sage-5.9.beta2 with the three patches applied:
[bjones@cabbage:devel/sage]% ../../sage ---------------------------------------------------------------------- | Sage Version 5.9.beta2, Release Date: 2013-03-28 | | Type "notebook()" for the browser-based notebook interface. | | Type "help()" for help. | ---------------------------------------------------------------------- ********************************************************************** * * * Warning: this is a prerelease version, and it may be unstable. * * * ********************************************************************** sage: y = function('y', x) sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 sage: f = desolve(diffeq, y, [1, 1, 1]); f (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))
which is exactly what the doctest expects.
Now doctesting sage/functions/bessel.py
with sage in the exact same state:
[bjones@cabbage:devel/sage]% ../../sage -t sage/functions/bessel.py Running doctests with ID 2013-04-01-21-31-17-4bc246be. Doctesting 1 file. sage -t sage/functions/bessel.py ********************************************************************** File "sage/functions/bessel.py", line 977, in sage.functions.bessel.Bessel Failed example: f = desolve(diffeq, y, [1, 1, 1]); f Expected: (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) - (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) Got: (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_j(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1)) - (bessel_j(0, 1) + bessel_j(1, 1))*bessel_Y(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1))
The lower cased bessel_j
is what Maxima returns if you don't register a conversion to bessel_J
, the Sage symbolic version of J that these patches introduce.
What is also strange is that the above output does actually represent a different form of the correct solution (modulo replacing bessel_j by bessel_J).
So what is going on here?
FWIW, sage-5.8 with same three patches applied does not exhibit this behavior.
comment:42 Changed 6 years ago by
Maxima is the same version in sage-5.8 and 5.9.beta2, so thats not it.
The new doctesting framework was introduced in sage-5.9.beta, thats a likely suspect. It uses fork, so external interfaces need to be shut down / restarted. Perhaps the restarted maxima process is not set up correctly? Though just using @fork
seems to work.
comment:43 Changed 6 years ago by
benjaminfjones: I don't understand the problem. It seems to me that the doctest should simply be changed to match the "correct" output.
comment:44 Changed 6 years ago by
So the correct output is:
(bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) - (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1))
that is both mathematically correct and what you get when you apply the patches, run sage, and execute the commands manually.
The problem is that when you doctest the same exact commands in the same exact sage, the output that comes back is different in two ways (explained in the comment above).
comment:45 follow-up: ↓ 46 Changed 6 years ago by
More clues: (?)
Here are the three commands in question
sage: y = function('y', x) sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 sage: f = desolve(diffeq, y, [1, 1, 1]); f
- running the commands in an interactive sage session directly after startup, output is:
(bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))
- running the same commands as included in a doctest in
sage/functions/bessel.py
, output is:
(bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_j(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1)) - (bessel_j(0, 1) + bessel_j(1, 1))*bessel_Y(0, x)/(bessel_j(0, 1)*bessel_Y(1, 1) - bessel_j(1, 1)*bessel_Y(0, 1))
- putting the three commands in a docstring alone in a new file
foo.py
and doctesting that, output is:
(bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))
note: same output as in (1)
So it looks like there is some hidden state affecting the doctesting in scenario (2).
comment:46 in reply to: ↑ 45 Changed 6 years ago by
Replying to benjaminfjones:
So it looks like there is some hidden state affecting the doctesting in scenario (2).
The doctests which are run before that test are extra "state".
comment:47 Changed 6 years ago by
See also #2516.
comment:48 Changed 6 years ago by
- Description modified (diff)
I still have not been able to track down the doctesting issue (described starting at comment:41). I've tried bisecting the doctest state (delta debugging) to no avail. The problem seems to be intimately linked to the doctesting environment. I'm unable to reproduce the doctest failure when running all the prior doctests which come before in bessel.py
, in the right order, by hand, and then running the offending one.
The only thing I can tell for sure is that at some point as the doctests are run, one of the pynac symbol tables get's modified. This causes the bessel_j
vs bessel_J
problem.
My inclination is to remove the offending doctest and file a new ticket to resolve the problem. Meanwhile, the bessel function code can be reviewed and become useful. I've posted a patch trac_symbolic_bessel_v6
which removes the offending doctest block.
Comments?
Patchbot, apply trac_symbolic_bessel_v6.patch trac_symbolic_bessel_doctests.patch trac_symbolic_bessel_deprecation.patch
comment:49 follow-up: ↓ 50 Changed 6 years ago by
Tests are passing on 5.10.beta4. The startup time plugin is failing, however. Should we lazy import here? Also how about lazy import for all the special functions (the idea being that relatively few users will need any given special function)?
comment:50 in reply to: ↑ 49 Changed 6 years ago by
Replying to benjaminfjones:
Tests are passing on 5.10.beta4. The startup time plugin is failing, however. Should we lazy import here? Also how about lazy import for all the special functions (the idea being that relatively few users will need any given special function)?
+100 to lazy import for all special functions, but in a separate ticket.
It would be great if you can use lazy import for the functions defined here to silence the startup time plugin.
About the doctesting framework issue: I briefly tried to debug this and didn't get anywhere either. I agree that we should move this to a new ticket and not hold this one up longer. Instead of removing the test, can you mark it # not tested
and mention the ticket number?
Thanks a lot for your work on this.
comment:51 Changed 6 years ago by
I'm trying to understand how the interfacing works. Is this the intended output?
sage: maxima(bessel_I).sage() bessel_i
comment:52 Changed 6 years ago by
maxima_function()
, used in sage.functions.bessel._bessel_J
, overwrites the symbol table entry:
sage: from sage.functions.special import maxima_function sage: type(sage.symbolic.pynac.symbol_table['maxima']['bessel_j']) <class 'sage.functions.bessel.Function_Bessel_J'> sage: maxima_function('bessel_j') bessel_j sage: type(sage.symbolic.pynac.symbol_table['maxima']['bessel_j']) <class 'sage.functions.special.NewMaximaFunction'>
comment:53 follow-up: ↓ 54 Changed 6 years ago by
@burcin: Okay! Good idea. I'll lazy import the Bessel stuff and create a new ticket for the rest.
@vbraun: Thanks! That's good detective work. How did you figure that out (I wasted a lot of time trying to find the culprit)?
I'll look into addressing the symbol table issue and report back with updated patches.
comment:54 in reply to: ↑ 53 Changed 6 years ago by
Replying to benjaminfjones:
@vbraun: Thanks! That's good detective work. How did you figure that out
You better sit down for this one ;-)
$ grep bessel_j sage/functions/*
comment:55 Changed 6 years ago by
Damn.
comment:56 Changed 6 years ago by
- Cc eviatarbach added
comment:57 Changed 6 years ago by
- Status changed from needs_review to needs_work
Just updating status due to the last few comments. This would be really good to have ready to review at Sage Days 48/Edu Days 5 next week.
comment:58 Changed 6 years ago by
This would be nice for my GSoC project as well.
comment:59 Changed 6 years ago by
OK, I'll try to set aside some time for this tomorrow (or weekend).
comment:60 Changed 6 years ago by
Thank you! That would be great :)
comment:61 Changed 6 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
OK, uploaded two new patches which:
- address the problem raised in comment:41, solution is to replace construction of a new
MaximaFunction
object (which alters the symbol table) to construcing and evaluating the function directly inside Maxima usingmaxima.function()
. - change Bessel imports in
sage/functions/all.py
to lazy imports
Doctests in all the touched files pass, I'll wait for the patchbot to see about the rest.
Hope y'all at Sage Days have a chance to review the patches.
Patchbot, apply trac_symbolic_bessel_v7.patch trac_symbolic_bessel_v7-doctests.patch
comment:62 Changed 6 years ago by
- Status changed from needs_review to needs_work
There is another problem in sage/calculus/desolvers.py
, I'm looking into it now.
comment:63 Changed 6 years ago by
When a symbolic function is lazily imported, its constructor isn't called until the function is referenced and the constructor is responsible for populating the symbol tables that make conversion work (e.g. maxima <-> sage). So with lazy importing here, if I start up Sage and call desolve
with Bessel's diffeq, we don't properly convert the result from Maxima to Sage because the Bessel functions are not in the symbol table (they haven't been constructed yet).
Any thoughts on this?
One idea I had is that we can register the conversions somewhere else (not in the constructors). But that would mean maintaining each symbolic function's code and its conversions in separate places (seems like a bad idea).
Another idea: import from sage.functions.bessel strictly (i.e. not lazily), but in bessel.py use lazy imports for everything that's not needed in the constructors.
comment:64 Changed 6 years ago by
I think my preference would be to put 'lazy import symbolic functions' off to another ticket. We can discuss the issues involved there. Meanwhile, this can be reviewed.
Changed 6 years ago by
comment:65 Changed 6 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I removed the lazy import.
Ready for review!
Patchbot apply trac_symbolic_bessel_v7.2.patch trac_symbolic_bessel_v7-doctests.patch
comment:66 Changed 6 years ago by
Thanks, I'm looking at it now.
Can I just remove all the assertions in Bessel
? Some of them are redundant and I'm not sure why it's asserting that the order is real.
comment:67 Changed 6 years ago by
Here's what I changed:
- Added documentation to the manuals
- Fixed a
bessel_K
identity that was incorrect - Added SymPy? conversions
- Minor formatting changes
- Removed assertions (I'll add them back if they were necessary)
- "for an arbitrary real number
\nu
(the order)" > "for an arbitrary complex number\nu
(the order)"
Looks good otherwise!
Changed 6 years ago by
comment:68 Changed 6 years ago by
Thanks, those changes all look fine. The assertions aren't necessary, I should have removed them.
Patchbot apply trac_symbolic_bessel_v7.2.patch trac_symbolic_bessel_v7-doctests.patch bessel_2.patch
comment:69 Changed 6 years ago by
Great. I just noticed that it was already in the manuals though. New patch.
Patchbot apply trac_symbolic_bessel_v7.2.patch trac_symbolic_bessel_v7-doctests.patch bessel_2.2.patch
Changed 6 years ago by
comment:70 Changed 6 years ago by
- Description modified (diff)
- Keywords sd48 added
- Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Burcin Erocal
All patches look good to me. We can switch this to positive review once the patchbot gives it a green light.
comment:71 Changed 6 years ago by
Nice work! A lot of the thing Eviatar fixed were exactly the things I was thinking about on the train. Good stuff. Patchbot doesn't seem to be using the right patch, so:
Patchbot apply trac_symbolic_bessel_v7.2.patch trac_symbolic_bessel_v7-doctests.patch bessel_2.2.patch
comment:72 Changed 6 years ago by
- Status changed from needs_review to needs_work
This passes all tests for me.
comment:73 Changed 6 years ago by
- Status changed from needs_work to positive_review
comment:74 Changed 6 years ago by
- Status changed from positive_review to needs_work
I get several doctest failures:
sage -t --long devel/sage/sage/functions/bessel.py ********************************************************************** File "devel/sage/sage/functions/bessel.py", line 101, in sage.functions.bessel Failed example: bessel_J(0, x).diff(x) Expected: 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x) Got: -1/2*bessel_J(1, x) + 1/2*bessel_J(-1, x) ********************************************************************** File "devel/sage/sage/functions/bessel.py", line 248, in sage.functions.bessel.Function_Bessel_J Failed example: f.diff(x) Expected: 1/2*bessel_J(1, x) - 1/2*bessel_J(3, x) Got: -1/2*bessel_J(3, x) + 1/2*bessel_J(1, x) ********************************************************************** File "devel/sage/sage/functions/bessel.py", line 357, in sage.functions.bessel.Function_Bessel_J._derivative_ Failed example: derivative(f, z) Expected: z |--> 1/2*bessel_J(9, z) - 1/2*bessel_J(11, z) Got: z |--> -1/2*bessel_J(11, z) + 1/2*bessel_J(9, z) **********************************************************************
(and various similar failures)
Changed 6 years ago by
comment:76 Changed 6 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I uploaded a new patch to fix doctest failures caused by #9880.
comment:77 Changed 6 years ago by
- Status changed from needs_review to positive_review
Looks good and passes tests.
comment:78 Changed 6 years ago by
- Status changed from positive_review to needs_work
On silius
(ia64):
sage -t --long devel/sage/sage/functions/bessel.py ********************************************************************** File "devel/sage/sage/functions/bessel.py", line 258, in sage.functions.bessel.Function_Bessel_J Failed example: A[0] Expected: 0.44005058574493355 Got: 0.44005058574493366 **********************************************************************
Changed 6 years ago by
comment:79 Changed 6 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
New patch, adding a tolerance for the integral which is higher than the maximum error given by GSL.
comment:81 Changed 6 years ago by
- Merged in set to sage-5.11.rc0
- Resolution set to fixed
- Status changed from positive_review to closed
comment:82 Changed 6 years ago by
This is wrong:
sage: var('nu z') (nu, z) sage: bessel_J(nu, z).diff(nu) -1/2*bessel_J(nu + 1, z) + 1/2*bessel_J(nu - 1, z) sage: bessel_J(nu, z).diff(z) -1/2*bessel_J(nu + 1, z) + 1/2*bessel_J(nu - 1, z)
Changed 6 years ago by
comment:83 follow-ups: ↓ 84 ↓ 86 Changed 6 years ago by
New patch fixes this issue. Can this be merged in 5.11?
comment:84 in reply to: ↑ 83 Changed 6 years ago by
Replying to eviatarbach:
New patch fixes this issue. Can this be merged in 5.11?
No, you should make a follow-up ticket for this.
comment:85 Changed 6 years ago by
Oh okay, then can this ticket be removed from 5.11? I'm just worried about having mathematically incorrect answers in the release.
comment:86 in reply to: ↑ 83 Changed 6 years ago by
Replying to eviatarbach:
New patch fixes this issue. Can this be merged in 5.11?
Sorry, I really meant: perhaps yes it can be in sage-5.11, but in any case there needs to be a follow-up ticket (make it milestone: sage-5.11 and priority: blocker) and the new patch needs to be reviewed.
comment:87 Changed 6 years ago by
Also, add a doctest for the NotImplementedError
(like the tests from 82)
comment:88 Changed 6 years ago by
Okay, thank you. This is now #15019.
Orthogonal polynomials seem to work fine; you can plot the
legendre_P
,hermite
, andjacobi_P
polynomials withplot(legendre_P(4, x), (x,-1,1))
and so on. So whatever magic they are using might work for the Bessel functions and other special functions.