From c35c0d9e9a3a4ae4cef3a999813fb7c199acd22d Mon Sep 17 00:00:00 2001
From: Matthew Brett <matthew.brett@gmail.com>
Date: Fri, 7 Jul 2017 17:59:15 +0100
Subject: [PATCH] BF: new numpy protects polynomial coefficients

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.
---
 nipy/algorithms/statistics/rft.py            | 14 +++++++-------
 nipy/algorithms/statistics/tests/test_rft.py |  6 +++---
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/nipy/algorithms/statistics/rft.py b/nipy/algorithms/statistics/rft.py
index bb6e5af8ee..f1cd0e5160 100644
--- a/nipy/algorithms/statistics/rft.py
+++ b/nipy/algorithms/statistics/rft.py
@@ -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):
diff --git a/nipy/algorithms/statistics/tests/test_rft.py b/nipy/algorithms/statistics/tests/test_rft.py
index f52e69f36e..0cf1cb034b 100644
--- a/nipy/algorithms/statistics/tests/test_rft.py
+++ b/nipy/algorithms/statistics/tests/test_rft.py
@@ -400,7 +400,7 @@ 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]:
@@ -408,8 +408,8 @@ def test_hotelling2():
             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