Fix crash and remove pexpect-Maxima usage in Y(m,n)

[rws]

sage: spherical_harmonic(3, 2, 1, 2*pi/3)
...
TypeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.

The eval function of SphericalHarmonic uses pexpect-Maxima. Since this needs only a formula the call to Maxima could be replaced by a few lines of Python.

​http://dlmf.nist.gov/14.30.E1

Replacing pexpect-Maxima is part of #17753.

[rws]

This depends on the associated Legendre polynomials being implemented, e,g, #16813.

[rws]

It won't be dependent on P(l,m,x) if the P special value formula is directly used in _eval_. The Y floating point evaluation is done via mpmath anyway.

[rws]

New commits:

​f233afa

20939: Fix crash and remove pexpect-Maxima usage in Y(m,n)

[tscrim]

LGTM modulo

         Check that :trac:`20939` is fixed::
 
             sage: spherical_harmonic(3,2,1,2*pi/3)
             1/240*sqrt(30)*(-15*I*sqrt(7)*sqrt(3) -
-            ...15*sqrt(7))*cos(1)*sin(1)^2/sqrt(pi)
+            ... + 15*sqrt(7))*cos(1)*sin(1)^2/sqrt(pi)
         """
         if n in ZZ and m in ZZ and n > -1:
             if abs(m) > n:
                 return ZZ(0)
-            if (m == 0 and theta.is_zero()):
+            if m == 0 and theta.is_zero():
                 return sqrt((2*n+1)/4/pi)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1a5cf6b

20939: cosmetics

[rws]

That doesn't work as given. First the plus is a minus, and then the space after the ellipsis will make it fail too.

I'll assume the last commit is still what you wanted and set positive. Just revert if not. Thanks for the review.

[tscrim]

So the ... is actually unneeded; I was thinking there was more than 2 terms. Please check my changes.

---

New commits:

​a080e1e

Removing uncessary ... in doctest.

[vbraun]

Documentation doesn't build, see patchbot

[rws]

Interesting. I get this

[dochtml] [plot3d   ] /home/ralf/sage/src/doc/en/reference/plot3d/sage/plot/plot3d/plot3d.rst:1937: WARNING: Exception occurred in plotting plot3d-35
[dochtml] [plot3d   ] from /home/ralf/sage/src/doc/en/reference/plot3d/sage/plot/plot3d/plot3d.rst:
[dochtml] [plot3d   ] Traceback (most recent call last):
[dochtml] [plot3d   ] File "/home/ralf/sage/local/lib/python2.7/site-packages/matplotlib-1.5.1-py2.7-linux-x86_64.egg/matplotlib/sphinxext/plot_directive.py", line 517, in run_code

I don't know how to get to the plot3d-35 plot commands but:

[dochtml] [plot3d   ] result[0] = interp_rdf(args
[dochtml] [plot3d   ] File "sage/symbolic/function.pyx", line 821, in sage.symbolic.function.GinacFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:9409)
[dochtml] [plot3d   ] res = super(GinacFunction, self).__call__(*args, **kwds)
[dochtml] [plot3d   ] File "sage/symbolic/function.pyx", line 972, in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:11034)
[dochtml] [plot3d   ] res = super(BuiltinFunction, self).__call__(
[dochtml] [plot3d   ] File "sage/symbolic/function.pyx", line 488, in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6886)
[dochtml] [plot3d   ] res = g_function_eval2(self._serial, (<Expression>args[0])._gobj,
[dochtml] [plot3d   ] RuntimeError: arctan2_eval(): arctan2(0,0) encountered

The only calls to atan2(x,y) are in revolution_plot3d() and if I remove all its doctests it still crashes on build at this point, so I no longer know where to look. Either the 3d plot number 35 or the uncompiled plot3d.rst would be necessary.

Apart from the bug at hand I think Pynac's atan2_eval should return unsigned infinity instead of throwing an exception. There is a reason atan2(0,0) is called however.

[rws]

This is the trigger:

sage: phi, theta = var('phi, theta')
sage: Y = spherical_harmonic(2, 1, theta, phi)
sage: rea = spherical_plot3d(abs(real(Y)), (phi,0,2*pi), (theta,0,pi), color='bl
....: ue', opacity=0.6)
sage: ima = spherical_plot3d(abs(imag(Y)), (phi,0,2*pi), (theta,0,pi), color='re
....: d', opacity=0.6)
sage:  (rea + ima).show(aspect_ratio=1)
---------------------------------------------------------------------------
RuntimeError: arctan2_eval(): arctan2(0,0) encountered

[rws]

The new code makes Y expand immediately, and the real/imag calls then really blow it up.

sage: Y = 1/4*sqrt(6)*sqrt(5)*sqrt(sin(theta)^2)*cos(theta)*e^(I*phi)/sqrt(pi)
sage: abs(real(Y))
abs(1/4*sqrt(6)*sqrt(5)*sqrt(abs(sin(theta))^2)*cos(1/2*\
arctan2(2*cos(real_part(theta))*cosh(imag_part(theta))*sin(real_part(theta))*sinh(imag_part(theta)),
cosh(imag_part(theta))^2*sin(real_part(theta))^2-cos(real_part(theta))^2*sinh(imag_part(theta))^2))
*cos(real_part(phi))*cos(real_part(theta))*cosh(imag_part(theta))*e^(-imag_part(phi))/sqrt(pi) - 
1/4*sqrt(6)*sqrt(5)*sqrt(abs(sin(theta))^2)*cos(real_part(theta))*cosh(imag_part(theta))*e^(-imag_part(phi))*sin(1/2*arctan2(2*cos(real_part(theta))*cosh(imag_part(theta))*sin(real_part(theta))*sinh(imag_part(theta)), cosh(imag_part(theta))^2*sin(real_part(theta))^2 - cos(real_part(theta))^2*sinh(imag_part(theta))^2))*sin(real_part(phi))/sqrt(pi) + 1/4*sqrt(6)*sqrt(5)*sqrt(abs(sin(theta))^2)*cos(real_part(phi))*e^(-imag_part(phi))*sin(1/2*arctan2(2*cos(real_part(theta))*cosh(imag_part(theta))*sin(real_part(theta))*sinh(imag_part(theta)), cosh(imag_part(theta))^2*sin(real_part(theta))^2 - cos(real_part(theta))^2*sinh(imag_part(theta))^2))*sin(real_part(theta))*sinh(imag_part(theta))/sqrt(pi) + 1/4*sqrt(6)*sqrt(5)*sqrt(abs(sin(theta))^2)*cos(1/2*arctan2(2*cos(real_part(theta))*cosh(imag_part(theta))*sin(real_part(theta))*sinh(imag_part(theta)), cosh(imag_part(theta))^2*sin(real_part(theta))^2 - cos(real_part(theta))^2*sinh(imag_part(theta))^2))*e^(-imag_part(phi))*sin(real_part(phi))*sin(real_part(theta))*sinh(imag_part(theta))/sqrt(pi))

[rws]

sage: x=abs(real(Y))
sage: x.subs(theta==0)
...
RuntimeError: arctan2_eval(): arctan2(0,0) encountered
sage: y=abs(imag(Y))
sage: y.subs(theta==0)
...
RuntimeError: arctan2_eval(): arctan2(0,0) encountered
sage: Y
1/4*sqrt(6)*sqrt(5)*sqrt(sin(theta)^2)*cos(theta)*e^(I*phi)/sqrt(pi)
sage: Y.subs(theta==0)
0

So it seems the real/imag expansions produce results in some cases that are not defined at zero, so are not correct.

[rws]

So to summarize the code here is good, we just depend on #21614.

[vbraun]

File "src/sage/plot/plot3d/plot3d.py", line 1195, in sage.plot.plot3d.plot3d.spherical_plot3d
Failed example:
    (rea + ima).show(aspect_ratio=1)  # long time (4s on sage.math, 2011)
Exception raised:
    Traceback (most recent call last):
      File "/home/chapoton/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 498, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/chapoton/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 861, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.plot.plot3d.plot3d.spherical_plot3d[13]>", line 1, in <module>
        (rea + ima).show(aspect_ratio=Integer(1))  # long time (4s on sage.math, 2011)
      File "sage/plot/plot3d/base.pyx", line 1379, in sage.plot.plot3d.base.Graphics3d.show (build/cythonized/sage/plot/plot3d/base.c:17925)
        dm.display_immediately(self, **kwds)
      File "/home/chapoton/sage/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.py", line 791, in display_immediately
        plain_text, rich_output = self._rich_output_formatter(obj, rich_repr_kwds)
      File "/home/chapoton/sage/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.py", line 623, in _rich_output_formatter
        rich_output = self._call_rich_repr(obj, rich_repr_kwds)
      File "/home/chapoton/sage/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.py", line 581, in _call_rich_repr
        return obj._rich_repr_(self, **rich_repr_kwds)
      File "sage/plot/plot3d/base.pyx", line 130, in sage.plot.plot3d.base.Graphics3d._rich_repr_ (build/cythonized/sage/plot/plot3d/base.c:4422)
        opts = self._process_viewing_options(kwds)
      File "sage/plot/plot3d/base.pyx", line 1257, in sage.plot.plot3d.base.Graphics3d._process_viewing_options (build/cythonized/sage/plot/plot3d/base.c:17397)
        self._determine_frame_aspect_ratio(opts['aspect_ratio'])
      File "sage/plot/plot3d/base.pyx", line 472, in sage.plot.plot3d.base.Graphics3d._determine_frame_aspect_ratio (build/cythonized/sage/plot/plot3d/base.c:9247)
        a_min, a_max = self._safe_bounding_box()
      File "sage/plot/plot3d/base.pyx", line 490, in sage.plot.plot3d.base.Graphics3d._safe_bounding_box (build/cythonized/sage/plot/plot3d/base.c:9420)
        a_min, a_max = self.bounding_box()
      File "sage/plot/plot3d/base.pyx", line 1831, in sage.plot.plot3d.base.Graphics3dGroup.bounding_box (build/cythonized/sage/plot/plot3d/base.c:24445)
        v = [obj.bounding_box() for obj in self.all]
      File "sage/plot/plot3d/parametric_surface.pyx", line 379, in sage.plot.plot3d.parametric_surface.ParametricSurface.bounding_box (build/cythonized/sage/plot/plot3d/parametric_surface.c:5079)
        self.triangulate()
      File "sage/plot/plot3d/parametric_surface.pyx", line 427, in sage.plot.plot3d.parametric_surface.ParametricSurface.triangulate (build/cythonized/sage/plot/plot3d/parametric_surface.c:5826)
        raise
      File "sage/plot/plot3d/parametric_surface.pyx", line 422, in sage.plot.plot3d.parametric_surface.ParametricSurface.triangulate (build/cythonized/sage/plot/plot3d/parametric_surface.c:5746)
        self.eval_grid(urange, vrange)
      File "sage/plot/plot3d/parametric_surface.pyx", line 601, in sage.plot.plot3d.parametric_surface.ParametricSurface.eval_grid (build/cythonized/sage/plot/plot3d/parametric_surface.c:7437)
        (<Wrapper_rdf>fx).call_c(uv, &self.vs[i*n+j].x)
      File "sage/ext/interpreters/wrapper_rdf.pyx", line 90, in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.call_c (build/cythonized/sage/ext/interpreters/wrapper_rdf.c:2163)
        result[0] = interp_rdf(args
      File "sage/symbolic/expression.pyx", line 1410, in sage.symbolic.expression.Expression.__float__ (build/cythonized/sage/symbolic/expression.cpp:11032)
        ret = float(self._eval_self(float))
      File "sage/symbolic/expression.pyx", line 1197, in sage.symbolic.expression.Expression._eval_self (build/cythonized/sage/symbolic/expression.cpp:9741)
        res = self._gobj.evalf(0, {'parent':R})
      File "sage/libs/pynac/pynac.pyx", line 2160, in sage.libs.pynac.pynac.py_eval_constant (build/cythonized/sage/libs/pynac/pynac.cpp:25827)
        constant = constants_table[serial]
    KeyError: 3

[fbissey]

I should say I saw it too. plot3d-35 is still problematic with this here.

[rws]

New commits:

​55deed7	Merge branch 'develop' into t/20939/remove_pexpect_maxima_in_ymn-20939
​31e6c64	20939: fix NaN handling

[tscrim]

Would it be possible to have a test to make sure this stays fixed?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​04f79cd

20939: make sure atan2(0,0) returns always NaN

[rws]

This commit contains a test for the specific NaN evaluation problem causing the second crash (the fix for the first crash is tested with #21614). The commit also fixes an inconsistency. As to a general test of the second problem I opened ​https://github.com/pynac/pynac/issues/211

[tscrim]

Thanks. Now back to the buildbots.

[fbissey]

I don know what Volker will report if anything but now I have a segfault in constant.py

sage: sig_on_count() # check sig_on/off pairings (virtual doctest) ## line 688 ##
0
sage: loads(dumps(NaN)) ## line 704 ##
------------------------------------------------------------------------
/usr/lib64/python2.7/site-packages/cysignals/signals.so(+0x4556)[0x7f4f180eb556]
/usr/lib64/python2.7/site-packages/cysignals/signals.so(+0x4a6c)[0x7f4f180eba6c]
/usr/lib64/python2.7/site-packages/cysignals/signals.so(+0x6d6f)[0x7f4f180edd6f]
/lib64/libpthread.so.0(+0x10d70)[0x7f4f1ce53d70]
/usr/lib64/python2.7/site-packages/sage/libs/pynac/pynac.so(+0x2940d)[0x7f4d3ece640d]
/usr/lib64/libpynac.so.8(_ZN5GiNaC8constant9unarchiveERKNS_12archive_nodeERNS_9containerINSt7__cxx114listEEE+0x196)[0x7f4d3e955dcc]
/usr/lib64/libpynac.so.8(_ZNK5GiNaC12archive_node9unarchiveERNS_9containerINSt7__cxx114listEEE+0xf4)[0x7f4d3e936836]
/usr/lib64/libpynac.so.8(_ZNK5GiNaC7archive12unarchive_exERKNS_9containerINSt7__cxx114listEEEj+0x130)[0x7f4d3e936a42]
/usr/lib64/python2.7/site-packages/sage/symbolic/expression.so(+0xec87f)[0x7f4d3e63787f]
/usr/lib64/python2.7/site-packages/sage/symbolic/expression.so(+0xecc2a)[0x7f4d3e637c2a]
/usr/lib64/libpython2.7.so.1.0(PyCFunction_Call+0xf1)[0x7f4f1d0ddbd9]
/usr/lib64/libpython2.7.so.1.0(PyObject_Call+0x5c)[0x7f4f1d0a8123]

[rws]

Confirmed. I'm sorry. With Pynac changes I have a testing procedure in place here, but with Python I was relying on patchbot. This issue is also somewhat unusual as to the surprising problems presented.

[vbraun]

Same here...

[fbissey]

Positive stuff for me, it made me realise that some bits of cysignals are not currently working properly on Gentoo - even so the tests pass :(

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​bb4630f

Merge branch 'develop' into t/20939/remove_pexpect_maxima_in_ymn-20939

[tscrim]

Let's try again.