BUG: Fix ellipj when m is near 1.0 (#8480) #8548

FormerPhysicist · 2018-03-12T03:30:23Z

This pull request fixes ellipj when m is near 1.0. Previously, ellipj was using equations 16.15.1-16.15.4 from Abramowitz and Stegun's Handbook of Mathematical Functions to calculate sn, cn, dn, and phi when m was near 1.0.

However, these formulas were only valid when phi < pi/2 (or equivalently when u<K, where K=ellipk(m) is the quarter-period) and so ellipj failed when u>K.

To fix this issue, this pull request will instead use formula 16.15.4 to calculate the "principal" part of the amplitude phi, or the remainder when dphi=phi/(pi/2). The total phi will be set to dphi+n*pi/2 where n=floor(u/K) (basically the number of quarter periods in u) and this total phi will be used to correctly compute sn=sin(phi) and cn=cos(phi). With this fix, ellipj will produce the right values for all u, when m is very near 1.0.

Closes #8480

larsoner

Changes look reasonable to me, I've found Abramowitz and Stegun useful in the past, and the tests appear to validate the functionality. But it would be better if someone with more knowledge in this domain could look!

larsoner · 2018-03-12T14:02:55Z

scipy/special/add_newdocs.py

+    except when `m` is within 1e-9 of 0 or 1.  
+
+    In the latter case, when m is near 1, the functions `sn`, `cn`, and `dn`
+    are found using an approximation for `phi` given by 16.15.4 in [2]. 


[2]_ to make it link properly.

See CircleCI output for example:

https://5997-1460385-gh.circle-artifacts.com/0/html-scipyorg/generated/scipy.special.ellipj.html#scipy.special.ellipj

@larsoner fixed

larsoner · 2018-03-12T14:03:22Z

scipy/special/add_newdocs.py

+    are found using an approximation for `phi` given by 16.15.4 in [2]. 
+    As this approximation is only valid when `phi < pi/2` (i.e.,
+    when `u<K`), the computation only uses 16.15.4 to find the "principal" 
+    part of `phi` or the remainder `dphi=phi%(pi/2)`. The returned value 


double-backticks on code for it to render properly

Which parts? Do you just mean the dphi=phi%(pi/2) part? Or everywhere I've used one set of backticks (e.g., should phi just be italicized (one pair of ) or as code (two pairs of )?)

larsoner · 2018-03-12T14:08:28Z

scipy/special/cephes/ellpj.c

+        /* Evaluate the complete elliptic integral using the
+           1-m<<1 approximation in ellpk.c */
+        K = ellpk(1.0-m);
+	uabs=fabs(u);


indentation off by one

larsoner · 2018-03-12T14:09:51Z

scipy/special/tests/test_mpmath.py


-    FuncData(dn, dataset, (0, 1), 2, rtol=1e-10).check()
+    # Add the points that will use the m near 1 method
+    m = 0.99999999999


it's probably worth going a step further. In addition to the one right at the threshold (0.9999999999), how about one just outside it (0.9999999998) and beyond it (0.99999999999)? It helps validate the choice of threshold.

Is it better to write 1. - 1e-10 for readability of the digits?

Sounds good to me.

ilayn · 2018-03-13T13:45:08Z

Other than the PEP8 whitespace after comma issues in test_mpmath.py also looks good to me. But indeed a domain expert look would be better for the contents.

person142 · 2018-03-24T15:57:04Z

Before I review this in earnest:

The tab -> spaces changes in ellpj.c are reasonable, but better done in a separate PR. That can get merged quickly and then the code changes here will be clearer.
Please remove the unrelated whitespace changes in test_mpmath.py

FormerPhysicist · 2018-03-26T03:34:32Z

@person142 Done

person142 · 2018-03-27T14:03:17Z

Can you comment on your general approach versus e.g. just doing what Boost does (reading their source code is fine):

http://www.boost.org/doc/libs/1_52_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/jacobi/jacobi_elliptic.html

In general argument reduction by irrational numbers (here pi/2) is inaccurate and should be avoided if possible.

ev-br · 2018-03-31T23:52:12Z

scipy/special/add_newdocs.py

+    In the latter case, when m is near 1, the functions `sn`, `cn`, and `dn`
+    are found using an approximation for `phi` given by 16.15.4 in [2]_. 
+    As this approximation is only valid when `phi < pi/2` (i.e.,
+    when `u<K`), the computation only uses 16.15.4 to find the "principal" 


ReST weirdness: single backticks are for references, double backticks for code

Can you elaborate what you mean by 'reference'? Looking at the other function docstrings, it appears that they use a single backtick for variables and math equations (sometimes with the ReST role(?) keyword :math:). Should I change the existing text for ellpj to use double backticks?

…8480

FormerPhysicist · 2018-07-26T23:27:28Z

I've finally gotten around to looking at the source code for boost version 1.67 implementation of ellpj and there appear to be two main differences with scipy's current implementation:

They handle the m>1 case using equations 16.11 from Abramowitz. Currently, scipy doesn't attempt to calculate ellpj at all when m>1 and just returns NAN.
They don't provide any special handling when m is close to unity. They just continue to use the arithmetic-geometric mean algorithm* (AGM). Scipy on the other hand tries to use the approximations in equations 16.15 of Abramowitz which are valid for m close to 1 and small u. My pull request would fix our use of these approximations so that u can be large.

However, what's also interesting, is that boost had a section that attempted to treat the m close to unity case, but it's been commented out with the comments:

/* Can't get this to work to adequate precision - disabled for now...
//
// Asymptotic forms from A&S 16.15:

Skimming their code, it looks like they were trying the same naive implementation of 16.15 that scipy does and I wonder if they were getting the same error that my pull request addresses.

So the TL;DR: either we stick with AGM for m close to unity, or we can use the modified version of 16.15 proposed by vasilevgmaslov

(*Note that AGM seems to be a standard way of getting ellpj for m<1; scipy uses it when m is < 1 and not close to unity.)

tylerjereddy · 2020-05-19T23:42:43Z

@ilayn @person142 @larsoner Any realistic chance this is ready and merged in ~1 week for branching? No reviewer activity in almost two years.

person142 · 2020-05-20T01:35:16Z

@tylerjereddy I removed the milestone.

FormerPhysicist added 3 commits March 11, 2018 20:05

BUG: fixing Jacobian elliptic functions when m is near 1

42684da

DOC: updated docstring for ellipj

7410c27

TST: adding new tests for ellipj

6b8bbf1

larsoner reviewed Mar 12, 2018

View reviewed changes

FormerPhysicist added 3 commits March 12, 2018 09:29

DOC: fixed link in docstring for ellipj

8ac55a2

TST: adding points outside and inside m-near-1 boundary

c709127

BUG: small spacing fixes to ellpj.c

1000fc1

TST: fix E231 missing whitespace after ',' error

ff2d2b8

rgommers added scipy.special maintenance Items related to regular maintenance tasks labels Mar 17, 2018

FormerPhysicist added 2 commits March 25, 2018 19:53

TST: reverting unrelated whitespace fixes

ec3db66

BUG: reverting spacing issues in ellpj.c

f9440e5

ev-br reviewed Mar 31, 2018

View reviewed changes

rgommers added the needs-work Items that are pending response from the author label Jun 10, 2018

Merge branch 'master' of github.com:FormerPhysicist/scipy into issue-…

7e0abf8

…8480

Merging master.

e1df379

tylerjereddy added this to the 1.5.0 milestone Jan 13, 2020

FormerPhysicist added 2 commits January 13, 2020 18:03

Merge branch 'master' of github.com:scipy/scipy into issue-8480

7bbbff0

Merge branch 'master' of github.com:scipy/scipy into issue-8480

7ad5f97

person142 removed this from the 1.5.0 milestone May 20, 2020

FormerPhysicist added 3 commits June 27, 2020 11:51

Merge branch 'master' of github.com:scipy/scipy into issue-8480

56ef8f7

Merge branch 'master' of github.com:scipy/scipy into issue-8480

c18bb02

MAINT: fix ellipj algorithm for m near 1.

f23cbca

FormerPhysicist mentioned this pull request Jan 1, 2021

BUG: adding Buhring's analytic continuation method to special.hyp2f1 (redux) #13310

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: Fix ellipj when m is near 1.0 (#8480) #8548

BUG: Fix ellipj when m is near 1.0 (#8480) #8548

FormerPhysicist commented Mar 12, 2018

larsoner left a comment

larsoner Mar 12, 2018

larsoner Mar 12, 2018

FormerPhysicist Mar 12, 2018

larsoner Mar 12, 2018

FormerPhysicist Mar 12, 2018

larsoner Mar 12, 2018

FormerPhysicist Mar 12, 2018

larsoner Mar 12, 2018

ilayn Mar 12, 2018

larsoner Mar 12, 2018

FormerPhysicist Mar 12, 2018

ilayn commented Mar 13, 2018

person142 commented Mar 24, 2018

FormerPhysicist commented Mar 26, 2018

person142 commented Mar 27, 2018

ev-br Mar 31, 2018

FormerPhysicist Jul 26, 2018

FormerPhysicist commented Jul 26, 2018

tylerjereddy commented May 19, 2020

person142 commented May 20, 2020

BUG: Fix ellipj when m is near 1.0 (#8480) #8548

Are you sure you want to change the base?

BUG: Fix ellipj when m is near 1.0 (#8480) #8548

Conversation

FormerPhysicist commented Mar 12, 2018

larsoner left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ilayn commented Mar 13, 2018

person142 commented Mar 24, 2018

FormerPhysicist commented Mar 26, 2018

person142 commented Mar 27, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FormerPhysicist commented Jul 26, 2018

tylerjereddy commented May 19, 2020

person142 commented May 20, 2020