search

sagemath / sage

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

Sage is missing the lambert_w function #11888

Closed benjaminfjones closed 12 years ago

benjaminfjones commented 12 years ago

Maxima returns solutions to some exponential equations in terms of the lambert_w function. Sage is missing a conversion for this function:


sage: solve(e^(5*x)+x==0, x, to_poly_solve=True)
[x == -1/5*lambert_w(5)]
sage: S = solve(e^(5*x)+x==0, x, to_poly_solve=True)
sage: z = S[0].rhs()
sage: z
-1/5*lambert_w(5)
sage: N(z)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/Users/jonesbe/sage/sage-4.7.2.alpha2/devel/sage-test/sage/<ipython console> in <module>()

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/misc/functional.pyc in numerical_approx(x, prec, digits)
   1264             prec = int((digits+1) * 3.32192) + 1
   1265     try:
-> 1266         return x._numerical_approx(prec)
   1267     except AttributeError:
   1268         from sage.rings.complex_double import is_ComplexDoubleElement

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:17950)()

TypeError: cannot evaluate symbolic expression numerically
sage: lambert_w(5)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)

/Users/jonesbe/sage/sage-4.7.2.alpha2/devel/sage-test/sage/<ipython console> in <module>()

NameError: name 'lambert_w' is not defined
sage:

mpmath can evaluate the lambert_w function, so it should be easy to add a new symbolic function to Sage that will fix this issue.

Apply:

CC: @kcrisman @sagetrac-ktkohl

Component: symbolics

Keywords: lambert_w symbolics conversion maxima sd35.5 sd40.5

Author: Benjamin Jones

Reviewer: Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson, Burcin Erocal, Douglas McNeil, William Stein

Merged: sage-5.1.beta4

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

benjaminfjones commented 12 years ago

Attachment: trac_11888.patch.gz

add lambert_w symbolic function

benjaminfjones commented 12 years ago

Changed keywords from lambert_w symbolics conversion maxima to lambert_w symbolics conversion maxima sd35.5

benjaminfjones commented 12 years ago

Author: Benjamin Jones

benjaminfjones commented 12 years ago
comment:2

Preliminary patch needs review. The function has been added using the template developed as part of #11143. The issue described in the description is addressed in one of the doctests.

kini commented 12 years ago

apply to $SAGE_ROOT/devel/sage

kini commented 12 years ago
comment:3

Attachment: trac_11888-doctests.patch.gz

Running make ptestlong now. I fixed a couple of doctests that broke, and fixed some typos and rST syntax problems in your docstring.

kini commented 12 years ago

Description changed:

--- 
+++ 
@@ -35,3 +35,10 @@

mpmath can evaluate the lambert_w function, so it should be easy to add a new symbolic function to Sage that will fix this issue. + +--- + +Apply: + +1. attachment: trac_11888.patch to $SAGE_ROOT/devel/sage +2. attachment: trac_11888-doctests.patch to $SAGE_ROOT/devel/sage

kini commented 12 years ago
comment:4

All tests pass.

benjaminfjones commented 12 years ago
comment:6

Thanks for the fixes, kini. I've run make ptestlong with the patches applied and verified that all tests pass. Maybe I can get @kcrisman to finish a review this afternoon.

kcrisman commented 12 years ago
comment:7

I don't see any obvious problems, but the random expression usually doesn't change much with these new functions and this one is really different.

It's also spread across many lines, and I'm not sure if this is appropriate (just in this one case, of course).

kini commented 12 years ago
comment:8

I spread it across lines because 1) I was trying to keep within the recommended PEP 8 guidelines for line length, and 2) because of this

[2012-01-10 22:54:53] <kini> while I was fixing the second doctest, some weird
stuff started happening to vim
[2012-01-10 22:55:02] <kini> I thought my terminal had frozen or something
[2012-01-10 22:56:02] <kini> but it turns out that apparently opening a new line
after a line with a 1800-character-long Sage symbolic expression on it causes
vim to take a full 12 seconds to compute the correct indentation level for the
next line
[2012-01-10 22:56:20] <benjaminfjones> ha!
[2012-01-10 22:56:30] <kini> on a 4.5 GHz Core i5-2500K and utilizing three
cores!
[2012-01-10 22:56:39] <benjaminfjones> wow

What is inappropriate about adding line breaks?

As for the length of the expression, it seems to be a fluke. With the patches applied, starting with random seeds other than 2 gives expressions of a more "normal" length.

benjaminfjones commented 12 years ago
comment:9

I agree, it looks like a fluke that the expression grows so large. I did some testing of random_expr and found that it "normally" produces output around 200 - 400 character long, but occasionally the outputs can be 10 times that (I saw a few around 2500 characters long!)

