Numerical Recipes (was tagging 0.7rc1 this weekend?)
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.
On Wed, Dec 3, 2008 at 14:35, 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?
There's a difference between referencing NR's descriptions of algorithms and translating its code. I suspect that the functions in stats.py, originally written by Gary Perlman, probably fall into the latter category, but the others probably don't. Before pulling out functions, it would be worthwhile to compare our implementations with NR's to see if they in fact do look translated. But do keep in mind that the algorithms are not copyrightable, and the fact that there are only so many ways to implement said algorithms correctly is fairly substantial protection for us. IANAL. TINLA. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On Wed, Dec 3, 2008 at 12:48 PM, Robert Kern <robert.kern@gmail.com> wrote:
There's a difference between referencing NR's descriptions of algorithms and translating its code. I suspect that the functions in stats.py, originally written by Gary Perlman, probably fall into the latter category, but the others probably don't. Before pulling out functions, it would be worthwhile to compare our implementations with NR's to see if they in fact do look translated.
But do keep in mind that the algorithms are not copyrightable, and the fact that there are only so many ways to implement said algorithms correctly is fairly substantial protection for us.
Agreed. I don't necessarily think that we need to rip out the functions nor simply remove the reference to Numerical Recipes. If it turns out that we all ready have (or could easily have with a quick rewrite) code that is based on the algorithms described in Numerical Recipes (but not *derived* from their code), we should leave them in and make sure to clarify in the comments that this is the case. -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
On Wed, Dec 3, 2008 at 1:53 PM, Jarrod Millman <millman@berkeley.edu> wrote:
On Wed, Dec 3, 2008 at 12:48 PM, Robert Kern <robert.kern@gmail.com> wrote:
There's a difference between referencing NR's descriptions of algorithms and translating its code. I suspect that the functions in stats.py, originally written by Gary Perlman, probably fall into the latter category, but the others probably don't. Before pulling out functions, it would be worthwhile to compare our implementations with NR's to see if they in fact do look translated.
But do keep in mind that the algorithms are not copyrightable, and the fact that there are only so many ways to implement said algorithms correctly is fairly substantial protection for us.
Agreed. I don't necessarily think that we need to rip out the functions nor simply remove the reference to Numerical Recipes. If it turns out that we all ready have (or could easily have with a quick rewrite) code that is based on the algorithms described in Numerical Recipes (but not *derived* from their code), we should leave them in and make sure to clarify in the comments that this is the case.
Especially since much of the code in NR is often a slight rewrite of other code. For instance, Brent's original was written in Algol-68 and published in a book among other places. Moler did a nice Fortran version with lots of gotos, still the easiest to understand IMHO. I did structured version myself in the 80's, a poor fit, actually, as the algorithm is best regarded as a finite state machine. NR uses a structured version but I have no idea where it comes from; their formulas for inverse quadratic extrapolation are Brent's and not the versions I used, etc., etc.... These things really have to be checked on a case by case basis and judgement brought to bear. There is very little that is original in NR except the collection into a book with explanatory text. That's a major accomplishment but it doesn't render the code theirs. IIRC, their code for such things as computing elliptic integrals looked very much "copied" in early editions. Chuck
As for `ridder` and `brentq`, I hardly think NR would have a hope of copyrighting anything but their actual code implementing these. A classic reference should be given in any case: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a Real Continuous Function." IEEE Trans. Circuits Systems 26, 979-980, 1979. Alan Isaac
I looked at ttest_rel and ttest_ind, both are simple t-tests for equality of two samples, either in mean or of a twice ("related") observed random variable. step 1: calculate difference between samples (means) step 2: divide by appropriate standard deviation measure step 3: look up distribution to report p-value The doc strings don't explain well what the test is useful for, but it looks a straight forward implementation of the statistical formulas. The only messy part is to make it more general so that it also applies for higher dimensional arrays, and that looks all numpy to me. ttest_ind is identical to http://en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test ttest_rel is described in http://en.wikipedia.org/wiki/T-test#Dependent_t-test So, if there is any resemblance left to NR (which I don't know) then it's purely because it calculates the same thing. I think the two Wikipedia references are a lot better than NR, since Wikipedia also explains a bit the background and usage. The only part that I wouldn't have immediately implemented, is handling of zerodivision problems when the variance (of both samples or of difference of samples) is zero, which is an unusual boundary case. Josef
josef.pktd@gmail.com wrote:
I looked at ttest_rel and ttest_ind, both are simple t-tests for equality of two samples, either in mean or of a twice ("related") observed random variable.
step 1: calculate difference between samples (means) step 2: divide by appropriate standard deviation measure step 3: look up distribution to report p-value
The doc strings don't explain well what the test is useful for, but it looks a straight forward implementation of the statistical formulas. The only messy part is to make it more general so that it also applies for higher dimensional arrays, and that looks all numpy to me.
ttest_ind is identical to http://en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test ttest_rel is described in http://en.wikipedia.org/wiki/T-test#Dependent_t-test
So, if there is any resemblance left to NR (which I don't know) then it's purely because it calculates the same thing. I think the two Wikipedia references are a lot better than NR, since Wikipedia also explains a bit the background and usage.
The only part that I wouldn't have immediately implemented, is handling of zerodivision problems when the variance (of both samples or of difference of samples) is zero, which is an unusual boundary case.
Josef _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
Hi, Could you also please replace the use of 'betai' with the appropriate call to the t-distribution to the the probabilities? Numerical Recipes uses betai and also it is more understandable to use the actual t-distribution. Bruce
Hi, Could you also please replace the use of 'betai' with the appropriate call to the t-distribution to the the probabilities? Numerical Recipes uses betai and also it is more understandable to use the actual t-distribution.
Bruce
betai is used in several places in scipy.stats, until now I never touched scipy.special calls unless there was a real bug. There is some duplication in scipy.special, but I have no idea about the relative advantages and disadvantages of the different implementations. For some functions, I know that they don't work correctly over some parameter range. So, until now I'm reluctant to change any of the calls to special. Usually I don't change anything without having some tests first, and neither ttest_ind nor ttest_rand have any tests in the test suite and I'm running out of time to write tests before 0.70. but this replacement looks simple enough:
df=5;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2 0.6382988716409288 0.63829887164092902 df=500;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2 0.61729505009465102 0.61729505009466568 df=5000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2 0.61709708085302761 0.61709708085403392 df=50000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2 0.61707727784367095 0.61707727785345856 df=2;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2 0.66666666666666596 0.66666666666666652
I will make the change tonight. Josef
josef.pktd@gmail.com wrote:
Hi, Could you also please replace the use of 'betai' with the appropriate call to the t-distribution to the the probabilities? Numerical Recipes uses betai and also it is more understandable to use the actual t-distribution.
Bruce
betai is used in several places in scipy.stats, until now I never touched scipy.special calls unless there was a real bug.
There is some duplication in scipy.special, but I have no idea about the relative advantages and disadvantages of the different implementations. For some functions, I know that they don't work correctly over some parameter range. So, until now I'm reluctant to change any of the calls to special.
Usually I don't change anything without having some tests first, and neither ttest_ind nor ttest_rand have any tests in the test suite and I'm running out of time to write tests before 0.70.
but this replacement looks simple enough:
df=5;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
0.6382988716409288 0.63829887164092902
df=500;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
0.61729505009465102 0.61729505009466568
df=5000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
0.61709708085302761 0.61709708085403392
df=50000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
0.61707727784367095 0.61707727785345856
df=2;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
0.66666666666666596 0.66666666666666652
I will make the change tonight.
Josef _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
The probability density function can be written as an incomplete beta function see for example: http://en.wikipedia.org/wiki/Student%27s_t-distribution So it is probably just numerical error between division of two gamma functions and using the inverse of a beta function. While not an expert on this, it looks a like numerical precision of a double since the difference starts around 14 decimal place. Bruce
On Wed, Dec 3, 2008 at 1:59 PM, Alan G Isaac <aisaac@american.edu> wrote:
As for `ridder` and `brentq`, I hardly think NR would have a hope of copyrighting anything but their actual code implementing these. A classic reference should be given in any case:
Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a Real Continuous Function." IEEE Trans. Circuits Systems 26, 979-980, 1979.
Hello Alan, Would you be willing to take care of `ridder` and `brentq`? If you are willing, I would appreciate it if you could verify that our implementations are not derived from the Numerical Recipes code. The docstrings for these two functions seem to indicate that the code differs from that in Numerical Recipes, so you may be able to find references for our implementations. Hopefully you could do somethings like Josef did for ks2_samp: http://projects.scipy.org/scipy/scipy/changeset/5223 Since we need to make sure we aren't violating the Numerical Recipes' license (now that we have become aware of it) before releasing 0.7.0, I hope that you will be able to take care of this ASAP. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
On Fri, Dec 5, 2008 at 9:18 PM, Jarrod Millman <millman@berkeley.edu> wrote:
On Wed, Dec 3, 2008 at 1:59 PM, Alan G Isaac <aisaac@american.edu> wrote:
As for `ridder` and `brentq`, I hardly think NR would have a hope of copyrighting anything but their actual code implementing these. A classic reference should be given in any case:
Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a Real Continuous Function." IEEE Trans. Circuits Systems 26, 979-980, 1979.
Hello Alan,
Would you be willing to take care of `ridder` and `brentq`? If you are willing, I would appreciate it if you could verify that our implementations are not derived from the Numerical Recipes code. The docstrings for these two functions seem to indicate that the code differs from that in Numerical Recipes, so you may be able to find references for our implementations.
Hopefully you could do somethings like Josef did for ks2_samp: http://projects.scipy.org/scipy/scipy/changeset/5223
Since we need to make sure we aren't violating the Numerical Recipes' license (now that we have become aware of it) before releasing 0.7.0, I hope that you will be able to take care of this ASAP.
Hi Jarrod, I wrote the routines. They are similar too, but not the same as the NR versions. Chuck
On Fri, Dec 5, 2008 at 9:18 PM, Jarrod Millman <millman@berkeley.edu wrote: Would you be willing to take care of `ridder` and `brentq`? If you are willing, I would appreciate it if you could verify that our implementations are not derived from the Numerical Recipes code. The docstrings for these two functions seem to indicate that the code differs from that in Numerical Recipes, so you may be able to find references for our implementations. Hopefully you could do somethings like Josef did for ks2_samp: http://projects.scipy.org/scipy/scipy/changeset/5223
On 12/6/2008 3:19 AM Charles R Harris apparently wrote:
I wrote the routines. They are similar too, but not the same as the NR versions.
So I take it then that is needed is a docstring revision? I have attempted this for `ridder`: http://docs.scipy.org/scipy/docs/scipy.optimize.zeros.ridder/ (Note that the parameter `disp` remains undocumented as I do not know its function.) If this is what is desired, I am happy to change the brentq documentation similarly. (http://docs.scipy.org/scipy/docs/scipy.optimize.zeros.brentq/) Cheers, Alan
On 12/6/2008 11:14 AM Alan G Isaac apparently wrote:
I am happy to change the brentq documentation similarly. (http://docs.scipy.org/scipy/docs/scipy.optimize.zeros.brentq/)
Well anyway, I did this. Let me know if other action is desired. Alan
On Sat, Dec 6, 2008 at 12:19 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
I wrote the routines. They are similar too, but not the same as the NR versions.
Hey Chuck, That is good to know. It is fine if the code is similar as long as it isn't derived from their code. I assume that the similarity is do to the fact that the two codes are using similar algorithms. Could you confirm that this is indeed the case? Also could you review the docstring edits that Alan made: http://docs.scipy.org/scipy/docs/scipy.optimize.zeros.ridder/ http://docs.scipy.org/scipy/docs/scipy.optimize.zeros.brentq/ Finally, Alan wasn't able to document the parameter `disp` for ridder. Could you please provide Alan with an explanation for `disp`? Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
On 12/7/2008 5:08 AM Jarrod Millman apparently wrote:
Finally, Alan wasn't able to document the parameter `disp` for ridder. Could you please provide Alan with an explanation for `disp`?
I would think it would be: disp : number non-zero to print convergence messages. But in testing I did not see any effect ... Alan
Sun, 07 Dec 2008 09:11:03 -0500, Alan G Isaac wrote:
On 12/7/2008 5:08 AM Jarrod Millman apparently wrote:
Finally, Alan wasn't able to document the parameter `disp` for ridder. Could you please provide Alan with an explanation for `disp`?
I would think it would be:
disp : number non-zero to print convergence messages.
But in testing I did not see any effect ...
Looking at the source code, output is printed only if it fails to converge. -- Pauli Virtanen
On Sun, Dec 7, 2008 at 10:40 AM, Pauli Virtanen <pav@iki.fi> wrote:
Sun, 07 Dec 2008 09:11:03 -0500, Alan G Isaac wrote:
On 12/7/2008 5:08 AM Jarrod Millman apparently wrote:
Finally, Alan wasn't able to document the parameter `disp` for ridder. Could you please provide Alan with an explanation for `disp`?
I would think it would be:
disp : number non-zero to print convergence messages.
But in testing I did not see any effect ...
Looking at the source code, output is printed only if it fails to converge.
I'll fix it to raise an error. The question is whether is should go into the current release... Chuck
On Sun, Dec 7, 2008 at 12:44 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
I'll fix it to raise an error. The question is whether is should go into the current release...
Yes, please update the trunk so that this will make the upcoming release. It would be good if you could also add a test. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
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.
On Mon, Dec 15, 2008 at 11:59 PM, Jarrod Millman <millman@berkeley.edu>wrote:
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
The first mention I see of this method of calculating the sample points and weights for Gaussian quadrature is in this 1967 technical report: ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/67/81/CS-TR-67-81.pdf. It doesn't look to me like the NR version of the Golub-Welsch algorithm, or numpy's for that matter, differs in any significant way. Note that the original either misspells Welsch or a citation mutation occurred down the line and the name has been spelled Welsh by all those who copied it without looking up the original. A more official publication is G.H. Golub and J.H.Welsch, Computation of Gauss quadrature rules, Math. Comp., 23(1969), pp. 221-230 Chuck
On 12/16/2008 1:59 AM Jarrod Millman apparently wrote:
3. scipy/special/orthogonal.py
Unresolved.
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
Travis O. is listed as writing this orthogonal.py in 2000 and doing bug fixes in 2003. Note that the reference to Numerical Recipes is one of two references documenting an algebraic recursion relation and is not referring to an implementation. http://www.fizyka.umk.pl/nrbook/c4-5.pdf http://www-lsp.ujf-grenoble.fr/recherche/a3t2/a3t2a2/bahram/biblio/abramowit... I do not see any problem with the code. I think this one should also be considered resolved. Alan
Howdy, On Tue, Dec 16, 2008 at 7:34 AM, Alan G Isaac <aisaac@american.edu> wrote:
See Numerical Recipies in C, page 156 and Abramowitz and Stegun p. 774, 782
Travis O. is listed as writing this orthogonal.py in 2000 and doing bug fixes in 2003. Note that the reference to Numerical Recipes is one of two references documenting an algebraic recursion relation and is not referring to an implementation.
http://www.fizyka.umk.pl/nrbook/c4-5.pdf
http://www-lsp.ujf-grenoble.fr/recherche/a3t2/a3t2a2/bahram/biblio/abramowit...
I do not see any problem with the code. I think this one should also be considered resolved.
I concur. I also had a look at the Golub&Welsch reference pointed to by Chuck, and I don't see any problems. The basic algorithm may be the same, but the actual implementation in those chapters of NR is fairly different from what's in special/orthogonal.py. The Scipy code is structured as a bunch of routines that provide quite a few different orthogonal polyinomial objects, while NR only lists code for a few specific cases where the weights and abscissas are being computed (gauleg, gaulag, gauher, gaujac), plus the gaucof() that implements the Golub-Welsch algorithm for generic abscissas and weights. In orthogonal.py, the closest pairing is between NR's gaucof() and gen_roots_and_weights(), but their actual codes are fairly different to the extent possible. That is, it's a short code that does a very specific thing, but the scipy code does it in a very python/numpy specific way, while the NR code is very 'C'-ish. Other than the fact that these codes both use the same piece of mathematics, I see no problem here whatsoever. It might be worth clarifying in the docstring that reads """ See Numerical Recipies in C, page 156 and Abramowitz and Stegun p. 774, 782 """ that the NR reference is to the formulas only, not to any actual implementation. If nothing else, to prevent any questions in the future... Cheers, f
On Tue, Dec 16, 2008 at 7:34 AM, Alan G Isaac <aisaac@american.edu> wrote:
Travis O. is listed as writing this orthogonal.py in 2000 and doing bug fixes in 2003. Note that the reference to Numerical Recipes is one of two references documenting an algebraic recursion relation and is not referring to an implementation. http://www.fizyka.umk.pl/nrbook/c4-5.pdf http://www-lsp.ujf-grenoble.fr/recherche/a3t2/a3t2a2/bahram/biblio/abramowit... I do not see any problem with the code.
Fernando Perez wrote:
Other than the fact that these codes both use the same piece of mathematics, I see no problem here whatsoever. It might be worth clarifying in the docstring
Am I correct that the documentation project does not give access to the module-level docstrings? I would have been happy to make this change in orthogonal.py but did not see how to add documentation via the project. Anyway, here is the reST for the GW cite and what I think the AS citation info is. (Not sure about the address.) fwiw, Alan Isaac For the mathematical background, see [golub.welsch-1969-mathcomp]_ and [abramowitz.stegun-1965]_. .. [golub.welsch-1969-mathcomp] Golub, Gene H, and John H Welsch. 1969. Calculation of Gauss Quadrature Rules. *Mathematics of Computation* 23, 221-230+s1--s10. .. [abramowitz.stegun-1965] Abramowitz, Milton, and Irene A Stegun. (1965) *Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables*. Gaithersburg, MD: National Bureau of Standards. http://www.math.sfu.ca/~cbm/aands/
Hi Alan, On Sun, Jan 4, 2009 at 2:20 PM, Alan G Isaac <aisaac@american.edu> wrote:
Anyway, here is the reST for the GW cite and what I think the AS citation info is. (Not sure about the address.)
I just took the liberty of committing this to trunk so it doesn't fall through the cracks, with due credit. Jarrod will backport it now to the release branch. Thanks! f
Sun, 04 Jan 2009 17:20:00 -0500, Alan G Isaac wrote:
On Tue, Dec 16, 2008 at 7:34 AM, Alan G Isaac <aisaac@american.edu> wrote:
Travis O. is listed as writing this orthogonal.py in 2000 and doing bug fixes in 2003. Note that the reference to Numerical Recipes is one of two references documenting an algebraic recursion relation and is not referring to an implementation. http://www.fizyka.umk.pl/nrbook/c4-5.pdf http://www-lsp.ujf-grenoble.fr/recherche/a3t2/a3t2a2/bahram/biblio/ abramowitz_and_stegun/page_774.htm I do not see any problem with the code.
Fernando Perez wrote:
Other than the fact that these codes both use the same piece of mathematics, I see no problem here whatsoever. It might be worth clarifying in the docstring
Am I correct that the documentation project does not give access to the module-level docstrings? I would have been happy to make this change in orthogonal.py but did not see how to add documentation via the project.
It does, here: http://docs.scipy.org/scipy/docs/scipy.special.orthogonal/ -- Pauli Virtanen
On 1/11/2009 12:42 PM Pauli Virtanen apparently wrote:
OK, thanks. I cleaned that up a bit. Alan
On Tue, Dec 16, 2008 at 01:59, Jarrod Millman <millman@berkeley.edu> wrote:
2. quickselect in scipy/signal/medianfilter.c
Unresolved.
I looked at the NR code for select (2nd ed, page 342). It was perhaps a little too close for comfort, but someone else might disagree. I rewrote the select algorithm(s) from scratch, referring to the Wikipedia page describing quickselect, and tested it against the existing functions on a large number of random cases (small and large buffers, all the same values, etc.). It is perhaps a little faster than the current version. It's also written as a single macro that expands for the float/double/byte cases, rather than replicating the code. Someone should review this code, and if it looks good, I'll submit a patch for medianfilter.c. For 0.8, perhaps scipy should use a variant of constant time median filtering (http://nomis80.org/ctmf.html) Ray Jones
On Tue, Dec 16, 2008 at 12:41 PM, Thouis (Ray) Jones <thouis@broad.mit.edu> wrote:
On Tue, Dec 16, 2008 at 01:59, Jarrod Millman <millman@berkeley.edu> wrote:
2. quickselect in scipy/signal/medianfilter.c
Unresolved.
I looked at the NR code for select (2nd ed, page 342). It was perhaps a little too close for comfort, but someone else might disagree.
I rewrote the select algorithm(s) from scratch, referring to the Wikipedia page describing quickselect, and tested it against the existing functions on a large number of random cases (small and large buffers, all the same values, etc.). It is perhaps a little faster than the current version. It's also written as a single macro that expands for the float/double/byte cases, rather than replicating the code.
Thanks very much. This was the code I was most concerned with as well. Hopefully, someone will have a chance to review your code over the next few days. -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
Here's a slightly better version (clearer, about 10% faster). Ray Jones
On Tue, Dec 16, 2008 at 6:50 PM, Thouis (Ray) Jones <thouis@broad.mit.edu> wrote:
Here's a slightly better version (clearer, about 10% faster).
It looks good to me. Did you take a look at whether medfilter2 was from NR as well? If not, I would appreciate it if you would take a quick look. It would also be nice if you would be willing to rewrite the medfilter2s as a single macro which expands for the float/double/byte cases like you did for the quick select. Could you also provide some tests for medfilt2 in test_signaltools.py: trunk/scipy/signal/tests/test_signaltools.py Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
On Wed, Dec 17, 2008 at 16:25, Jarrod Millman <millman@berkeley.edu> wrote:
On Tue, Dec 16, 2008 at 6:50 PM, Thouis (Ray) Jones <thouis@broad.mit.edu> wrote:
Here's a slightly better version (clearer, about 10% faster).
It looks good to me. Did you take a look at whether medfilter2 was from NR as well? If not, I would appreciate it if you would take a quick look. It would also be nice if you would be willing to rewrite the medfilter2s as a single macro which expands for the float/double/byte cases like you did for the quick select.
Could you also provide some tests for medfilt2 in test_signaltools.py: trunk/scipy/signal/tests/test_signaltools.py
The medfilter2 is not from NR (NR doesn't include median filtering, as far as I can tell). I've attached the patch against the latest SVN. It includes the following: * new quickselect (from scratch) in medianfilter.c * new macros for {f,d,b}_medfilt2 in medianfilter.c * modified test for medfilt and medfilt2 (larger array, non-square filter). The N-d order filter (called by medfilt) works by filling a buffer, then calling qsort, then extracting the desired index. It could easily take advantage of a qselect() function based on the one in medianfilter.c (except without . I'll file an enhancement request, as I doubt you want such a change in 0.7. Also, medfilt2d should fall back to medfilt (nD version) when handed a type not in its special cases (float, double, unsigned char). Ray Jones
On Thu, Dec 18, 2008 at 9:20 AM, Thouis (Ray) Jones <thouis@broad.mit.edu> wrote:
The medfilter2 is not from NR (NR doesn't include median filtering, as far as I can tell).
I've attached the patch against the latest SVN. It includes the following: * new quickselect (from scratch) in medianfilter.c * new macros for {f,d,b}_medfilt2 in medianfilter.c * modified test for medfilt and medfilt2 (larger array, non-square filter).
Thank you very much for looking into this and supplying a patch. I went ahead and applied the patch: http://projects.scipy.org/scipy/scipy/changeset/5279 Thanks for the enhancement tickets as well: http://projects.scipy.org/scipy/scipy/ticket/818 http://projects.scipy.org/scipy/scipy/ticket/819 -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/
2008/12/16 Jarrod Millman <millman@berkeley.edu>:
4. 3 functions in scipy/stats/stats.py.
Partially resolved. One down (Josef rewrote ks_2samp), two left (ttest_ind and ttest_rel).
From the NR license:
"""Copyright does not protect ideas, but only the expression of those ideas in a particular form. In the case of a computer program, the ideas consist of the program's methodology and algorithm, including the necessary sequence of steps adopted by the programmer. The expression of those ideas is the program source code (particularly any arbitrary or stylistic choices embodied in it), its derived object code, and any other derivative work. If you analyze the ideas contained in a program, and then express those ideas in your own completely different implementation, then that new program implementation belongs to you. That is what we have done for those programs in this book that are not entirely of our own devising. When programs in this book are said to be "based" on programs published in copyright sources, we mean that the ideas are the same. The expression of these ideas as source code is our own. We believe that no material in this book infringes on an existing copyright.""" I am not convinced that their interpretation of copyright law is correct, but either way I think ttest_ind and ttest_rel are safe (according to *their* rules). The code makes use of constructs such as "where" and broadcasting, that aren't available in C. The author therefore took an idea from NR and reimplemented it in a novel way. Regards Stéfan
On Sat, 20 Dec 2008 00:41:20 +0200 "Stéfan van der Walt" <stefan@sun.ac.za> wrote:
2008/12/16 Jarrod Millman <millman@berkeley.edu>:
4. 3 functions in scipy/stats/stats.py.
Partially resolved. One down (Josef rewrote ks_2samp), two left (ttest_ind and ttest_rel).
From the NR license:
"""Copyright does not protect ideas, but only the expression of those ideas in a particular form. In the case of a computer program, the ideas consist of the program's methodology and algorithm, including the necessary sequence of steps adopted by the programmer. The expression of those ideas is the program source code (particularly any arbitrary or stylistic choices embodied in it), its derived object code, and any other derivative work.
I am not convinced that their interpretation of copyright law is correct, but either way I think ttest_ind and ttest_rel are safe (according to *their* rules). The code makes use of constructs such as "where" and broadcasting, that aren't available in C. The author therefore took an idea from NR and reimplemented it in a novel way.
Most of what I know about copyright law is : - it is complex and changes depending on the object being discussed - publishers like to claim many more rights than they actually have For example, for cookbooks, the recipe itself is explicitly not copyrightable, only the layout in a book and the description. There is also a paper I read a few years ago about music copyright showing examples of publishers trying to claim copyright for Bach sonatas! I would not be surprised if their interpretation of copyright law for software is substantially correct. That is probably why people try to patent software now (one-click anyone?) - that way they get greater (though shorter term) protection. -- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
Quoting from http://www.iusmentis.com/copyright/software/protection/ An implementation of a standard algorithm may require very little creativity and if so, is likely not protected by copyright. The function of a work can only be protected by a patent. Even when providing a largely functional implementation, creativity can still be found in things like function and variable names, code optimizations and even the comments accompanying the code. Copying such an implementation wholesale then still infringes on copyright. fwiw, Alan
On Sat, Dec 20, 2008 at 10:18 AM, Alan G Isaac <aisaac@american.edu> wrote:
Quoting from http://www.iusmentis.com/copyright/software/protection/
An implementation of a standard algorithm may require very little creativity and if so, is likely not protected by copyright. The function of a work can only be protected by a patent.
Even when providing a largely functional implementation, creativity can still be found in things like function and variable names, code optimizations and even the comments accompanying the code. Copying such an implementation wholesale then still infringes on copyright.
The problem, though, is that terms like "requires little creativity" are subjective. You may think the code required no creativity, but the publisher of NR may disagree. And then it takes a court to decide who is right. And at that point even if SciPy is vindicated, the SciPy community still loses because of all the lost time an money that has to go into such a fight. It's easier to just make sure you steer well clear of any potential infringement from the beginning. IANAL, but I think that especially once the issue is raised in a public forum like this -- now you've got a record that the developers were aware of a potential problem. If nothing is done about it, it might be regarded as "willful infringement". But if an honest attempt is made to replace the code, then you've at least got a public record that you were attempting to do your "due dilligence" to eliminate infringements. So, silly as it may seem, I think that Jarrod's plan of action here is correct. The NR guys may be wrong about copyright law, but being right is not a prerequisite for filing a lawsuit. Merely thinking you are is generally enough. And sometimes not even that... --bb
On Fri, Dec 19, 2008 at 8:49 PM, Bill Baxter <wbaxter@gmail.com> wrote:
On Sat, Dec 20, 2008 at 10:18 AM, Alan G Isaac <aisaac@american.edu> wrote:
Quoting from http://www.iusmentis.com/copyright/software/protection/
An implementation of a standard algorithm may require very little creativity and if so, is likely not protected by copyright. The function of a work can only be protected by a patent.
Even when providing a largely functional implementation, creativity can still be found in things like function and variable names, code optimizations and even the comments accompanying the code. Copying such an implementation wholesale then still infringes on copyright.
The problem, though, is that terms like "requires little creativity" are subjective. You may think the code required no creativity, but the publisher of NR may disagree. And then it takes a court to decide who is right. And at that point even if SciPy is vindicated, the SciPy community still loses because of all the lost time an money that has to go into such a fight. It's easier to just make sure you steer well clear of any potential infringement from the beginning.
IANAL, but I think that especially once the issue is raised in a public forum like this -- now you've got a record that the developers were aware of a potential problem. If nothing is done about it, it might be regarded as "willful infringement". But if an honest attempt is made to replace the code, then you've at least got a public record that you were attempting to do your "due dilligence" to eliminate infringements.
So, silly as it may seem, I think that Jarrod's plan of action here is correct. The NR guys may be wrong about copyright law, but being right is not a prerequisite for filing a lawsuit. Merely thinking you are is generally enough. And sometimes not even that...
--bb _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
In these cases, I don't see any copyrightable creativity. I replaced the use of betai with the t-distribution on recommendation of Bruce Southey. The formula for `ttest_ind` is exactly the same as wikipedia and I guess in any number of textbooks, http://en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test section: Unequal sample sizes, equal variance Variable names are not very creative and the shape and axis handling is standard numpy style. The handling of the zero div problem seems to be added later and is not in the current version of Gary Strangmans code So I don't see what can be changed here, except for cosmetic changes. def ttest_ind(a, b, axis=0): a, b, axis = _chk2_asarray(a, b, axis) x1 = mean(a,axis) x2 = mean(b,axis) v1 = var(a,axis) v2 = var(b,axis) n1 = a.shape[axis] n2 = b.shape[axis] df = n1+n2-2 svar = ((n1-1)*v1+(n2-1)*v2) / float(df) zerodivproblem = svar == 0 t = (x1-x2)/np.sqrt(svar*(1.0/n1 + 1.0/n2)) # N-D COMPUTATION HERE!!!!!! t = np.where(zerodivproblem, 1.0, t) # replace NaN t-values with 1.0 #probs = betai(0.5*df,0.5,float(df)/(df+t*t)) probs = distributions.t.sf(np.abs(t),df)*2 if not np.isscalar(t): probs = probs.reshape(t.shape) if not np.isscalar(probs) and len(probs) == 1: probs = probs[0] return t, probs The second test, looks a little bit more complicated in the current version, for example it defines several variables that are never used. I am currently rewriting the main formula to follow exactly the wikipedia formula, which is much easier to read. For a Monte Carlo with 10000 replication and sample size 500, I get at maximum difference of the new version to the old version of 3.9968028886505635e-015 which should be close enough. Rejection rates for alpha = 1%, 5% and 10% are pretty close, so the test is correctly sized. This is the new code: def ttest_rel(a,b,axis=None): a, b, axis = _chk2_asarray(a, b, axis) if len(a)!=len(b): raise ValueError, 'unequal length arrays' n = a.shape[axis] df = float(n-1) d = (a-b).astype('d') #denom is just sqrt(var(d)/df), Note: denominator for denom is 1/n * 1/(n-1) denom = np.sqrt(np.var(d,axis,ddof=1)/float(n)) zerodivproblem = denom == 0 t = np.mean(d, axis) / denom t = np.where(zerodivproblem, 1.0, t) # replace NaN t-values with 1.0 probs = distributions.t.sf(np.abs(t),df)*2 if not np.isscalar(t): probs = np.reshape(probs, t.shape) if not np.isscalar(probs) and len(probs) == 1: probs = probs[0] return t, probs As an aside, is there any difference between np.sum() and np.add.reduce(), or are they equivalent or exactly the same? Josef
Bill Baxter wrote:
So, silly as it may seem, I think that Jarrod's plan of action here is correct. The NR guys may be wrong about copyright law, but being right is not a prerequisite for filing a lawsuit. Merely thinking you are is generally enough. And sometimes not even that...
I think there is a very simple argument to be made here: we are all (hopefully) better at programming than at understanding law. So we should rather spent our time reimplementing the code rather than trying to understand how the various countries implement copyrights :) David
2008/12/16 Jarrod Millman <millman@berkeley.edu>:
I am not convinced that their interpretation of copyright law is correct, but either way I think ttest_ind and ttest_rel are safe (according to *their* rules). The code makes use of constructs such as "where" and broadcasting, that aren't available in C. The author therefore took an idea from NR and reimplemented it in a novel way.
I haven't looked at the actual code in question here. Just beware that there are Fortran 90 code on the NR CD-ROM and web site. In 1996 the "NR in Fortran 90" book was published. Their Fortran 90 code does make use of where constructs, broadcasting and array slices. Sturla Molden
I am not convinced that their interpretation of copyright law is correct, but either way I think ttest_ind and ttest_rel are safe (according to *their* rules). The code makes use of constructs such as "where" and broadcasting, that aren't available in C. The author therefore took an idea from NR and reimplemented it in a novel way.
I haven't looked at the actual code in question here. Just beware that there are Fortran 90 code on the NR CD-ROM and web site. In 1996 the "NR in Fortran 90" book was published. Their Fortran 90 code does make use of where constructs, broadcasting and array slices.
Sturla Molden
Come on ;) All these algo are very standard one. You can refer to this great book to *explain* these pieces of math but that's it. If you want to code let's say a heapshort or a golden ratio based search, you may refer to one another great book to explain how is it supposed to work. Fine. That said, there is no many different ways to code such algo in fortran (it is of course the same whatever the language could be). If you think there is a need to obfuscate standard pieces of math well ....that is your choice. If someone try to claim something about maths, you can just remove all the ref to this strange guy and put a ref to any classical paper written by someone reasonable (or event to wikipedia...) Xavier
participants (15)
-
Alan G Isaac -
Alan Jackson -
Bill Baxter -
Bruce Southey -
Charles R Harris -
David Cournapeau -
Fernando Perez -
Jarrod Millman -
josef.pktd＠gmail.com -
Pauli Virtanen -
Robert Kern -
Sturla Molden -
Stéfan van der Walt -
Thouis (Ray) Jones -
Xavier Gnata