Merge pull request #249 from VirgileFritsch/master

matthew-brett · matthew-brett · commit f3323676d914 · 2013-01-22T11:57:24.000-08:00
BF: Fix bug in multiple_fast_inv

Some computations were done outside the loop, but they depend on values changing at each iteration. This resulted in wrong inverse values because of the use of the wrong piv variable. Also the check of info was done only for the last decomposition.
diff --git a/nipy/algorithms/statistics/tests/test_utils.py b/nipy/algorithms/statistics/tests/test_utils.py
@@ -3,7 +3,7 @@
 import numpy as np
 from scipy.stats import norm
 
-from ..utils import multiple_mahalanobis, z_score
+from ..utils import multiple_mahalanobis, z_score, multiple_fast_inv
 from nose.tools import assert_true
 from numpy.testing import assert_almost_equal, assert_array_almost_equal
 
@@ -31,6 +31,16 @@ def test_mahalanobis2():
     f_mah = (multiple_mahalanobis(x, Aa))[i]
     assert_true(np.allclose(mah, f_mah))
 
+def test_multiple_fast_inv():
+    shape = (10, 20, 20)
+    X = np.random.randn(shape[0], shape[1], shape[2])
+    X_inv_ref = np.zeros(shape)
+    for i in range(shape[0]):
+        X[i] = np.dot(X[i], X[i].T)
+        X_inv_ref[i] = np.linalg.inv(X[i])
+    X_inv = multiple_fast_inv(X)
+    assert_almost_equal(X_inv_ref, X_inv)
+
 
 if __name__ == "__main__":
     import nose
diff --git a/nipy/algorithms/statistics/utils.py b/nipy/algorithms/statistics/utils.py
@@ -43,32 +43,32 @@ def multiple_fast_inv(a):
     from scipy.linalg import calc_lwork
     from scipy.linalg.lapack import  get_lapack_funcs
     a1, n = a[0], a.shape[0]
-    getrf, getri = get_lapack_funcs(('getrf','getri'), (a1,))
-    if getrf.module_name[:7] == 'clapack' != getri.module_name[:7]:
-        # ATLAS 3.2.1 has getrf but not getri.
-        for i in range(n):
+    getrf, getri = get_lapack_funcs(('getrf', 'getri'), (a1,))
+    for i in range(n):
+        if (getrf.module_name[:7] == 'clapack'
+            and getri.module_name[:7] != 'clapack'):
+            # ATLAS 3.2.1 has getrf but not getri.
             lu, piv, info = getrf(np.transpose(a[i]), rowmajor=0,
                                   overwrite_a=True)
             a[i] = np.transpose(lu)
-    else:
-        for i in range(n):
+        else:
             a[i], piv, info = getrf(a[i], overwrite_a=True)
-    if info == 0:
-        if getri.module_name[:7] == 'flapack':
-            lwork = calc_lwork.getri(getri.prefix, a1.shape[0])
-            lwork = lwork[1]
-            # XXX: the following line fixes curious SEGFAULT when
-            # benchmarking 500x500 matrix inverse. This seems to
-            # be a bug in LAPACK ?getri routine because if lwork is
-            # minimal (when using lwork[0] instead of lwork[1]) then
-            # all tests pass. Further investigation is required if
-            # more such SEGFAULTs occur.
-            lwork = int(1.01 * lwork)
-            for i in range(n):
+        if info == 0:
+            if getri.module_name[:7] == 'flapack':
+                lwork = calc_lwork.getri(getri.prefix, a1.shape[0])
+                lwork = lwork[1]
+                # XXX: the following line fixes curious SEGFAULT when
+                # benchmarking 500x500 matrix inverse. This seems to
+                # be a bug in LAPACK ?getri routine because if lwork is
+                # minimal (when using lwork[0] instead of lwork[1]) then
+                # all tests pass. Further investigation is required if
+                # more such SEGFAULTs occur.
+                lwork = int(1.01 * lwork)
                 a[i], _ = getri(a[i], piv, lwork=lwork, overwrite_lu=1)
-        else: # clapack
-            for i in range(n):
+            else:  # clapack
                 a[i], _ = getri(a[i], piv, overwrite_lu=1)
+        else:
+            raise ValueError('Matrix LU decomposition failed')
     return a
 
 
