Sparse Fascicle Model

.travis.yml

@@ -23,6 +23,9 @@ matrix:
       env:
         - COVERAGE=1
         - DEPENDS="cython==0.18 numpy==1.6.0 scipy==0.9.0 nibabel==1.2.0"
+    - python: 2.7
+      env:
+        - DEPENDS="cython numpy scipy matplotlib h5py nibabel cvxopt scikit_learn"
 before_install:
     - virtualenv --python=python venv
     - source venv/bin/activate

dipy/core/optimize.py

@@ -3,11 +3,13 @@
 Only L-BFGS-B and Powell is supported in this class for versions of
 Scipy < 0.12. All optimizers are available for scipy >= 0.12.
 """
-
+import abc
 from distutils.version import LooseVersion
 import numpy as np
 import scipy
 import scipy.sparse as sps
+import scipy.optimize as opt
+from dipy.utils.six import with_metaclass
 
 SCIPY_LESS_0_12 = LooseVersion(scipy.__version__) < '0.12'
 
@@ -218,7 +220,6 @@ def __init__(self, fun,  x0, args=(), method='L-BFGS-B', jac=None,
 
                 def history_of_x(kx):
                     self._evol_kx.append(kx)
-
                 res = minimize(fun, x0, args, method, jac, hess, hessp, bounds,
                                constraints, tol, callback=history_of_x,
                                options=options)
@@ -293,24 +294,6 @@ def spdot(A, B):
         return np.dot(A, B)
 
 
-def rsq(ss_residuals, ss_residuals_to_mean):
-    """
-    Calculate: $R^2 = \frac{1-SSE}{\sigma^2}$
-
-    Parameters
-    ----------
-    ss_residuals : array
-        Model fit errors relative to the data
-    ss_residuals_to_mean : array
-        Residuals of the data relative to the mean of the data (variance)
-
-    Returns
-    -------
-    rsq : the variance explained.
-    """
-    return 100 * (1 - ss_residuals/ss_residuals_to_mean)
-
-
 def sparse_nnls(y, X,
                 momentum=1,
                 step_size=0.01,
@@ -357,17 +340,14 @@ def sparse_nnls(y, X,
     h_best : The best estimate of the parameters.
 
     """
-    num_data = y.shape[0]
     num_regressors = X.shape[1]
     # Initialize the parameters at the origin:
     h = np.zeros(num_regressors)
     # If nothing good happens, we'll return that:
     h_best = h
     gradient = np.zeros(num_regressors)
     iteration = 1
-    count = 1
     ss_residuals_min = np.inf  # This will keep track of the best solution
-    ss_residuals_to_mean = np.sum((y - np.mean(y)) ** 2)  # The variance of y
     sse_best = np.inf   # This will keep track of the best performance so far
     count_bad = 0  # Number of times estimation error has gone up.
     error_checks = 0  # How many error checks have we done so far
@@ -396,21 +376,73 @@ def sparse_nnls(y, X,
             if sse < ss_residuals_min:
                 # Update your expectations about the minimum error:
                 ss_residuals_min = sse
-                n_iterations = iteration  # This holds the number of iterations
-                                          # for the best solution so far.
                 h_best = h  # This holds the best params we have so far
-
                 # Are we generally (over iterations) converging on
                 # sufficient improvement in r-squared?
                 if sse < converge_on_sse * sse_best:
                     sse_best = sse
                     count_bad = 0
                 else:
-                    count_bad +=1 
+                    count_bad += 1
             else:
                 count_bad += 1
 
             if count_bad >= max_error_checks:
                 return h_best
             error_checks += 1
         iteration += 1
+
+
+class SKLearnLinearSolver(with_metaclass(abc.ABCMeta, object)):
+    """
+    Provide a sklearn-like uniform interface to algorithms that solve problems
+    of the form: $y = Ax$ for $x$
+
+    Sub-classes of SKLearnLinearSolver should provide a 'fit' method that have
+    the following signature: `SKLearnLinearSolver.fit(X, y)`, which would set
+    an attribute `SKLearnLinearSolver.coef_`, with the shape (X.shape[1],),
+    such that an estimate of y can be calculated as:
+    `y_hat = np.dot(X, SKLearnLinearSolver.coef_.T)`
+    """
+    def __init__(self, *args, **kwargs):
+        self._args = args
+        self._kwargs = kwargs
+
+    @abc.abstractmethod
+    def fit(self, X, y):
+        """Implement for all derived classes """
+
+    def predict(self, X):
+        """
+        Predict using the result of the model
+
+        Parameters
+        ----------
+        X : array-like (n_samples, n_features)
+            Samples.
+
+        Returns
+        -------
+        C : array, shape = (n_samples,)
+            Predicted values.
+        """
+        X = np.asarray(X)
+        return np.dot(X, self.coef_.T)
+
+
+class NonNegativeLeastSquares(SKLearnLinearSolver):
+    """
+    A sklearn-like interface to scipy.optimize.nnls
+
+    """
+    def fit(self, X, y):
+        """
+        Fit the NonNegativeLeastSquares linear model to data
+
+        Parameters
+        ----------
+
+        """
+        coef, rnorm = opt.nnls(X, y)
+        self.coef_ = coef
+        return self

dipy/core/tests/test_optimize.py

@@ -1,50 +1,43 @@
 import numpy as np
 import scipy.sparse as sps
-from numpy.testing import (assert_equal,
-                           assert_almost_equal,
-                           assert_array_almost_equal,
-                           assert_array_equal,
-                           run_module_suite)
 
 import numpy.testing as npt
 from dipy.core.optimize import Optimizer, SCIPY_LESS_0_12, sparse_nnls, spdot
+import dipy.core.optimize as opt
 
 
 def func(x):
-
     return x[0]**2 + x[1]**2 + x[2]**2
 
 
 def func2(x):
-
     return x[0]**2 + 0.5 * x[1]**2 + 0.2 * x[2]**2 + 0.2 * x[3]**2
 
 
 @npt.dec.skipif(SCIPY_LESS_0_12)
 def test_optimize_new_scipy():
-
     opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='Powell')
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='L-BFGS-B',
                     options={'maxcor': 10, 'ftol': 1e-7,
                              'gtol': 1e-5, 'eps': 1e-8})
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
-    assert_equal(opt.evolution, None)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
+    npt.assert_equal(opt.evolution, None)
 
