Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(771)

Issue 20078: Fixed Bessel I implementation for Scipy (Closed)

Can't Edit
Can't Publish+Mail
Start Review
Created:
5 years, 1 month ago by pav
Modified:
4 years, 9 months ago
Visibility:
Public.

Description

Aims to address
http://scipy.org/scipy/scipy/ticket/849
http://scipy.org/scipy/scipy/ticket/503

It's a reimplementation of Iv, based on code converted to C from Boost. Licenses
should be OK.

I haven't checked if the netlib.org/specfun (not the Specfun in Scipy!) codes
would be better than this. Do they have license issues?

Patch Set 1 #

Total comments: 29

Patch Set 2 : Correct C style #

Unified diffs Side-by-side diffs Delta from patch set Stats (+699 lines, -139 lines) Patch
M scipy/special/cephes/const.c View 1 1 chunk +2 lines, -0 lines 0 comments Download
M scipy/special/cephes/hyperg.c View 4 chunks +13 lines, -6 lines 0 comments Download
M scipy/special/cephes/iv.c View 1 chunk +0 lines, -119 lines 0 comments Download
A scipy/special/cephes/scipy_iv.c View 1 1 chunk +648 lines, -0 lines 0 comments Download
M scipy/special/tests/test_basic.py View 4 chunks +36 lines, -14 lines 0 comments Download

Messages

Total messages: 6
njs
Read through, didn't see any issues.
5 years, 1 month ago #1
charlesr.harris
Do you want to leave the fixed up function in cephes? If we are going ...
5 years, 1 month ago #2
charlesr.harris
http://codereview.appspot.com/20078/diff/1/2 File scipy/special/cephes/const.c (right): http://codereview.appspot.com/20078/diff/1/2#newcode64 Line 64: double EULER = 0.57721566490153286061; /* Euler constant */ ...
5 years, 1 month ago #3
charlesr.harris
http://codereview.appspot.com/20078/diff/1/3 File scipy/special/cephes/hyperg.c (left): http://codereview.appspot.com/20078/diff/1/3#oldcode87 Line 87: double hyperg( a, b, x) The names hyperg ...
5 years, 1 month ago #4
stefanv
Hey all, Charles, a question for you. Pauli, wondering about how to check for nan's. ...
5 years, 1 month ago #5
pav
5 years, 1 month ago #6
Charles wrote:
> Do you want to leave the fixed up function in cephes? If 
> we are going to have a verified set of functions it might
> be better to move it to a separate directory.

It uses the error handling, constants, and gamma function from Cephes, so we'd
have to copy those to the core library, too. I changed the file name to indicate
that it's a significantly revised version, though.

I think the 'core' library should be introduced separately later, as there are
build issues etc. to think about.

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (left):

http://codereview.appspot.com/20078/diff/1/3#oldcode87
Line 87: double hyperg( a, b, x)
On 2009/02/26 04:19:35, charlesr.harris wrote:
> The names hyperg and hy1f1p are a bit cryptic. I don't know what we want to do
> about that, though.

These are the Cephes names. I'd rather not go changing them.

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (right):

http://codereview.appspot.com/20078/diff/1/3#newcode207
Line 207: *err = 1.0;
On 2009/02/26 22:15:46, Stefan wrote:
> Looks like the whole file needs some re-indentation.  According to the Python
3k
> standard that we follow, 4 spaces per indent.

Yes, the Cephes indentation is horrible. But I wanted to make any
re-indentations as a separate SVN commit, so that it's more clear what has been
changed and what hasn't.

http://codereview.appspot.com/20078/diff/1/3#newcode288
Line 288: if (acanc != acanc)
> but the recommended way to test for NaN is with the isnan

isnan is a C99 extension, so no other alternatives than testing for (x == x).

http://codereview.appspot.com/20078/diff/1/5
File scipy/special/cephes/scipy_iv.c (right):

http://codereview.appspot.com/20078/diff/1/5#newcode96
Line 96: v = -v;		/* symmetry */
On 2009/02/26 22:15:46, Stefan wrote:
> On 2009/02/26 04:13:12, charlesr.harris wrote:
> > Trailing comments tend to be fragile. I'm not sure this one is needed, but
> > prefer inline comments where possible.
> 
> Hey Charles, out of curiosity, what do you mean by "fragile"?

I'd rather not spend time on changing these to inline comments, unless there's a
good reason to do it.

http://codereview.appspot.com/20078/diff/1/5#newcode107
Line 107: if (v != 2.0 * floor(v / 2.0))
On 2009/02/26 04:13:12, charlesr.harris wrote:
> This is a style question. I like all the statements in if statements to be
> bracketed as it makes them stand out more. Folks differ about this, but we
> should settle on a style. 

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode118
Line 118: } else
On 2009/02/26 04:13:12, charlesr.harris wrote:
> When one branch of an if statement is bracketed, the other should be also,
even
> when it is only a single statement.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode124
Line 124: /* Uniform asymptotic expansion for large orders.
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Multiline comments like this:
> 
> /*
>  * blah
>  * blah
>  */
> 
> I'm more or less following linux here, with modifications for python.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode130
Line 130: } else {
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Another one way or the other thingie. I prefer
> 
> }
> else {
> 
> But in any case we should settle on a style.

I think that's ugly, but changed nevertheless.

http://codereview.appspot.com/20078/diff/1/5#newcode159
Line 159: factor = (mu - (2 * k - 1) * (2 * k - 1)) / (8 * x) / k;
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Do we want to put spaces around * and / ?

This is what 'indent -kr' gave. Certainly I'd put a space between (...) / (...),
but possibly not in 2 * k.

http://codereview.appspot.com/20078/diff/1/5#newcode175
Line 175: * Computed with:
On 2009/02/26 04:13:12, charlesr.harris wrote:
> I wonder if we shouldn't consider sympy or something with arbitrary precision
> for these sort of things so that they can be extended beyond double precision.
> Not complaining, just suggesting.

I think we could then just use Numpy's longdoubles. But it's probably not a big
trouble to rewrite this in Sympy/mpmath.

http://codereview.appspot.com/20078/diff/1/5#newcode398
Line 398: /* |x| <= |v|, CF1_ik converges rapidly
On 2009/02/26 04:13:12, charlesr.harris wrote:
> See above on multiline comments.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode445
Line 445: BOOST_ASSERT(fabs(x) > 1);
On 2009/02/26 04:13:12, charlesr.harris wrote:
> BOOST_ASSERT is nice, but boost specific. Is it defined somewhere?

Defined in this file, left in place since it's boost code.

http://codereview.appspot.com/20078/diff/1/5#newcode490
Line 490: need_i = 1,
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Ummm.. Ok. Why not hex as these are used as flags? Maybe a comment to that
> effect also. Static const int need_i = 0x1 might be an alternative.

Done.
Sign in to reply to this message.

Powered by Google App Engine
RSS Feeds Recent Issues | This issue
This is Rietveld 1278:e6ce13d99bf5