sagemath / sage

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

Make Airy functions symbolic #12455

Closed 38968367-17c9-42b1-b82d-c1adf20431c2 closed 9 years ago

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago

As discussed in sage-support.

Currently sage can evaluate airy functions numerically:

sage: airy_ai(1.4)
0.0820380498076

but it doesn't work symbolically

sage: airy_ai(x)
...
TypeError: Cannot evaluate symbolic expression to a numeric value.

We should make it symbolical for both airy_ai and airy_bi, as well as their derivatives.

Depends on #12289 Depends on #17130

CC: @kcrisman @burcin @benjaminfjones @fredrik-johansson @eviatarbach

Component: symbolics

Keywords: Airy functions sd40.5 sd48

Author: Oscar Gerardo Lazo Arjona, Benjamin Jones, Douglas McNeil, Eviatar Bach, Ralf Stephan

Branch: 2f6945a

Reviewer: Eviatar Bach, Karl-Dieter Crisman, Burcin Erocal, Ralf Stephan, Jeroen Demeyer, Marc Mezzarobba

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

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago

Description changed:

--- 
+++ 
@@ -1,3 +1,5 @@
+As discussed in [sage-support](http://groups.google.com/group/sage-support/browse_thread/thread/f458d0f9cfd89c9a).
+
 Currently sage can evaluate airy functions numerically:
38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:4

Attachment: trac_12455-symbolic_airy.patch.gz

I've added a patch, which should do the job, but it has a few shortcomings:

1.-The resulting symbolic functions seem to remain on hold:

sage: airy_ai(1.0)
airy_ai(1.00000000000000)

You need to force it to evaluate:

sage: airy_ai(1.0).n()
0.135292416313

2.- This doesn't work:

sage: airy_ai(2.0).n(digits=100)
0.0349241304233

3.- There is no evaluation for airy_ai_prime or airy_bi_prime

4.- I'm not sure about how should the functions be called, some possible schemes are

{ai,bi,aip,bip} {ai,bai,aip,baip} {airy_ai,airy_bi,airy_ai_prime,airy_bi_prime}

And also whether the latex representation should be capitalized or not. I chose the third scheme, and capitalized typesetting.

kcrisman commented 12 years ago
comment:6

Good start. A few responses.

Thanks for getting us a start on this! That's great.

fredrik-johansson commented 12 years ago
comment:7

mpmath has the derivatives (as well as integrals) -- just use the optional 'derivative' parameter.

kcrisman commented 12 years ago
comment:8

mpmath has the derivatives (as well as integrals) -- just use the optional 'derivative' parameter.

Ok, I didn't realize that was what that parameter was about. Though in retrospect it seems obvious!

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:9

On second thought, I think it would be better to use the airy equation to calculate derivatives or order higher than 1. Like

sage: airy_ai(2,x)
x*airy_ai(x)
sage: airy_ai(3,x)
airy_ai(x)+x*airy_ai_prime(x)
sage: diff(airy_ai(x),x,2)
x*airy_ai(x)
sage: diff(airy_ai(x),x,3)
airy_ai(x)+x*airy_ai_prime(x)

which is very likey to be the way mpmath calculates higher order derivatives. Integrals however, would be returned as:

sage: airy_ai(-1,x)
airy_ai(-1,x)
sage: integral(airy_ai(x),x)
airy_ai(-1,x)

what do you think?

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:10

Attachment: trac_12455-symbolic_airy2.patch.gz

I just added a new patch. The new version includes generalized derivatives, evaluation with mpmath, and special values of the functions and their derivatives. Just a one doubt, the coverage is:

oscar@oscar-netbook:~$ sage -coverage airy.py
SCORE /home/oscar/sage/my_patches/airy.py: 76% (20 of 26)

Missing documentation:
     * __init__(self):
     * __init__(self):
     * __init__(self):
     * __init__(self):
     * __init__(self):
     * __init__(self):

Possibly wrong (function name doesn't occur in doctests):
     * _derivative_(self, x, diff_param=None):
     * _eval_(self, x):
     * _evalf_(self, x, parent=None):
     * _derivative_(self, x, diff_param=None):
     * _eval_(self, x):
     * _evalf_(self, x, parent=None):
     * _derivative_(self, alpha, *args, **kwds):
     * _eval_(self, alpha, *args):
     * _evalf_(self, alpha, x, parent=None):
     * _derivative_(self, x, diff_param=None):
     * _eval_(self, x):
     * _evalf_(self, x, parent=None):
     * _derivative_(self, x, diff_param=None):
     * _eval_(self, x):
     * _evalf_(self, x, parent=None):
     * _derivative_(self, alpha, *args, **kwds):
     * _eval_(self, alpha, *args):
     * _evalf_(self, alpha, x, parent=None):

Is that important?

kcrisman commented 12 years ago
comment:11

Is that important?

The second part isn't, as you are implicitly testing these 'hidden' functions. By the way,

Examples::

should be

EXAMPLES::

However, the first part (the initialization methods) is important for our coverage. You can just do a couple from the big examples that you have.

Do these pass doctests? I'm surprised that

sage: airy_ai_simple= FunctionAiryAiSimple()

would work when you run tests, assuming you didn't import FunctionAiryAiSimple into the global namespace.

Anyway, overall on a quick glance this looks great; unfortunately, I won't have time to review it properly soon - hopefully someone else will, because you put a lot of great work into it and we should totally have these symbolic. Thanks!

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:12

Do these pass doctests? I'm surprised that

sage: airy_ai_simple= FunctionAiryAiSimple()

would work when you run tests, assuming you didn't import FunctionAiryAiSimple into the global namespace.

It looks like they do:

oscar@oscar-netbook:~$ sage -coverage airy.py
sage -t  "/home/oscar/sage/my_patches/airy.py"              
     [51.5 s]

----------------------------------------------------------------------
All tests passed!
Total time for all tests: 51.7 seconds

Why shouldn't they? I'll take care of the __init__ functions.

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:13

I'm sorry, that should be:

oscar@oscar-netbook:~$ sage -t airy.py
sage -t  "/home/oscar/sage/my_patches/airy.py"              
     [51.5 s]

----------------------------------------------------------------------
All tests passed!
Total time for all tests: 51.7 seconds

(the command was sage -t)

kcrisman commented 12 years ago
comment:14
sage: airy_ai_simple= FunctionAiryAiSimple()

would work when you run tests, assuming you didn't import FunctionAiryAiSimple into the global namespace.

I'm just surprised this doesn't cause trouble.

On another note, this apparently does something weird when applied to 5.0.beta4. Namely, Sage won't start:

---> 65 from sage.functions.all import sin, cos

ImportError: cannot import name sin
Error importing ipy_profile_sage - perhaps you should run %upgrade?
WARNING: Loading of ipy_profile_sage failed.

I'll spare you the rest of the traceback. I think there is some kind of circular import thing going on here, but I don't understand imports that well. Perhaps moving the import of airy to further down (after trig functions) would help - that's a totally uninformed guess, of course.

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago
comment:15

Attachment: trac_12455-symbolic_airy3.patch.gz

I just added a new patch, now with 100 % test coverage, but I still get this message:

oscar@oscar-netbook:~$ sage -coverage /home/oscar/sage/my_patches/airy.py ----------------------------------------------------------------------
/home/oscar/sage/my_patches/airy.py
ERROR: Please add a `TestSuite(s).run()` doctest.
SCORE /home/oscar/sage/my_patches/airy.py: 100% (26 of 26)
----------------------------------------------------------------------

What is this TestSuite thing? I also changed the order in which functions are initialized, so that airy_ai_general is initialized first. Doing it last made some symbolic operations, such as simplify not work (It said something about airy_ai requiring 2 arguments).

38968367-17c9-42b1-b82d-c1adf20431c2 commented 12 years ago

Author: Oscar Gerardo Lazo Arjona

benjaminfjones commented 12 years ago

Reviewer: Benjamin Jones

benjaminfjones commented 12 years ago
comment:16

Hi Oscar, I haven't looked over your patch very closely yet, but I plan to. One comment on first glance: in several places you call

return mpmath_utils.call(airy_bi_mpmath, x, derivative=1, parent=RR) 

but I think we want to preserve the parent of x instead of forcing it to be RR. Importing the parent function from sage.structure.coerce using a different name like sage_parent (because as burcin pointed out in one of my tickets, parent is the name of the parameter in _evalf_) and doing something like

from sage.structure.coerce import parent as sage_parent
R = parent or sage_parent(x)
return mpmath_utils.call(airy_bi_mpmath, x, derivative=1, parent=R)

would do the trick I think.

Other comment: can we add conversions to Maxima along with the Mma conversion you included? Here is a link to the appropriate chapter in the manual: http://maxima.sourceforge.net/docs/manual/en/maxima_15.html#SEC80

benjaminfjones commented 12 years ago
comment:17

The patch applies to 5.0.beta6 cleanly, but upon sage -br I get an import error ending with:

/home/jonesbe/sage/sage-5.0.beta6/local/lib/python2.7/site-packages/sage/gsl/dft.py in <module>()
     63 from sage.rings.real_mpfr import RR
     64 from sage.rings.all import I
---> 65 from sage.functions.all import sin, cos
     66 from sage.gsl.fft import FastFourierTransform
     67 from sage.gsl.dwt import WaveletTransform

ImportError: cannot import name sin
Error importing ipy_profile_sage - perhaps you should run %upgrade?
WARNING: Loading of ipy_profile_sage failed.

Seems like the same thing kcrisman ran into with the previous patch applied to 5.0.beta4.

benjaminfjones commented 12 years ago
comment:18

I found the cause of the problem in comment [comment:17]. There is a circular import going on because

from airy import airy_ai, airy_bi

occurs at the top of sage/functions/all.py. If you move it to the bottom of the file, sage starts just fine.

There a several doctests that fail, however. Oscar, I think that running sage -coverage only checks to see if there are doctest for every function, it doesn't actually run the tests. Use

sage -t devel/sage/sage/functions/airy.py

to run the tests and get feedback.

benjaminfjones commented 12 years ago

Work Issues: circular import, doctest failures

benjaminfjones commented 12 years ago
comment:19

Here at SD 40.5, we're going to pick this ticket up, do the work, and get it into sage.

benjaminfjones commented 12 years ago

Changed reviewer from Benjamin Jones to none

benjaminfjones commented 12 years ago

Changed author from Oscar Gerardo Lazo Arjona to Oscar Gerardo Lazo Arjona, Benjamin Jones

benjaminfjones commented 12 years ago

Changed keywords from Airy functions to Airy functions sd40.5

56048686-9665-4c5e-973e-6c3add3aa805 commented 12 years ago
comment:21

Okay, I'm posting this for future reference. I've attempted to preserve as much as possible of the original code while switching to the "new-style" dictionary-argument .n() approach we're working on in #12289, but I'm setting it to "needs work" myself while we think through some of the consequences of the switch. Some of the _evalf_ stuff violates DRY and I'm not yet sure of the best way to deal with it.

As well, the patch currently segfaults in testing when it computes beta(0.5, 0.5) due to some problem with the current #12289 package, which I suspect is unrelated. But I think the new approach is beginning to take shape.

56048686-9665-4c5e-973e-6c3add3aa805 commented 12 years ago

Attachment: trac_12455_newstyle_airy.patch.gz

newstyle rewrite, version 1

56048686-9665-4c5e-973e-6c3add3aa805 commented 12 years ago

Dependencies: #12289

56048686-9665-4c5e-973e-6c3add3aa805 commented 12 years ago
comment:23

Oh, yeah: and #13028 was a significant nuisance while working on this.

eviatarbach commented 11 years ago

Attachment: trac_12455_newstyle_airy2.patch.gz

eviatarbach commented 11 years ago
comment:25

New patch, which should pass all doctests once #12289 is merged.

It makes the following changes:

Patchbot apply trac_12455_newstyle_airy.patch trac_12455_newstyle_airy2.patch

eviatarbach commented 11 years ago

Changed keywords from Airy functions sd40.5 to Airy functions sd40.5 sd48

kcrisman commented 11 years ago
comment:27

Needs review, then? Note that I had to rebase #12289.

kcrisman commented 11 years ago

Description changed:

--- 
+++ 
@@ -35,3 +35,7 @@

We should make it symbolical for both airy_ai and airy_bi, as well as their derivatives. + +--- + +Apply attachment: trac_12455_newstyle_airy.patch and attachment: trac_12455_newstyle_airy2.patch.

kcrisman commented 11 years ago

Changed author from Oscar Gerardo Lazo Arjona, Benjamin Jones to Oscar Gerardo Lazo Arjona, Benjamin Jones, Eviatar Bach

kcrisman commented 11 years ago

Description changed:

--- 
+++ 
@@ -38,4 +38,4 @@

 ---

-Apply [attachment: trac_12455_newstyle_airy.patch](https://github.com/sagemath/sage-prod/files/10654776/trac_12455_newstyle_airy.patch.gz) and [attachment: trac_12455_newstyle_airy2.patch](https://github.com/sagemath/sage-prod/files/10654777/trac_12455_newstyle_airy2.patch.gz).
+Apply [attachment: trac_12455-newstyle-airy-rebase.patch](https://github.com/sagemath/sage-prod/files/10654778/trac_12455-newstyle-airy-rebase.patch.gz) and [attachment: trac_12455_newstyle_airy2.patch](https://github.com/sagemath/sage-prod/files/10654777/trac_12455_newstyle_airy2.patch.gz).
kcrisman commented 11 years ago
comment:28

Patchbot apply trac_12455-newstyle-airy-rebase.patch and trac_12455_newstyle_airy2.patch

I don't know that either of the needs work issues are still there...

kcrisman commented 11 years ago

Reviewer: Eviatar Bach, Karl-Dieter Crisman

kcrisman commented 11 years ago
comment:29

This needs a trivial extra patch to fix


sage -t sage/functions/special.py
**********************************************************************
File "sage/functions/special.py", line 389, in sage.functions.special._init
Failed example:
    spherical_harmonic(3,2,1,2)
Expected:
    15/4*sqrt(7/30)*e^(4*I)*sin(1)^2*cos(1)/sqrt(pi)
Got:
    15/4*sqrt(7/30)*cos(1)*e^(4*I)*sin(1)^2/sqrt(pi)
**********************************************************************

which is really a rebase issue off of #9880, no worries.

kcrisman commented 11 years ago
comment:30

Okay, a little more rebasing for #12289. Patches coming up.

kcrisman commented 11 years ago

Description changed:

--- 
+++ 
@@ -38,4 +38,4 @@

 ---

-Apply [attachment: trac_12455-newstyle-airy-rebase.patch](https://github.com/sagemath/sage-prod/files/10654778/trac_12455-newstyle-airy-rebase.patch.gz) and [attachment: trac_12455_newstyle_airy2.patch](https://github.com/sagemath/sage-prod/files/10654777/trac_12455_newstyle_airy2.patch.gz).
+Apply [attachment: trac_12455-newstyle-airy-rebase.patch](https://github.com/sagemath/sage-prod/files/10654778/trac_12455-newstyle-airy-rebase.patch.gz) and [attachment: trac_12455-newstyle-airy2-rebase.patch](https://github.com/sagemath/sage-prod/files/10654779/trac_12455-newstyle-airy2-rebase.patch.gz).
kcrisman commented 11 years ago
comment:31

Apply trac_12455-newstyle-airy-rebase.patch and trac_12455-newstyle-airy2-rebase.patch

kcrisman commented 11 years ago
comment:32

Eviatar, how much would you say still needs to be reviewed? (In the sense that you have not yet given it positive review.) Just your patch?

(I'm running doctests now.)

kcrisman commented 11 years ago
comment:33

I made some mistake while rebasing - one moment.

eviatarbach commented 11 years ago
comment:34

I think I was fairly thorough, but I haven't positive-reviewed since my patch hasn't been looked at by someone else.

kcrisman commented 11 years ago
comment:35

Attachment: trac_12455-newstyle-airy-rebase.patch.gz

Okay, now I have to fix something in the other one... sigh. Final version coming soon, it does pass tests!

kcrisman commented 11 years ago

Attachment: trac_12455-newstyle-airy2-rebase.patch.gz

kcrisman commented 11 years ago
comment:36

Patchbot, apply trac_12455-newstyle-airy-rebase.patch and trac_12455-newstyle-airy2-rebase.patch

kcrisman commented 11 years ago
comment:37

(for some reason mpmath is able to evaluate numerically with non-integral n)

But this is even in the documentation you wrote!

Return the `\alpha`-th order fractional derivative with respect to `z`.

"Fractional", right? So ... ? Plus, you do not make this change in airy_bi, only airy_ai, for some reason.

Another question:

sage: (plot(airy_bi(x), (x, -10, 5)) +\ 

Are we allowing this kind of continuation in doctests still? I can't remember but seem to recall something about this being a problem.

Otherwise this all seems fine.

eviatarbach commented 11 years ago
comment:38

I didn't write it. But looking at the mpmath documentation, it should work. It didn't work for most values before my patch anyway, but I'll see what I can do.

kcrisman commented 11 years ago
comment:39

I didn't write it.

Good point, sorry! In fact, it's just copied from the mpmath doc. But now you have taken charge :-) Unfortunately, I can't find a lot of independent implementations of this...

burcin commented 11 years ago

Attachment: trac_12455-airy_review.patch.gz

burcin commented 11 years ago

Description changed:

--- 
+++ 
@@ -38,4 +38,8 @@

 ---

-Apply [attachment: trac_12455-newstyle-airy-rebase.patch](https://github.com/sagemath/sage-prod/files/10654778/trac_12455-newstyle-airy-rebase.patch.gz) and [attachment: trac_12455-newstyle-airy2-rebase.patch](https://github.com/sagemath/sage-prod/files/10654779/trac_12455-newstyle-airy2-rebase.patch.gz).
+Apply 
+
+* [attachment: trac_12455-newstyle-airy-rebase.patch](https://github.com/sagemath/sage-prod/files/10654778/trac_12455-newstyle-airy-rebase.patch.gz)
+* [attachment: trac_12455-newstyle-airy2-rebase.patch](https://github.com/sagemath/sage-prod/files/10654779/trac_12455-newstyle-airy2-rebase.patch.gz)
+* [attachment: trac_12455-airy_review.patch]