Opened 3 years ago

Closed 3 years ago

#20939 closed defect (fixed)

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

Reported by: rws Owned by:
Priority: major Milestone: sage-7.5
Component: symbolics Keywords:
Cc: Merged in:
Authors: Ralf Stephan Reviewers: Travis Scrimshaw
Report Upstream: Fixed upstream, in a later stable release. Work issues:
Branch: bb4630f (Commits) Commit: bb4630fa7a2ac6b9b0a1d98486378877fe60f97c
Dependencies: #21614, #21963 Stopgaps:

Description (last modified by 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.

Change History (35)

comment:1 Changed 3 years ago by rws

  • Dependencies set to 16813

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

comment:2 Changed 3 years ago by rws

  • Dependencies changed from 16813 to #16813

comment:3 Changed 3 years ago by 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.

comment:4 Changed 3 years ago by rws

  • Dependencies #16813 deleted
  • Description modified (diff)
  • Summary changed from Remove pexpect-Maxima usage in Y(m,n) to Fix crash and remove pexpect-Maxima usage in Y(m,n)
  • Type changed from enhancement to defect

comment:5 Changed 3 years ago by rws

  • Branch set to u/rws/remove_pexpect_maxima_usage_in_y_m_n_

comment:6 Changed 3 years ago by rws

  • Authors set to Ralf Stephan
  • Commit set to f233afa35058ab40a1f6f33e6d635e6eac515f97
  • Milestone changed from sage-7.3 to sage-7.4
  • Status changed from new to needs_review

New commits:

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

comment:7 Changed 3 years ago by tscrim

  • Reviewers set to Travis Scrimshaw

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)

comment:8 Changed 3 years ago by git

  • Commit changed from f233afa35058ab40a1f6f33e6d635e6eac515f97 to 1a5cf6b9dbf75041e634a5c0d27b0a2a09932834

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

1a5cf6b20939: cosmetics

comment:9 Changed 3 years ago by rws

  • Status changed from needs_review to positive_review

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.

comment:10 Changed 3 years ago by tscrim

  • Branch changed from u/rws/remove_pexpect_maxima_usage_in_y_m_n_ to u/tscrim/remove_pexpect_maxima_in_ymn-20939
  • Commit changed from 1a5cf6b9dbf75041e634a5c0d27b0a2a09932834 to a080e1e8824bf56a760d54cbd7e6c9b79c351cef
  • Status changed from positive_review to needs_review

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


New commits:

a080e1eRemoving uncessary ... in doctest.

comment:11 Changed 3 years ago by rws

  • Status changed from needs_review to positive_review

comment:12 Changed 3 years ago by vbraun

  • Status changed from positive_review to needs_work

Documentation doesn't build, see patchbot

comment:13 Changed 3 years ago by 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.

comment:14 Changed 3 years ago by 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

comment:15 Changed 3 years ago by 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))

comment:16 Changed 3 years ago by 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.

Last edited 3 years ago by rws (previous) (diff)

comment:17 Changed 3 years ago by rws

  • Dependencies set to #21614

comment:18 Changed 3 years ago by rws

  • Status changed from needs_work to positive_review

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

comment:19 Changed 3 years ago by vbraun

  • Status changed from positive_review to needs_work
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

comment:20 Changed 3 years ago by fbissey

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

comment:21 Changed 3 years ago by rws

  • Branch changed from u/tscrim/remove_pexpect_maxima_in_ymn-20939 to u/rws/remove_pexpect_maxima_in_ymn-20939

comment:22 Changed 3 years ago by rws

  • Commit changed from a080e1e8824bf56a760d54cbd7e6c9b79c351cef to 31e6c649c586363c621f2daabf805283be10775a
  • Status changed from needs_work to needs_review

New commits:

55deed7Merge branch 'develop' into t/20939/remove_pexpect_maxima_in_ymn-20939
31e6c6420939: fix NaN handling

comment:23 Changed 3 years ago by tscrim

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

comment:24 Changed 3 years ago by git

  • Commit changed from 31e6c649c586363c621f2daabf805283be10775a to 04f79cdf8ed0c887b7f7d1797ce5023f782d0a50

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

04f79cd20939: make sure atan2(0,0) returns always NaN

comment:25 Changed 3 years ago by 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

comment:26 Changed 3 years ago by tscrim

  • Status changed from needs_review to positive_review

Thanks. Now back to the buildbots.

comment:27 Changed 3 years ago by 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]

comment:28 Changed 3 years ago by rws

  • Status changed from positive_review to needs_work

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.

comment:29 Changed 3 years ago by vbraun

Same here...

comment:30 Changed 3 years ago by 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 :(

comment:31 Changed 3 years ago by rws

  • Dependencies changed from #21614 to #21614, pynac-0.7.1
  • Report Upstream changed from N/A to Fixed upstream, in a later stable release.

comment:32 Changed 3 years ago by rws

  • Dependencies changed from #21614, pynac-0.7.1 to #21614, #21963
  • Status changed from needs_work to needs_review

comment:33 Changed 3 years ago by git

  • Commit changed from 04f79cdf8ed0c887b7f7d1797ce5023f782d0a50 to bb4630fa7a2ac6b9b0a1d98486378877fe60f97c

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

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

comment:34 Changed 3 years ago by tscrim

  • Milestone changed from sage-7.4 to sage-7.5
  • Status changed from needs_review to positive_review

Let's try again.

comment:35 Changed 3 years ago by vbraun

  • Branch changed from u/rws/remove_pexpect_maxima_in_ymn-20939 to bb4630fa7a2ac6b9b0a1d98486378877fe60f97c
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.