Sage is missing the lambert_w function

[benjaminfjones]

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:

• attachment: trac_11888_v8.patch to $SAGE_ROOT/devel/sage

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]

Attachment: trac_11888.patch.gz

add lambert_w symbolic function

[benjaminfjones]

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

[benjaminfjones]

Author: Benjamin Jones

[benjaminfjones]

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]

apply to $SAGE_ROOT/devel/sage

[kini]

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]

comment:4

All tests pass.

[benjaminfjones]

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]

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]

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]

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]

comment:10

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

[kcrisman]

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]

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]

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]

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]

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]

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]

Work Issues: add second parameter, RDF/CDF stuff

[kcrisman]

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

[benjaminfjones]

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]

comment:51

I can make a quick spelling fix patch, but I'll wait and see if anyone else has changes to suggest.

[sagetrac-dsm]

comment:52

After some real-world discussions, I'd prefer to avoid falling into Python ints for (some of) the special values:

sage: parent(lambert_w(0))
Integer Ring
sage: parent(lambert_w(e))
<type 'int'>
sage: parent(lambert_w(-1/e))
<type 'int'>
sage: parent(lambert_w(SR(-1/e)))
<type 'int'>
sage: parent(lambert_w(SR(0)))
Integer Ring

Mysteriously enough, instrumenting it reveals that _eval_ is actually returning SR(1) which then in some part of the code I don't understand becomes int(1) before we get it back. If we explicitly return Integer(1) then it seems to stay as Integer(1). This isn't the biggest deal in the world but there have been several bug reports caused by something falling out of Sagespace into Pythonspace.

[sagetrac-dsm]

comment:53

It looks like whatever happens after _eval_ might "dereference" the SR; the call seems to give Integer(1) if _eval_ returns SR(Integer(1)). Probably someone who actually knows what's going on could explain it in one line.

[benjaminfjones]

comment:54

One solution is to change the return statements in these special cases where the automatic simplification returns an integer to just explicitly return Integer(1), etc..

How does that sound?

[benjaminfjones]

comment:55

New patch is ready for review. I hope this can be the final revision!

Patchbot: apply attachment: trac_11888_v8.patch to $SAGE_ROOT/devel/sage

[sagetrac-dsm]

comment:56

Looking at it now.

[benjaminfjones]

fixed spelling / grammer mistakes, returned parent(Integer(...)) for special values

[sagetrac-dsm]

comment:57

Attachment: trac_11888_v8.patch.gz

Okay, this looks good. Two copyedits, an extra doc describing the behaviours of the derivative function, some tests making sure we can't differentiate with respect to the branch number, and the addition of lambert_W(-pi/2) = pi/2*I as a special value.

I give positive review to the preexisting parts of v8; if the new bits of v9 look okay I think we're good to go.

[sagetrac-dsm]

Attachment: trac_11888_v9.patch.gz

post-review version of lambertw sf

[benjaminfjones]

comment:58

New changes look good; good catch about the derivative w.r.t. branch. All relavent tests pass for me on sage-5.0. I would say this is ready to go in! Thanks for the very thorough review, Doug.

[sagetrac-dsm]

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

[sagetrac-dsm]

comment:60

+1. Tx for the work!

[williamstein]

Changed reviewer from Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson, Burcin Erocal, Douglas McNeil to Keshav Kini, Karl-Dieter Crisman, Fredrik Johansson, Burcin Erocal, Douglas McNeil, William Stein

[williamstein]

comment:61

This patch is awesome! It's also a great example of how to make a well-documented new symbolic function that illustrates many issues. Here are a few trivial nitpicks:

• What is "simplication"?

When automatic simplication occurs, the parent of the output value should be 

• This docstring should start with r" since it contains a backslash:

        646	        """ 
 	647	        The derivative of `W_n(x)` is `W_n(x)/(x \cdot W_n(x) + x)`. 

(check for similar instances throughout).

• Don't use periods at the end of exceptions (also don't capitalize). Many instances of this being wrong, e.g.,

 	679	            raise ValueError("Derivative not defined with respect to the branch number.") 

Here's a good example of what an exception string should look like (built into python):

>>> 1/0
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: integer division or modulo by zero

[sagetrac-dsm]

comment:62

Updated version taking into account comments of was.

[sagetrac-dsm]

minor edits

[williamstein]

comment:64

Attachment: trac_11888_v10.patch.gz

[jdemeyer]

Merged: sage-5.1.beta4