Merge pull request #179 from ljchang/decomposition

Add decompose method to brain_data

Former-commit-id: b9b97b1

Author: ljchang

MANIFEST.in

@@ -0,0 +1 @@
+include requirements.txt

nltools/data/brain_data.py

@@ -46,7 +46,8 @@
                             glover_hrf,
                             attempt_to_import,
                             concatenate,
-                            _bootstrap_apply_func)
+                            _bootstrap_apply_func,
+                            set_decomposition_algorithm)
 from nltools.cross_validation import set_cv
 from nltools.plotting import (dist_from_hyperplane_plot,
                               scatterplot,
@@ -1426,6 +1427,30 @@ def bootstrap(self, function, n_samples=5000, save_weights=False,
         bootstrapped = Brain_Data(bootstrapped)
         return summarize_bootstrap(bootstrapped, save_weights=save_weights)
 
+    def decompose(self, algorithm='pca', n_components=None, *args, **kwargs):
+        ''' Decompose Brain_Data object
+
+        Args:
+            algorithm: (str) Algorithm to perform decomposition
+                        types=['pca','ica','nnmf','fa']
+            n_components: (int) number of components. If None then retain
+                        as many as possible.
+        Returns:
+            output: a dictionary of decomposition parameters
+        '''
+
+        out = {}
+        out['decomposition_object'] = set_decomposition_algorithm(
+                                                    algorithm=algorithm,
+                                                    n_components=n_components,
+                                                    *args, **kwargs)
+        out['decomposition_object'].fit(self.data.T)
+        out['components'] = self.empty()
+        out['components'].data = out['decomposition_object'].transform(
+                                                                self.data.T).T
+        out['weights'] = out['decomposition_object'].components_
+        return out
+
 class Groupby(object):
     def __init__(self, data, mask):
 

nltools/tests/test_data.py

@@ -239,6 +239,42 @@ def test_brain_data_3mm(tmpdir):
     diff = m2-m1
     assert np.sum(diff.data) == 0
 
+    # Test Bootstrap
+    masked = dat.apply_mask(create_sphere(radius=10, coordinates=[0, 0, 0]))
+    n_samples = 3
+    b = masked.bootstrap('mean', n_samples=n_samples)
+    assert isinstance(b['Z'], Brain_Data)
+    b = masked.bootstrap('std', n_samples=n_samples)
+    assert isinstance(b['Z'], Brain_Data)
+    b = masked.bootstrap('predict', n_samples=n_samples, plot=False)
+    assert isinstance(b['Z'], Brain_Data)
+    b = masked.bootstrap('predict', n_samples=n_samples,
+                    plot=False, cv_dict={'type':'kfolds','n_folds':3})
+    assert isinstance(b['Z'], Brain_Data)
+    b = masked.bootstrap('predict', n_samples=n_samples,
+                    save_weights=True, plot=False)
+    assert len(b['samples'])==n_samples
+    
+    # Test decompose
+    n_components = 3
+    stats = dat.decompose(algorithm='pca', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    stats = dat.decompose(algorithm='ica', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    dat.data = dat.data + 2
+    dat.data[dat.data<0] = 0
+    stats = dat.decompose(algorithm='nnmf', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    stats = dat.decompose(algorithm='fa', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
 def test_brain_data_2mm(tmpdir):
     MNI_Template["resolution"] = '2mm'
     sim = Simulator()
@@ -468,6 +504,26 @@ def test_brain_data_2mm(tmpdir):
                     save_weights=True, plot=False)
     assert len(b['samples'])==n_samples
 
+    # Test decompose
+    n_components = 3
+    stats = dat.decompose(algorithm='pca', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    stats = dat.decompose(algorithm='ica', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    dat.data = dat.data + 2
+    dat.data[dat.data<0] = 0
+    stats = dat.decompose(algorithm='nnmf', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
+    stats = dat.decompose(algorithm='fa', n_components=n_components)
+    assert n_components == len(stats['components'])
+    assert stats['weights'].shape == (n_components,len(dat))
+
 def test_adjacency(tmpdir):
     n = 10
     sim = np.random.multivariate_normal([0,0,0,0],[[1, 0.8, 0.1, 0.4],

nltools/utils.py

@@ -17,6 +17,7 @@
             'all_same',
             'concatenate',
             '_bootstrap_apply_func',
+            'set_decomposition_algorithm'
             ]
 __author__ = ["Luke Chang"]
 __license__ = "MIT"
@@ -41,7 +42,7 @@ def get_anatomical():
     """
     return nib.load(os.path.join(get_resource_path(),'MNI152_T1_2mm.nii.gz'))
 
-def set_algorithm(algorithm, **kwargs):
+def set_algorithm(algorithm, *args, **kwargs):
     """ Setup the algorithm to use in subsequent prediction analyses.
 
     Args:
@@ -90,11 +91,11 @@ def load_class(import_string):
     if algorithm in algs_classify.keys():
         predictor_settings['prediction_type'] = 'classification'
         alg = load_class(algs_classify[algorithm])
-        predictor_settings['predictor'] = alg(**kwargs)
+        predictor_settings['predictor'] = alg(*args, **kwargs)
     elif algorithm in algs_predict:
         predictor_settings['prediction_type'] = 'prediction'
         alg = load_class(algs_predict[algorithm])
-        predictor_settings['predictor'] = alg(**kwargs)
+        predictor_settings['predictor'] = alg(*args, **kwargs)
     elif algorithm == 'lassopcr':
         predictor_settings['prediction_type'] = 'prediction'
         from sklearn.linear_model import Lasso
@@ -121,6 +122,44 @@ def load_class(import_string):
 
     return predictor_settings
 
+def set_decomposition_algorithm(algorithm, n_components=None, *args, **kwargs):
+    """ Setup the algorithm to use in subsequent decomposition analyses.
+
+    Args:
+        algorithm: The decomposition algorithm to use. Either a string or an
+                    (uninitialized) scikit-learn decomposition object.
+                    If string must be one of 'pca','nnmf', ica','fa'
+        kwargs: Additional keyword arguments to pass onto the scikit-learn
+                clustering object.
+
+    Returns:
+        predictor_settings: dictionary of settings for prediction
+
+    """
+
+    # NOTE: function currently located here instead of analysis.py to avoid circular imports
+
+    def load_class(import_string):
+        class_data = import_string.split(".")
+        module_path = '.'.join(class_data[:-1])
+        class_str = class_data[-1]
+        module = importlib.import_module(module_path)
+        return getattr(module, class_str)
+
+    algs = {
+        'pca': 'sklearn.decomposition.PCA',
+        'ica': 'sklearn.decomposition.FastICA',
+        'nnmf': 'sklearn.decomposition.NMF',
+        'fa': 'sklearn.decomposition.FactorAnalysis'
+        }
+
+    if algorithm in algs.keys():
+        alg = load_class(algs[algorithm])
+        alg = alg(n_components, *args, **kwargs)
+    else:
+        raise ValueError("""Invalid prediction/classification algorithm name.
+            Valid options are 'pca','ica', 'nnmf', 'fa'""")
+    return alg
 
 # The following are nipy source code implementations of the hemodynamic response function HRF
 # See the included nipy license file for use permission.

nltools/version.py

@@ -1,4 +1,4 @@
 """Specifies current version of nltools to be used by setup.py and __init__.py
 """
 
-__version__ = '0.3.4'
+__version__ = '0.3.5'

setup.py

@@ -1,30 +1,32 @@
-# from nltools.version import __version__
 from setuptools import setup, find_packages
 
-__version__ = '0.3.4'
+version = {}
+with open("nltools/version.py") as f:
+    exec(f.read(), version)
+
+with open('requirements.txt') as f:
+    requirements = f.read().splitlines()
 
 extra_setuptools_args = dict(
     tests_require=['pytest']
 )
 
-with open('requirements.txt') as f:
-    requirements = f.read().splitlines()
-
 setup(
-    name='nltools',
-    version=__version__,
-    author='Luke Chang',
-    author_email='[email protected]',
-    url='http://neurolearn.readthedocs.org/en/latest/',
-    install_requires=requirements,
-    extras_require={
+    name = 'nltools',
+    version = version['__version__'],
+    author = 'Luke Chang',
+    author_email = '[email protected]',
+    url = 'http://neurolearn.readthedocs.org/en/latest/',
+    install_requires = requirements,
+    extras_require = {
     'ibrainViewer':['ipywidgets>=5.2.2']
     },
-    packages=find_packages(exclude=['nltools/tests']),
-    package_data={'nltools': ['resources/*']},
-    license='LICENSE.txt',
-    description='A Python package to analyze neuroimaging data',
-    long_description='nltools is a collection of python tools to perform '
+    packages = find_packages(exclude=['nltools/tests']),
+    package_data = {'nltools': ['resources/*']},
+    include_package_data = True,
+    license = 'LICENSE.txt',
+    description = 'A Python package to analyze neuroimaging data',
+    long_description = 'nltools is a collection of python tools to perform '
                      'preprocessing, univariate GLMs, and predictive '
                      'multivariate modeling of neuroimaging data. It is the '
                      'analysis engine powering www.neuro-learn.org.',