[SciPy-dev] Numerical Recipes (was tagging 0.7rc1 this weekend?)

Jarrod Millman millman@berkeley....
Tue Dec 16 00:59:51 CST 2008


We are in better shape regarding the NR issue.  3 out of the 7 issues
have been resolved.  Before branching and tagging the release
candidate for 0.7.0, we have to deal with the remaining 4 issues.

Unfortunately, I don't have time to even look into this right now.  We
are really close and I hope someone will pitch in and resolve the
remaining issues.  Unless someone can say I wrote this code and I
didn't derive it from the NR code.  Since no one has stepped up to say
this, I assume we are going to have to review the code and make a
determination about this.  Basically, someone needs compare our code
with the NR code.  Check to see if the code is identical.  If the code
isn't identical, does it look like it is the NR code with
modifications (e.g., the same call signature, the same loops, the same
variable names, etc.).  If the code is derived from NR, we need to
remove it and if we want to include the functionality someone would
need to rewrite the code from scratch.  If our code isn't derived from
NR code, we need to correctly document this fact so that we don't
revisit this issue in a year.  Despite how annoying this may be, the
NR copyright specifically states that we don't have the right to port
their code to another language.  So having taken their code and
rewritten it in Python isn't allowed.  We are allowed to read their
descriptions of the algorithms and write code based on that
description; we just can't use their code to write the code.

Here is a summary of where I think we are now (please let me know if I
am incorrect):

1.  2 functions (ridder, brentq) in scipy/optimize/zeros.py

Resolved.

2.  quickselect in scipy/signal/medianfilter.c

Unresolved.

3.  scipy/special/orthogonal.py

Unresolved.

4.  3 functions in scipy/stats/stats.py.

Partially resolved.  One down (Josef rewrote ks_2samp), two left
(ttest_ind and ttest_rel).

On Wed, Dec 3, 2008 at 12:35 PM, Jarrod Millman <millman@berkeley.edu> wrote:
> Given the discussion about Numerical Recipes derived-code, I did a
> little grepping and found 7 instances where code is referencing
> Numerical Recipes.  Before branching for 0.7, I would like to see this
> resolved.  I am not sure what we should do, but the license for using
> Numerical Recipes derived code clearly does not permit us to include
> it in scipy:
>  http://www.nr.com/licenses/redistribute.html
>
> Below I include all the instances I could easily find.  Thoughts?
>
> --------------------------------------------------
>
> In scipy/optimize/zeros.py:
>
> def ridder(f, a, b, args=(),
>           xtol=_xtol, rtol=_rtol, maxiter=_iter,
>           full_output=False, disp=False):
>    """Find root of f in [a,b]
>
>    Ridder routine to find a zero of the function
>    f between the arguments a and b. f(a) and f(b) can not
>    have the same signs. Faster than bisection, but not
>    generaly as fast as the brent rountines. A description
>    may be found in a recent edition of Numerical Recipes.
>    The routine here is modified a bit to be more careful
>    of tolerance.
>
> def brentq(f, a, b, args=(),
>           xtol=_xtol, rtol=_rtol, maxiter=_iter,
>           full_output=False, disp=False):
>    """Find root of f in [a,b]
>
>    The classic Brent routine to find a zero of the function
>    f between the arguments a and b. f(a) and f(b) can not
>    have the same signs. Generally the best of the routines here.
>    It is a safe version of the secant method that uses inverse
>    quadratic extrapolation. The version here is a slight
>    modification that uses a different formula in the extrapolation
>    step. A description may be found in Numerical Recipes, but the
>    code here is probably easier to understand.
>
> In scipy/signal/medianfilter.c:
>
> /*
>  *  This Quickselect routine is based on the algorithm described in
>  *  "Numerical recipies in C", Second Edition,
>  *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
>  */
>
> In  scipy/special/orthogonal.py:
>
> """
> A collection of functions to find the weights and abscissas for
> Gaussian Quadrature.
>
> These calculations are done by finding the eigenvalues of a
> tridiagonal matrix whose entries are dependent on the coefficients
> in the recursion formula for the orthogonal polynomials with the
> corresponding weighting function over the interval.
>
> Many recursion relations for orthogonal polynomials are given:
>
> a1n f_n+1 (x) = (a2n + a3n x ) f_n (x) - a4n f_n-1 (x)
>
> The recursion relation of interest is
>
> P_n+1 (x) = (x - A_n) P_n (x) - B_n P_n-1 (x)
>
> where P has a different normalization than f.
>
> The coefficients can be found as:
>
> A_n = -a2n / a3n
>
> B_n = ( a4n / a3n sqrt(h_n-1 / h_n))**2
>
>     where
>             h_n = int_a^b w(x) f_n(x)^2
> assume:
> P_0(x) = 1
> P_-1(x) == 0
>
> See Numerical Recipies in C, page 156 and
> Abramowitz and Stegun p. 774, 782
>
> In scipy/stats/stats.py:
>
> def ttest_ind(a, b, axis=0):
>    """Calculates the t-obtained T-test on TWO INDEPENDENT samples of scores
>    a, and b.  From Numerical Recipies, p.483. Axis can equal None (ravel
>    array first), or an integer (the axis over which to operate on a and b).
>
> def ttest_rel(a,b,axis=None):
>    """Calculates the t-obtained T-test on TWO RELATED samples of scores, a
>    and b.  From Numerical Recipies, p.483. Axis can equal None (ravel array
>    first), or an integer (the axis over which to operate on a and b).
>
> def ks_2samp(data1, data2):
>    """ Computes the Kolmogorov-Smirnof statistic on 2 samples.  Modified
>    from Numerical Recipies in C, page 493.  Returns KS D-value, prob.  Not
>    ufunc- like.


More information about the Scipy-dev mailing list