fredrik-johansson commented 12 years ago
comment:10

I strongly recommend implementing the general version of the Lambert W function (taking a branch parameter).

kcrisman commented 12 years ago
comment:11

I strongly recommend implementing the general version of the Lambert W function (taking a branch parameter).

Can you be more specific? (Is this standard with other multivalued functions in Sage?) Maybe this could be a separate ticket, unless the change was really easy.

benjaminfjones commented 12 years ago
comment:12

The change should be simple. mpmath implements the a branch W_k(z) for each integer k. It's just a matter of adding a second parameter to the wrapper and putting in some tests. I'm sitting on the train from Beverly MA to Logan airport now, I'll see if I can get it uploaded before the train stops (or my battery dies).

kcrisman commented 12 years ago
comment:13

Sweet, I didn't realize it was that quick. I love doing Sage development on that train :) There is also free wifi at Logan, I believe.

kcrisman commented 12 years ago
comment:14

Ping. I'd love to review this but sounds like Fredrik's point is good and if it's pretty easy for you to add that, we might as well.

fredrik-johansson commented 12 years ago
comment:15

Yes, it should be easy; just add an optional branch parameter, lambertw(z, branch=0).

Another suggestion is to use scipy.special.lambertw for evaluation over RDF and CDF. The SciPy implementation is a Cython translation of the double precision version in mpmath; it supports all branches and has excellent numerical stability, and runs quite a bit faster.

import scipy.special
import mpmath
timeit("mpmath.lambertw(-35.0r+4.6jr,2r)")
timeit("mpmath.fp.lambertw(-35.0r+4.6jr,2r)")
timeit("scipy.special.lambertw(-35.0r+4.6jr,2r)")
print repr(complex(mpmath.lambertw(-35.0r+4.6jr,2r)))
print repr(mpmath.fp.lambertw(-35.0r+4.6jr,2r))
print repr(scipy.special.lambertw(-35.0r+4.6jr,2r))

625 loops, best of 3: 301 µs per loop
625 loops, best of 3: 65.1 µs per loop
625 loops, best of 3: 6.75 µs per loop
(0.91763023745202721+14.071606637742889j)
(0.91763023745202721+14.071606637742889j)
(0.91763023745202721+14.071606637742889j)
kcrisman commented 12 years ago
comment:16

Nice; I wonder if there are more places we are beginning to default to mpmath where SciPy could be useful for the double fields.

kcrisman commented 12 years ago

Work Issues: add second parameter, RDF/CDF stuff

kcrisman commented 12 years ago

Reviewer: Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson

benjaminfjones commented 12 years ago
comment:17

Thanks for the ping. I'm still here (and I have a patch pretty much ready to go) I just got buried under teaching. I'll try to upload a patch this evening.

benjaminfjones commented 12 years ago
comment:18

After looking at this a bit more, it may be more involved than I initially envisioned to implement arbitrary branches of the lambert_w function in one symbolic function. Right now, the patch from SD 35.5 implements a subclass of !BuiltinFunction. The underlying assumptions about subclasses of BuiltinFunction include: (from sage/symbolic/function.pyx)

We assume that each subclass of this class will define one symbolic function.

One issue is that there isn't a way (as far as I can see) to pass a branch parameter to !BuiltinFunction's _call_ method. (Perhaps burcin or other authority on Sage symbolics can comment on this.)

Changing the evaluation numerical eval to use SciPy would be an easy change, that's for sure. I can do that quickly and upload a patch that implements the principle branch only.

Another idea I just had was to do something like what we have for the Bessel functions, in particular the !Bessel class in sage/functions/special.py which is just a basic python class returning one of the Bessel (I,J,Y) functions of a given order.

kcrisman commented 12 years ago
comment:19

Ok, that makes sense. I feel like there should be a way to do that nonetheless - see incomplete_gamma, with

        BuiltinFunction.__init__(self, "gamma", nargs=2, latex_name=r"\Gamma",
                conversions={'maxima':'gamma_incomplete', 'mathematica':'Gamma',
                    'maple':'GAMMA'})

and then use _eval_ and _evalf_, but I don't have time to try looking into whether that would work here now.

Based on Fredrik's comment, make sure to only use SciPy for RDF/CDF - hopefully there is a good model elsewhere to use for that.

benjaminfjones commented 12 years ago

proof of concept patch (not ready for review)

benjaminfjones commented 12 years ago
comment:20

Attachment: trac_11888_v2.patch.gz

OK, I made a second attempt. The patch isn't complete (I need to fix and add docstrings and do more testing) and is not ready for review, but if the reviewers will take a look at the basic implementation and give me feedback, I'd appreciate it.