-    assert_equal(opt.evolution, None)
+    npt.assert_equal(opt.evolution, None)
 
     opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='L-BFGS-B',
                     options={'maxcor': 10, 'ftol': 1e-7,
                              'gtol': 1e-5, 'eps': 1e-8},
                     evolution=False)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt.print_summary()
 
@@ -54,14 +47,14 @@ def test_optimize_new_scipy():
                              'gtol': 1e-5, 'eps': 1e-8},
                     evolution=True)
 
-    assert_equal(opt.evolution.shape, (opt.nit, 4))
+    npt.assert_equal(opt.evolution.shape, (opt.nit, 4))
 
     opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
                     method='Powell',
                     options={'xtol': 1e-6, 'ftol': 1e-6, 'maxiter': 1e6},
                     evolution=True)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
 
 
 @npt.dec.skipif(not SCIPY_LESS_0_12)
@@ -72,80 +65,105 @@ def test_optimize_old_scipy():
                     options={'maxcor': 10, 'ftol': 1e-7,
                              'gtol': 1e-5, 'eps': 1e-8})
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
                     method='Powell',
                     options={'xtol': 1e-6, 'ftol': 1e-6, 'maxiter': 1e6},
                     evolution=True)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
 
     opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]),
                     method='L-BFGS-B',
                     options={'maxcor': 10, 'eps': 1e-8})
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]),
                     method='L-BFGS-B',
                     options=None)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
                     method='L-BFGS-B',
                     options={'gtol': 1e-7, 'ftol': 1e-7, 'maxiter': 10000})
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]), 4)
-    assert_almost_equal(opt.fopt, 0)
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]), 4)
+    npt.assert_almost_equal(opt.fopt, 0)
 
     opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
                     method='Powell',
                     options={'maxiter': 1e6},
                     evolution=True)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
 
     opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
                     method='Powell',
                     options={'maxiter': 1e6},
                     evolution=True)
 
-    assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
+    npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))
+
+
+def test_sklearn_linear_solver():
+    class SillySolver(opt.SKLearnLinearSolver):
+        def fit(self, X, y):
+            self.coef_ = np.ones(X.shape[-1])
+
+    MySillySolver = SillySolver()
+    n_samples = 100
+    n_features = 20
+    y = np.random.rand(n_samples)
+    X = np.ones((n_samples, n_features))
+    MySillySolver.fit(X, y)
+    npt.assert_equal(MySillySolver.coef_, np.ones(n_features))
+    npt.assert_equal(MySillySolver.predict(X), np.ones(n_samples) * 20)
+
+
+def test_nonnegativeleastsquares():
+    n = 100
+    X = np.eye(n)
+    beta = np.random.rand(n)
+    y = np.dot(X, beta)
+    my_nnls = opt.NonNegativeLeastSquares()
+    my_nnls.fit(X, y)
+    npt.assert_equal(my_nnls.coef_, beta)
+    npt.assert_equal(my_nnls.predict(X), y)
 
 
 def test_spdot():
     n = 100
     m = 20
     k = 10
-    A = np.random.randn(n,m)
-    B = np.random.randn(m,k)
+    A = np.random.randn(n, m)
+    B = np.random.randn(m, k)
     A_sparse = sps.csr_matrix(A)
     B_sparse = sps.csr_matrix(B)
-
     dense_dot = np.dot(A, B)
     # Try all the different variations:
-    assert_array_almost_equal(dense_dot, spdot(A_sparse, B_sparse).todense())
-    assert_array_almost_equal(dense_dot, spdot(A, B_sparse))
-    assert_array_almost_equal(dense_dot, spdot(A_sparse, B))
+    npt.assert_array_almost_equal(dense_dot,
+                                  spdot(A_sparse, B_sparse).todense())
+    npt.assert_array_almost_equal(dense_dot, spdot(A, B_sparse))
+    npt.assert_array_almost_equal(dense_dot, spdot(A_sparse, B))
 
 
-def test_nnls():
+def test_sparse_nnls():
     # Set up the regression:
     beta = np.random.rand(10)
     X = np.random.randn(1000, 10)
     y = np.dot(X, beta)
     beta_hat = sparse_nnls(y, X)
     beta_hat_sparse = sparse_nnls(y, sps.csr_matrix(X))
     # We should be able to get back the right answer for this simple case
-    assert_array_almost_equal(beta, beta_hat, decimal=1)
-    assert_array_almost_equal(beta, beta_hat_sparse, decimal=1)
+    npt.assert_array_almost_equal(beta, beta_hat, decimal=1)
+    npt.assert_array_almost_equal(beta, beta_hat_sparse, decimal=1)
 
 
 if __name__ == '__main__':
-
-    run_module_suite()
+    npt.run_module_suite()