Skip to content

Commit

Permalink
BF: new numpy protects polynomial coefficients
Browse files Browse the repository at this point in the history
New numpy returns a copy of the polynomial coefficients from poly.c,
rather than the actual coefficent array, so when we were modifying the
polynomial coefficients, the polynomial object wasn't seeing our
changes, leading to errors in the RFT module - see e.g.

https://nipy.bic.berkeley.edu/builders/nipy-py2.7-osx-10.10/builds/52/steps/shell_9/logs/stdio

Numpy issue raised on numpy-discussion mailing list.
  • Loading branch information
matthew-brett committed Jul 7, 2017
1 parent d49e829 commit c35c0d9
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 10 deletions.
14 changes: 7 additions & 7 deletions nipy/algorithms/statistics/rft.py
Expand Up @@ -87,14 +87,14 @@ def Q(dim, dfd=np.inf):
j = dim
if j <= 0:
raise ValueError('Q defined only for dim > 0')
poly = hermitenorm(j-1)
poly = np.poly1d(np.around(poly.c))
coeffs = np.around(hermitenorm(j - 1).c)
if np.isfinite(m):
for l in range((j-1)//2+1):
f = np.exp(gammaln((m+1)/2.) - gammaln((m+2-j+2*l)/2.)
- 0.5*(j-1-2*l)*(np.log(m/2.)))
poly.c[2*l] *= f
return np.poly1d(poly.c)
for L in range((j - 1) // 2 + 1):
f = np.exp(gammaln((m + 1) / 2.)
- gammaln((m + 2 - j + 2 * L) / 2.)
- 0.5 * (j - 1 - 2 * L) * (np.log(m / 2.)))
coeffs[2 * L] *= f
return np.poly1d(coeffs)


class ECquasi(np.poly1d):
Expand Down
6 changes: 3 additions & 3 deletions nipy/algorithms/statistics/tests/test_rft.py
Expand Up @@ -400,16 +400,16 @@ def test_hotelling2():
chi = rft.ChiSquared(dfn=dfn)(x)
assert_almost_equal(h, chi)
chi2 = scipy.stats.chi2.sf(x, dfn)
yield assert_almost_equal, h, chi2
assert_almost_equal(h, chi2)
# XXX - p appears to be unused
p = rft.spherical_search(dfn)
for dfd in [40,50]:
fac = (dfd-dfn+1.)/(dfd*dfn)
h = rft.Hotelling(dfd=dfd,k=dfn)(x)
f = scipy.stats.f.sf(x*fac, dfn, dfd-dfn+1)
f2 = rft.FStat(dfd=dfd-dfn+1,dfn=dfn)(x*fac)
yield assert_almost_equal, f2, f
yield assert_almost_equal, h, f
assert_almost_equal(f2, f)
assert_almost_equal(h, f)


@dec.slow
Expand Down

0 comments on commit c35c0d9

Please sign in to comment.