In attachment: trac_11888_v2.patch there is a new symbolic function lambert_w_branch which takes two arguments, a complex number z and an integer branch n. This is implemented using scipy.special.lambertw for RDF/CDF arguments z and using mpmath otherwise.

There is also a wrapper function lambert_w that accepts either one or two arguments. For one argument it returns the principle branch lambert_w_branch(z,0), for two it returns lambert_w_branch(z,n). I still need to add the conversion from Maxima (by hand now, since lambert_w doesn't inherit from BuiltinFunction any more).

benjaminfjones commented 12 years ago
comment:21

I fixed the doctests and added lambert_w to the symbol table. I verified that all tests pass including the random_tests.py ones. The patch attachment: trac_11888_v3.patch is ready for review.

benjaminfjones commented 12 years ago

Changed work issues from add second parameter, RDF/CDF stuff to none

benjaminfjones commented 12 years ago

Description changed:

--- 
+++ 
@@ -40,5 +40,4 @@

 Apply:

-1. [attachment: trac_11888.patch](https://github.com/sagemath/sage-prod/files/10653796/trac_11888.patch.gz) to `$SAGE_ROOT/devel/sage`
-2. [attachment: trac_11888-doctests.patch](https://github.com/sagemath/sage-prod/files/10653797/trac_11888-doctests.patch.gz) to `$SAGE_ROOT/devel/sage`
+1. [attachment: trac_11888_v3.patch](https://github.com/sagemath/sage-prod/files/10653799/trac_11888_v3.patch.gz) to `$SAGE_ROOT/devel/sage`
fredrik-johansson commented 12 years ago
comment:23

principle -> principal, branchs -> branches

Otherwise, from looking at the patch, seems good.

benjaminfjones commented 12 years ago
comment:24

Thanks for looking at the patch, Fredrik. I've fixed the mistakes and replaced the latest patch.

benjaminfjones commented 12 years ago

adds lambert_w and lambert_w_branch functions

kcrisman commented 12 years ago
comment:25

Attachment: trac_11888_v3.patch.gz

  SciPy is used to evalute
benjaminfjones commented 12 years ago
comment:26
benjaminfjones commented 12 years ago

addressed reviewer issues, changed order of arguments to be consistant with Mma/Maple

benjaminfjones commented 12 years ago

Attachment: trac_11888_v4.patch.gz

fixes random tests after rebasing against #9130

benjaminfjones commented 12 years ago
comment:27

Attachment: trac_11888-random-tests.patch.gz

benjaminfjones commented 12 years ago

Description changed:

--- 
+++ 
@@ -40,4 +40,6 @@

 Apply:

-1. [attachment: trac_11888_v3.patch](https://github.com/sagemath/sage-prod/files/10653799/trac_11888_v3.patch.gz) to `$SAGE_ROOT/devel/sage`
+* Patches at #9130 
+* [attachment: trac_11888_v4.patch](https://github.com/sagemath/sage-prod/files/10653800/trac_11888_v4.patch.gz) to `$SAGE_ROOT/devel/sage`
+* [attachment: trac_11888-random-tests.patch](https://github.com/sagemath/sage-prod/files/10653801/trac_11888-random-tests.patch.gz) to `$SAGE_ROOT/devel/sage`
benjaminfjones commented 12 years ago

Dependencies: #9130

burcin commented 12 years ago

Changed reviewer from Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson to Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson, Burcin Erocal

burcin commented 12 years ago
comment:28

Do we really want to call this function lambert_w_branch()? Can we name it lambert_w()? I would even suggest to add custom printing methods (_print_() and _print_latex_()) to avoid printing the branch argument if it is 0.

If the function is named lambert_w, you can remove the wrapper function lambert_w() and the manual manipulation of the symbol table. In this case, a custom __call__() method would take the place of the wrapper method.

BTW, we should either open a new ticket to add known exact evaluations to _eval_() or do this here:

burcin commented 12 years ago
comment:29

I have one more comment. Sorry for multiple emails.

You should check if the parent is RDF or CDF using is, not ==. In this context, parent is an argument to the _evalf_() method, which overrides the parent() function imported from sage.structure.coerce. I suggest naming the argument parent_d instead of parent. Then you can do:

R = parent_d or parent(z)
if R is float or R is complex or R is RDF or R is CDF: 
    import scipy.special 
    return scipy.special.lambertw(z, n) 
else: 
    import mpmath 
    return mpmath_utils.call(mpmath.lambertw, z, n, parent=parent)
kcrisman commented 12 years ago
comment:30

Replying to @burcin:

Do we really want to call this function lambert_w_branch()? Can we name it lambert_w()? I would even suggest to add custom printing methods (_print_() and _print_latex_()) to avoid printing the branch argument if it is 0.

That's a great idea.

benjaminfjones commented 12 years ago

Changed dependencies from #9130 to #12507

benjaminfjones commented 12 years ago
comment:31

I've written a new patch that includes significant changes compared to the last one. I think I've included all of burcin's suggestions and I think it's much improved now. All tests pass with the patch applied on 5.0.beta4 + #12507.

One thing I haven't managed to figure out is how to get integration to work, e.g.

sage: integrate(lambert_w(x), x)
...
RuntimeError: ECL says: Error executing code in Maxima: lambert_w: wrong number of arguments.

I guess that's because there isn't a two-argument version of lambert_w defined in maxima. The conversion maxima -> Sage works (as shown in one of the doctests) but it looks like the other way doesn't. Another example:

sage: maxima(lambert_w(5))
Maxima ERROR:

lambert_w: wrong number of arguments.
 -- an error. To debug this try: debugmode(true);

Q: How do I get around this?

Numerical integration also fails unless I pass a lambda function:

sage: numerical_integral(lambert_w(x), 0, 1)
Exception TypeError: "function not supported for these types, and can't coerce safely to supported types" in 'sage.gsl.integration.c_ff' ignored
...
(0.0, 0.0)

but ....

sage: numerical_integral(lambda x: lambert_w(x), 0, 1)
(0.33036612476168054, 3.667800782666048e-15)

Q: How do I fix this?

benjaminfjones commented 12 years ago

adds lambert_w function

benjaminfjones commented 12 years ago

Description changed:

--- 
+++ 
@@ -40,6 +40,5 @@

 Apply:

-* Patches at #9130 
-* [attachment: trac_11888_v4.patch](https://github.com/sagemath/sage-prod/files/10653800/trac_11888_v4.patch.gz) to `$SAGE_ROOT/devel/sage`
-* [attachment: trac_11888-random-tests.patch](https://github.com/sagemath/sage-prod/files/10653801/trac_11888-random-tests.patch.gz) to `$SAGE_ROOT/devel/sage`
+* Patch at #12507
+* [attachment: trac_11888_v5.patch](https://github.com/sagemath/sage-prod/files/10653802/trac_11888_v5.patch.gz) to `$SAGE_ROOT/devel/sage`
benjaminfjones commented 12 years ago
comment:32

Attachment: trac_11888_v5.patch.gz

burcin commented 12 years ago
comment:33

Replying to @benjaminfjones:

I've written a new patch that includes significant changes compared to the last one. I think I've included all of burcin's suggestions and I think it's much improved now. All tests pass with the patch applied on 5.0.beta4 + #12507.

Thanks! The patch looks really good. When checking if the input is 0 in _eval_, you might want to return z instead of Integer(0) to preserve the type of the input. Similarly, we should return parent(z)(1) or parent(z)(-1) in the other branches.

> I guess that's because there isn't a two-argument version of lambert_w defined in maxima. The conversion maxima -> Sage works (as shown in one of the doctests) but it looks like the other way doesn't. Another example: > > ``` > sage: maxima(lambert_w(5)) > Maxima ERROR: > > lambert_w: wrong number of arguments. > -- an error. To debug this try: debugmode(true); > ``` > > Q: How do I get around this? You need to define `_maxima_init_evaled_()`. See line 895 of `sage/fuctions/other.py`: http://hg.sagemath.org/sage-main/file/c239be1054e0/sage/functions/other.py#l895
benjaminfjones commented 12 years ago
comment:34

Replying to @burcin:

You need to define _maxima_init_evaled_(). See line 895 of sage/fuctions/other.py:

http://hg.sagemath.org/sage-main/file/c239be1054e0/sage/functions/other.py#l895

It seems that adding _maxima_init_evaled_() solves one issue, converting to Maxima with _maxima_(),

sage: lambert_w(x)._maxima_()
lambert_w(x)
sage: lambert_w(1,x)._maxima_()
...
NotImplementedError: Non-principal branch lambert_w[1](x) is not implemented in Maxima

but integration still doesn't work (same error is raised as before). Looking closer it seems that the issue is here:

sage: z = lambert_w(x)
sage: z.operands()
[0, x]
sage: z.operator()
lambert_w

because when sr_to_max is called in the integration code, I get:

sage: from sage.interfaces.maxima_lib import sr_to_max
sage: sr_to_max(lambert_w(x))
<ECL: ((%LAMBERT_W) 0 $X)>
sage: sr_to_max(lambert_w(1, x))
<ECL: ((%LAMBERT_W) 1 $X)>

and Maxima barfs because it doesn't know what to do with ((%LAMBERT_W) 0 $X).