From c667a0f6518e1a6c1559a93c6b400e14d524e495 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 14:40:38 -0800 Subject: [PATCH 01/32] implemented RSRM --- brainiak/funcalign/rsrm.py | 515 +++++++++++++++++++++++++++++++++++++ 1 file changed, 515 insertions(+) create mode 100644 brainiak/funcalign/rsrm.py diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py new file mode 100644 index 000000000..ff31443e6 --- /dev/null +++ b/brainiak/funcalign/rsrm.py @@ -0,0 +1,515 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Robust Shared Response Model (RSRM) + +The implementation is based on the following publications: + +.. [Turek2017] "Capturing Shared and Individual Information in fMRI Data", + J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke + under review, 2017. +""" + +# Authors: Javier Turek (Intel Labs), 2017 + +import logging + +import numpy as np +from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin +from sklearn.utils import assert_all_finite +from sklearn.utils.validation import NotFittedError + +__all__ = [ + "RSRM" +] + +logger = logging.getLogger(__name__) + +class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): + """Robust Shared Response Model (RSRM) + + Given multi-subject data, factorize it as a shared response R among all + subjects, an orthogonal transform W per subject, and an outlying (individual) + sparse component S per subject: + + .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N + + Parameters + ---------- + + n_iter : int, default: 10 + Number of iterations to run the algorithm. + + features : int, default: 50 + Number of features to compute. + + gamma : float, default: 1.0 + Regularization parameter for the sparseness of the individual components. + + rand_seed : int, default: 0 + Seed for initializing the random number generator. + + + Attributes + ---------- + + w_ : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) for each subject. + + r_ : array, shape=[features, samples] + The shared response. + + s_ : list of array, element i has shape=[voxels_i, samples] + The individual components with outlying data for each subject. + + random_state_: `RandomState` + Random number generator initialized using rand_seed + + Note + ---- + + The number of voxels may be different between subjects. However, the + number of samples for the alignment data must be the same across + subjects. + + The Robust Shared Response Model is approximated using the + Block-Coordinate Descent (BCD) algorithm proposed in [Turek2017]_. + + This is a single node version. + """ + + def __init__(self, n_iter=10, features=50, gamma=1.0, rand_seed=0): + self.n_iter = n_iter + self.features = features + self.lam = gamma + self.rand_seed = rand_seed + + def fit(self, X, y=None): + """Compute the Robust Shared Response Model + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data of one subject. + + y : not used + """ + logger.info('Starting RSRM') + + # Check that the regularizer value is positive + if 0.0 >= self.lam: + raise ValueError("Gamma parameter should be positive.") + + # Check the number of subjects + if len(X) <= 1: + raise ValueError("There are not enough subjects in the input " + "data to train the model.") + + # Check for input data sizes + if X[0].shape[1] < self.features: + raise ValueError( + "There are not enough samples to train the model with " + "{0:d} features.".format(self.features)) + + # Check if all subjects have same number of TRs for alignment + number_trs = X[0].shape[1] + number_subjects = len(X) + for subject in range(number_subjects): + assert_all_finite(X[subject]) + if X[subject].shape[1] != number_trs: + raise ValueError("Different number of alignment samples " + "between subjects.") + + # Create a new random state + self.random_state_ = np.random.RandomState(self.rand_seed) + + # Run RSRM + self.w_, self.r_, self.s_ = self._rsrm(X) + + return self + + def transform(self, X, y=None): + """Use the model to transform new data to Shared Response space + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples_i] + Each element in the list contains the fMRI data of one subject. + + y : not used + + Returns + ------- + + r : list of 2D arrays, element i has shape=[features_i, samples_i] + Shared responses from input data (X) + + s : list of 2D arrays, element i has shape=[voxels_i, samples_i] + Outlying data obtained from fitting model to input data (X) + """ + # Check if the model exist + if hasattr(self, 'w_') is False: + raise NotFittedError("The model fit has not been run yet.") + + # Check the number of subjects + if len(X) != len(self.w_): + raise ValueError("The number of subjects does not match the one" + " in the model.") + + r = [None] * len(X) + s = [None] * len(X) + for subject in range(len(X)): + if X[subject] is not None: + r[subject], s[subject] = self._transform_new_data(X[subject], subject) + + return r, s + + def _transform_new_data(self, X, subject): + """Transform new data for a subjects by projecting to the shared subspace and + computing the outlying individual information. + + Parameters + ---------- + + X : array, shape=[voxels, samples] + The fMRI data of the subject. + + subject : int + The subject id. + + Returns + ------- + + R : array, shape=[features, samples] + Shared response from input data (X) + + S : array, shape=[voxels, samples] + Outlying data obtained from fitting model to input data (X) + """ + S = np.zeros_like(X) + R = None + for iter in range(self.n_iter): + R = self.w_[subject].T.dot(X - S) + S = self._shrink(X - self.w_[subject].dot(R), self.lam) + return R, S + + def _rsrm(self, X): + """Block-Coordinate Descent algorithm for fitting RSRM. + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data for alignment of + one subject. + + Returns + ------- + + W : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) :math:`W_i` for each subject. + + R : array, shape=[features, samples] + The shared response. + + S : list of array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject. + """ + subjs = len(X) + voxels = [X[i].shape[0] for i in range(subjs)] + TRs = X[0].shape[1] + features = self.features + + # Initialization + W = self._init_transforms(subjs, voxels, features, self.random_state_) + S = self._init_outliers(subjs, voxels, TRs) + R = self._update_shared_response(X, S, W, subjs, features, TRs) + + if logger.isEnabledFor(logging.INFO): + objective = self._objective_function(X, W, R, S, subjs, self.lam) + logger.info('Objective function %f' % objective) + + # Main loop + for iter in range(self.n_iter): + W = self._update_transforms(X, S, R, subjs) + S = self._update_outliers(X, W, R, subjs, self.lam) + R = self._update_shared_response(X, S, W, subjs, features, TRs) + # Print objective function every iteration + if logger.isEnabledFor(logging.INFO): + objective = self._objective_function(X, W, R, S, subjs, self.lam) + logger.info('Objective function %f' % objective) + + return W, R, S + + def _init_transforms(self, subjs, voxels, features, random_state): + """Initialize the mappings (Wi) with random orthogonal matrices. + Parameters + ---------- + + subjs : int + The number of subjects. + + voxels : list of int + A list with the number of voxels per subject. + + features : int + The number of features in the model. + + random_state : `RandomState` + A random state to draw the mappings. + + Returns + ------- + + W : list of array, element i has shape=[voxels_i, features] + The initialized orthogonal transforms (mappings) :math:`W_i` for each + subject. + + Note + ---- + Not thread safe. + """ + + # Init the Random seed generator + np.random.seed(self.rand_seed) + + # Draw a random W for each subject + W = [random_state.random_sample((voxels[i], features)) for i in range(subjs)] + for i in range(subjs): + W[i], _ = np.linalg.qr(W[i]) + + return W + + @staticmethod + def _objective_function(X, W, R, S, subjs, gamma): + """Evaluate the objective function. + + .. math:: \sum_{i=1}^{N} 1/2 \| X_i - W_i R - S_i \|_F^2 + .. math:: + /\gamma * \|S_i\|_1 + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data for alignment of + one subject. + + W : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) :math:`W_i` for each subject. + + R : array, shape=[features, samples] + The shared response. + + S : list of array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject. + + subjs : int + The number of subjects. + + gamma : float, default: 1.0 + Regularization parameter for the sparseness of the individual components. + + Returns + ------- + + func : float + The RSRM objective function evaluated on the parameters to this function. + """ + func = .0 + for i in range(subjs): + func += 0.5 * np.sum((X[i] - W[i].dot(R) - S[i])**2) \ + + gamma * np.sum(np.abs(S[i])) + return func + + @staticmethod + def _update_outliers(X, W, R, subjs, gamma): + """Updates the outlier terms `S_i`. + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data for alignment of + one subject. + + W : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) :math:`W_i` for each subject. + + R : array, shape=[features, samples] + The shared response. + + subjs : int + The number of subjects. + + gamma : float, default: 1.0 + Regularization parameter for the sparseness of the individual components. + + Returns + ------- + + S : list of array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject. + """ + S = [] + for i in range(subjs): + S.append(RSRM._shrink(X[i] - W[i].dot(R), gamma)) + return S + + @staticmethod + def _init_outliers(subjs, voxels, TRs): + """Initializes the outlier terms `S_i` to empty (all zeros). + + Parameters + ---------- + + subjs : int + The number of subjects. + + voxels : list of int + A list with the number of voxels per subject. + + TRs : int + The number of samples in the data. + + Returns + ------- + + S : list of 2D array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject initialized to zero. + """ + return [np.zeros((voxels[i], TRs)) for i in range(subjs)] + + @staticmethod + def _update_shared_response(X, S, W, subjs, features, TRs): + """Update the shared response `R`. + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data for alignment of + one subject. + + S : list of array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject. + + W : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) :math:`W_i` for each subject. + + subjs : int + The number of subjects. + + features : int + The number of features in the model. + + TRs : int + The number of samples in the data. + + Returns + ------- + + R : array, shape=[features, samples] + The updated shared response. + """ + R = np.zeros((features, TRs)) + for i in range(subjs): + R += W[i].T.dot(X[i]) + R -= W[i].T.dot(S[i]) + R /= subjs + return R + + @staticmethod + def _update_transform_subject(X, S, R): + """Updates the mappings `W_i` for one subject. + + Parameters + ---------- + + X : array, shape=[voxels, samples] + The fMRI data for aligning the subject. + + S : array, shape=[voxels, samples] + The individual component :math:`S_i` for the subject. + + R : array, shape=[features, samples] + The shared response. + + Returns + ------- + + W : array, shape=[voxels, features] + The orthogonal transform (mapping) :math:`W_i` for the subject. + """ + A = X.dot(R.T) + A -= S.dot(R.T) + U, _, V = np.linalg.svd(A, full_matrices=False) + return U.dot(V) + + @staticmethod + def _update_transforms(X, S, R, subjs): + """Updates the mappings `W_i` for each subject. + + Parameters + ---------- + + X : list of 2D arrays, element i has shape=[voxels_i, samples] + Each element in the list contains the fMRI data for alignment of + one subject. + + S : list of array, element i has shape=[voxels_i, samples] + The individual component :math:`S_i` for each subject. + + R : array, shape=[features, samples] + The shared response. + + subjs : int + The number of subjects. + + Returns + ------- + + W : list of array, element i has shape=[voxels_i, features] + The orthogonal transforms (mappings) :math:`W_i` for each subject. + """ + W = [] + for i in range(subjs): + W.append(RSRM._update_transform_subject(X[i], S[i], R)) + return W + + @staticmethod + def _shrink(v, gamma): + """Soft-shrinkage of an array with parameter gamma. + + Parameters + ---------- + + v : array + Array containing the values to be applied to the shrinkage operator + + gamma : float + Shrinkage parameter. + + Returns + ------- + + v : array + The same input array after the shrinkage operator was applied. + """ + pos = v > gamma + neg = v < -gamma + v[pos] -= gamma + v[neg] += gamma + v[np.logical_and(~pos, ~neg)] = .0 + return v From 5f9ebc6bb3dfe5fe5c882301a7f63c35c0c2c87e Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 14:48:14 -0800 Subject: [PATCH 02/32] Transform new subject and naming convention change --- brainiak/funcalign/rsrm.py | 40 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index ff31443e6..1571a38b4 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -200,11 +200,47 @@ def _transform_new_data(self, X, subject): """ S = np.zeros_like(X) R = None - for iter in range(self.n_iter): + for i in range(self.n_iter): R = self.w_[subject].T.dot(X - S) S = self._shrink(X - self.w_[subject].dot(R), self.lam) return R, S + def transform_subject(self, X, y=None): + """Transform a new subject using the existing model + + Parameters + ---------- + + X : 2D array, shape=[voxels, samples] + The fMRI data of the new subject. + + y : not used + + Returns + ------- + + w : 2D array, shape=[voxels, features] + Orthogonal mapping `W_{new}` for new subject + + s : 2D array, shape=[voxels, samples] + Individual term `S_{new}` for new subject + """ + # Check if the model exist + if hasattr(self, 'w_') is False: + raise NotFittedError("The model fit has not been run yet.") + + # Check the number of TRs in the subject + if X.shape[1] != self.r_.shape[1]: + raise ValueError("The number of samples(TRs) does not match the one" + " in the model.") + + s = np.zeros_like(X) + for i in range(self.n_iter): + w = self._update_transform_subject(X, s, self.r_) + s = self._shrink(X - w.dot(self.r_), self.lam) + + return w, s + def _rsrm(self, X): """Block-Coordinate Descent algorithm for fitting RSRM. @@ -242,7 +278,7 @@ def _rsrm(self, X): logger.info('Objective function %f' % objective) # Main loop - for iter in range(self.n_iter): + for i in range(self.n_iter): W = self._update_transforms(X, S, R, subjs) S = self._update_outliers(X, W, R, subjs, self.lam) R = self._update_shared_response(X, S, W, subjs, features, TRs) From 44dea7ed9e89388e43b932ab6857049c8048509a Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:08:55 -0800 Subject: [PATCH 03/32] added test for rsrm --- tests/funcalign/test_rsrm.py | 124 +++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 tests/funcalign/test_rsrm.py diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py new file mode 100644 index 000000000..b1150c775 --- /dev/null +++ b/tests/funcalign/test_rsrm.py @@ -0,0 +1,124 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from sklearn.exceptions import NotFittedError +import pytest + + +def test_can_instantiate(): + import brainiak.funcalign.rsrm + s = brainiak.funcalign.rsrm.RSRM() + assert s, "Invalid RSRM instance!" + + import numpy as np + np.random.seed(0) + + voxels = 100 + samples = 500 + subjects = 2 + features = 3 + + s = brainiak.funcalign.rsrm.RSRM(n_iter=5, features=features) + assert s, "Invalid RSRM instance!" + + # Create a Shared response S with K = 3 + theta = np.linspace(-4 * np.pi, 4 * np.pi, samples) + z = np.linspace(-2, 2, samples) + r = z**2 + 1 + x = r * np.sin(theta) + y = r * np.cos(theta) + + S = np.vstack((x, y, z)) + + X = [] + W = [] + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + + # Check that transform does NOT run before fitting the model + with pytest.raises(NotFittedError): + s.transform(X) + print("Test: transforming before fitting the model") + + # Check that it does NOT run with 1 subject + with pytest.raises(ValueError): + s.fit(X) + print("Test: running RSRM with 1 subject") + + for subject in range(1, subjects): + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + + # Check that runs with 2 subject + s.fit(X) + from pathlib import Path + sr_v0_4 = np.load(Path(__file__).parent / "sr_v0_4.npz")['sr'] + assert(np.allclose(sr_v0_4, s.s_)) + + assert len(s.w_) == subjects, ( + "Invalid computation of RSRM! (wrong # subjects in W)") + for subject in range(subjects): + assert s.s_[subject].shape[0] == voxels, ( + "Invalid computation of RSRM! (wrong # voxels in S)") + assert s.s_[subject].shape[1] == samples, ( + "Invalid computation of RSRM! (wrong # samples in S)") + assert s.w_[subject].shape[0] == voxels, ( + "Invalid computation of RSRM! (wrong # voxels in W)") + assert s.w_[subject].shape[1] == features, ( + "Invalid computation of RSRM! (wrong # features in W)") + ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) + - np.eye(s.w_[subject].shape[1]), + 'fro') + assert ortho < 1e-7, "A Wi mapping is not orthonormal in RSRM." + difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.r_), + 'fro') + datanorm = np.linalg.norm(X[subject], 'fro') + assert difference/datanorm < 1.0, "Model seems incorrectly computed." + assert s.r_.shape[0] == features, ( + "Invalid computation of RSRM! (wrong # features in R)") + assert s.r_.shape[1] == samples, ( + "Invalid computation of RSRM! (wrong # samples in R)") + + # Check that it does run to compute the shared response after the model + # computation + new_r, _ = s.transform(X) + + assert len(new_r) == subjects, ( + "Invalid computation of RSRM! (wrong # subjects after transform)") + for subject in range(subjects): + assert new_r[subject].shape[0] == features, ( + "Invalid computation of RSRM! (wrong # features after transform)") + assert new_r[subject].shape[1] == samples, ( + "Invalid computation of RSRM! (wrong # samples after transform)") + + # Check that it does NOT run with non-matching number of subjects + with pytest.raises(ValueError): + s.transform(X[1]) + print("Test: transforming with non-matching number of subjects") + + # Check that it does not run without enough samples (TRs). + with pytest.raises(ValueError): + s.set_params(features=(samples+1)) + s.fit(X) + print("Test: not enough samples") + + # Check that it does not run with different number of samples (TRs) + S2 = S[:, :-2] + X.append(Q.dot(S2)) + with pytest.raises(ValueError): + s.fit(X) + print("Test: different number of samples per subject") + + From a82ad1b2e7c4678ce204b18e8874b9967b2baf59 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:31:46 -0800 Subject: [PATCH 04/32] corrected some PEP issues --- brainiak/funcalign/rsrm.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 1571a38b4..99640a1d2 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -35,12 +35,13 @@ logger = logging.getLogger(__name__) + class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): """Robust Shared Response Model (RSRM) Given multi-subject data, factorize it as a shared response R among all - subjects, an orthogonal transform W per subject, and an outlying (individual) - sparse component S per subject: + subjects, an orthogonal transform W per subject, and an outlying + (individual) sparse component S per subject: .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N @@ -54,7 +55,8 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Number of features to compute. gamma : float, default: 1.0 - Regularization parameter for the sparseness of the individual components. + Regularization parameter for the sparseness of the individual + components. rand_seed : int, default: 0 Seed for initializing the random number generator. @@ -62,7 +64,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Attributes ---------- - + w_ : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) for each subject. @@ -77,7 +79,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Note ---- - + The number of voxels may be different between subjects. However, the number of samples for the alignment data must be the same across subjects. From decb785e3f9d1126ec8e31d7aae1ca0da4ade787 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:39:57 -0800 Subject: [PATCH 05/32] corrected more PEP issues --- brainiak/funcalign/rsrm.py | 89 ++++++++++++++++++++------------------ 1 file changed, 48 insertions(+), 41 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 99640a1d2..7945091e3 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -40,7 +40,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): """Robust Shared Response Model (RSRM) Given multi-subject data, factorize it as a shared response R among all - subjects, an orthogonal transform W per subject, and an outlying + subjects, an orthogonal transform W per subject, and an outlying (individual) sparse component S per subject: .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N @@ -55,7 +55,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Number of features to compute. gamma : float, default: 1.0 - Regularization parameter for the sparseness of the individual + Regularization parameter for the sparseness of the individual components. rand_seed : int, default: 0 @@ -64,7 +64,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Attributes ---------- - + w_ : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) for each subject. @@ -79,10 +79,10 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): Note ---- - + The number of voxels may be different between subjects. However, the number of samples for the alignment data must be the same across - subjects. + subjects. The Robust Shared Response Model is approximated using the Block-Coordinate Descent (BCD) algorithm proposed in [Turek2017]_. @@ -149,15 +149,15 @@ def transform(self, X, y=None): X : list of 2D arrays, element i has shape=[voxels_i, samples_i] Each element in the list contains the fMRI data of one subject. - - y : not used - + + y : not used + Returns ------- r : list of 2D arrays, element i has shape=[features_i, samples_i] Shared responses from input data (X) - + s : list of 2D arrays, element i has shape=[voxels_i, samples_i] Outlying data obtained from fitting model to input data (X) """ @@ -174,7 +174,8 @@ def transform(self, X, y=None): s = [None] * len(X) for subject in range(len(X)): if X[subject] is not None: - r[subject], s[subject] = self._transform_new_data(X[subject], subject) + r[subject], s[subject] = self._transform_new_data(X[subject], + subject) return r, s @@ -187,7 +188,7 @@ def _transform_new_data(self, X, subject): X : array, shape=[voxels, samples] The fMRI data of the subject. - + subject : int The subject id. @@ -222,7 +223,7 @@ def transform_subject(self, X, y=None): ------- w : 2D array, shape=[voxels, features] - Orthogonal mapping `W_{new}` for new subject + Orthogonal mapping `W_{new}` for new subject s : 2D array, shape=[voxels, samples] Individual term `S_{new}` for new subject @@ -233,8 +234,8 @@ def transform_subject(self, X, y=None): # Check the number of TRs in the subject if X.shape[1] != self.r_.shape[1]: - raise ValueError("The number of samples(TRs) does not match the one" - " in the model.") + raise ValueError("The number of samples(TRs) does not match the" + "one in the model.") s = np.zeros_like(X) for i in range(self.n_iter): @@ -251,7 +252,7 @@ def _rsrm(self, X): X : list of 2D arrays, element i has shape=[voxels_i, samples] Each element in the list contains the fMRI data for alignment of - one subject. + one subject. Returns ------- @@ -261,7 +262,7 @@ def _rsrm(self, X): R : array, shape=[features, samples] The shared response. - + S : list of array, element i has shape=[voxels_i, samples] The individual component :math:`S_i` for each subject. """ @@ -286,7 +287,8 @@ def _rsrm(self, X): R = self._update_shared_response(X, S, W, subjs, features, TRs) # Print objective function every iteration if logger.isEnabledFor(logging.INFO): - objective = self._objective_function(X, W, R, S, subjs, self.lam) + objective = self._objective_function(X, W, R, S, subjs, + self.lam) logger.info('Objective function %f' % objective) return W, R, S @@ -304,7 +306,7 @@ def _init_transforms(self, subjs, voxels, features, random_state): features : int The number of features in the model. - + random_state : `RandomState` A random state to draw the mappings. @@ -312,8 +314,8 @@ def _init_transforms(self, subjs, voxels, features, random_state): ------- W : list of array, element i has shape=[voxels_i, features] - The initialized orthogonal transforms (mappings) :math:`W_i` for each - subject. + The initialized orthogonal transforms (mappings) :math:`W_i` for + each subject. Note ---- @@ -324,7 +326,8 @@ def _init_transforms(self, subjs, voxels, features, random_state): np.random.seed(self.rand_seed) # Draw a random W for each subject - W = [random_state.random_sample((voxels[i], features)) for i in range(subjs)] + W = [random_state.random_sample((voxels[i], features)) + for i in range(subjs)] for i in range(subjs): W[i], _ = np.linalg.qr(W[i]) @@ -333,37 +336,39 @@ def _init_transforms(self, subjs, voxels, features, random_state): @staticmethod def _objective_function(X, W, R, S, subjs, gamma): """Evaluate the objective function. - + .. math:: \sum_{i=1}^{N} 1/2 \| X_i - W_i R - S_i \|_F^2 .. math:: + /\gamma * \|S_i\|_1 - + Parameters ---------- - + X : list of 2D arrays, element i has shape=[voxels_i, samples] Each element in the list contains the fMRI data for alignment of - one subject. + one subject. W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. R : array, shape=[features, samples] The shared response. - + S : list of array, element i has shape=[voxels_i, samples] The individual component :math:`S_i` for each subject. - + subjs : int The number of subjects. - + gamma : float, default: 1.0 - Regularization parameter for the sparseness of the individual components. - + Regularization parameter for the sparseness of the individual + components. + Returns ------- - + func : float - The RSRM objective function evaluated on the parameters to this function. + The RSRM objective function evaluated on the parameters to this + function. """ func = .0 for i in range(subjs): @@ -380,7 +385,7 @@ def _update_outliers(X, W, R, subjs, gamma): X : list of 2D arrays, element i has shape=[voxels_i, samples] Each element in the list contains the fMRI data for alignment of - one subject. + one subject. W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. @@ -392,7 +397,8 @@ def _update_outliers(X, W, R, subjs, gamma): The number of subjects. gamma : float, default: 1.0 - Regularization parameter for the sparseness of the individual components. + Regularization parameter for the sparseness of the individual + components. Returns ------- @@ -425,7 +431,8 @@ def _init_outliers(subjs, voxels, TRs): ------- S : list of 2D array, element i has shape=[voxels_i, samples] - The individual component :math:`S_i` for each subject initialized to zero. + The individual component :math:`S_i` for each subject initialized + to zero. """ return [np.zeros((voxels[i], TRs)) for i in range(subjs)] @@ -438,11 +445,11 @@ def _update_shared_response(X, S, W, subjs, features, TRs): X : list of 2D arrays, element i has shape=[voxels_i, samples] Each element in the list contains the fMRI data for alignment of - one subject. + one subject. S : list of array, element i has shape=[voxels_i, samples] The individual component :math:`S_i` for each subject. - + W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. @@ -451,7 +458,7 @@ def _update_shared_response(X, S, W, subjs, features, TRs): features : int The number of features in the model. - + TRs : int The number of samples in the data. @@ -476,7 +483,7 @@ def _update_transform_subject(X, S, R): ---------- X : array, shape=[voxels, samples] - The fMRI data for aligning the subject. + The fMRI data for aligning the subject. S : array, shape=[voxels, samples] The individual component :math:`S_i` for the subject. @@ -504,7 +511,7 @@ def _update_transforms(X, S, R, subjs): X : list of 2D arrays, element i has shape=[voxels_i, samples] Each element in the list contains the fMRI data for alignment of - one subject. + one subject.ß S : list of array, element i has shape=[voxels_i, samples] The individual component :math:`S_i` for each subject. @@ -535,7 +542,7 @@ def _shrink(v, gamma): v : array Array containing the values to be applied to the shrinkage operator - + gamma : float Shrinkage parameter. From 9ce3c748047b3bcc466750c9c4d7de53b11b10bf Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:41:59 -0800 Subject: [PATCH 06/32] corrected last PEP issues - luckily --- brainiak/funcalign/rsrm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 7945091e3..a0597aa2b 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -314,7 +314,7 @@ def _init_transforms(self, subjs, voxels, features, random_state): ------- W : list of array, element i has shape=[voxels_i, features] - The initialized orthogonal transforms (mappings) :math:`W_i` for + The initialized orthogonal transforms (mappings) :math:`W_i` for each subject. Note From 2d1df926bf34aa4f16a5342107069f7f436da742 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:44:13 -0800 Subject: [PATCH 07/32] corrected PEP issues for test --- tests/funcalign/test_rsrm.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py index b1150c775..2746528de 100644 --- a/tests/funcalign/test_rsrm.py +++ b/tests/funcalign/test_rsrm.py @@ -120,5 +120,3 @@ def test_can_instantiate(): with pytest.raises(ValueError): s.fit(X) print("Test: different number of samples per subject") - - From 9677ce87264439f022260d25a62757715d980011 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 15:54:23 -0800 Subject: [PATCH 08/32] corrected test for wrong result --- tests/funcalign/test_rsrm.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py index 2746528de..25a16c452 100644 --- a/tests/funcalign/test_rsrm.py +++ b/tests/funcalign/test_rsrm.py @@ -63,9 +63,6 @@ def test_can_instantiate(): # Check that runs with 2 subject s.fit(X) - from pathlib import Path - sr_v0_4 = np.load(Path(__file__).parent / "sr_v0_4.npz")['sr'] - assert(np.allclose(sr_v0_4, s.s_)) assert len(s.w_) == subjects, ( "Invalid computation of RSRM! (wrong # subjects in W)") From a3040c530a5be55d56d291723ae466edf676a3b6 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 16:42:16 -0800 Subject: [PATCH 09/32] additional funciton tested --- tests/funcalign/test_rsrm.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py index 25a16c452..8e010d164 100644 --- a/tests/funcalign/test_rsrm.py +++ b/tests/funcalign/test_rsrm.py @@ -100,6 +100,17 @@ def test_can_instantiate(): assert new_r[subject].shape[1] == samples, ( "Invalid computation of RSRM! (wrong # samples after transform)") + # Check that it does run to compute a new subject + new_w, new_s = s.transform_subject(X) + assert new_w.shape[1] == features, ( + "Invalid computation of RSRM! (wrong # features for new subject)") + assert new_s.shape[1] == samples, ( + "Invalid computation of RSRM! (wrong # samples for new subject)") + assert new_s.shape[0] == voxels, ( + "Invalid computation of RSRM! (wrong # voxels for new subject)") + assert new_w.shape[0] == voxels, ( + "Invalid computation of RSRM! (wrong # voxels for new subject)") + # Check that it does NOT run with non-matching number of subjects with pytest.raises(ValueError): s.transform(X[1]) From 23adcd2ae6b93e69a57b23f2a4142627f4545926 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 16:45:45 -0800 Subject: [PATCH 10/32] corrected test --- examples/funcalign/srm_image_prediction_example.ipynb | 6 +++--- tests/funcalign/test_rsrm.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/funcalign/srm_image_prediction_example.ipynb b/examples/funcalign/srm_image_prediction_example.ipynb index 52610f3d9..378a50173 100644 --- a/examples/funcalign/srm_image_prediction_example.ipynb +++ b/examples/funcalign/srm_image_prediction_example.ipynb @@ -404,7 +404,7 @@ " cm[subject] = confusion_matrix(test_labels, predicted_labels)\n", "\n", " # Normalize the confusion matrix\n", - " cm[subject] = cm[subject].astype('float') / cm[subject].sum(axis=1)[:, np.newaxis]\n" + " cm[subject] = cm[subject].astype('float') / cm[subject].sum(axis=1)[:, np.newaxis]" ] }, { @@ -436,7 +436,7 @@ "language_info": { "codemirror_mode": { "name": "ipython", - "version": 3 + "version": 3.0 }, "file_extension": ".py", "mimetype": "text/x-python", @@ -448,4 +448,4 @@ }, "nbformat": 4, "nbformat_minor": 0 -} +} \ No newline at end of file diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py index 8e010d164..b723a58c7 100644 --- a/tests/funcalign/test_rsrm.py +++ b/tests/funcalign/test_rsrm.py @@ -101,7 +101,7 @@ def test_can_instantiate(): "Invalid computation of RSRM! (wrong # samples after transform)") # Check that it does run to compute a new subject - new_w, new_s = s.transform_subject(X) + new_w, new_s = s.transform_subject(X[0]) assert new_w.shape[1] == features, ( "Invalid computation of RSRM! (wrong # features for new subject)") assert new_s.shape[1] == samples, ( From 88a6c9159be31c31a55141d3047c800ced49de0f Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 20:45:25 -0800 Subject: [PATCH 11/32] example with rsrm --- .../rsrm_synthetic_reconstruction.ipynb | 185 ++++++++++++++++++ 1 file changed, 185 insertions(+) create mode 100644 examples/funcalign/rsrm_synthetic_reconstruction.ipynb diff --git a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb new file mode 100644 index 000000000..728da32da --- /dev/null +++ b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb @@ -0,0 +1,185 @@ +{ + "cells": [ + { + "cell_type": "heading", + "metadata": {}, + "level": 1, + "source": [ + "This example replicates the synthetic experiment in the paper\n", + "\"Capturing Shared and Individual Information in fMRI Data\", J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Import some libraries that we will need" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import scipy.io\n", + "from scipy.stats import stats\n", + "import numpy as np\n", + "import brainiak.funcalign.srm.SRM as SRM\n", + "import brainiak.funcalign.rsrm.RSRM as RSRM\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Experiment setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "voxels = 100\n", + "samples = 200\n", + "subjects = 10\n", + "features = 3\n", + "snr = 20 # in dB" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we create some synthetic data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a Shared response R with K = 3\n", + "theta = np.linspace(-4 * np.pi, 4 * np.pi, samples)\n", + "z = np.linspace(-2, 2, samples)\n", + "r = z ** 2 + 1\n", + "x = r * np.sin(theta)\n", + "y = r * np.cos(theta)\n", + "curve = np.vstack((x, y, z))\n", + "print('Print curve max, min values:', np.max(curve), np.min(curve))\n", + "\n", + "# Create the subjects' data\n", + "data = [None] * subjects\n", + "W = [None] * subjects\n", + "noise_level = 0.0\n", + "for s in range(subjects):\n", + " R = curve\n", + " W[s], _ = np.linalg.qr(np.random.randn(voxels, 3))\n", + " data[s] = W[s].dot(R)\n", + " noise_level += np.sum(np.abs(data[s])**2)\n", + "\n", + "# Compute noise_sigma from desired SNR\n", + "noise_level = noise_level / (10 ** (snr / 10))\n", + "noise_level = np.sqrt(noise_level / subjects / voxels / samples)\n", + "print(noise_level)\n", + "\n", + "for s in range(subjects):\n", + " n = noise_level * np.random.randn(voxels, samples)\n", + " data[s] += n\n", + " print(20 * np.log10(np.sum(np.abs(data[s])) / np.sum(np.abs(n)) ))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we fit the algorithms, SRM and RSRM, to the synthetic data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "srm = SRM(features=3, n_iter=20)\n", + "srm.fit(data)\n", + "\n", + "\n", + "rsrm = RSRM(features=3, gamma=0.35, n_iter=20)\n", + "rsrm.fit(data)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following function finds the orthogonal transform to align the shared response to the original curve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def find_orthogonal_transform(shared_response):\n", + " u,_,vt = np.linalg.svd(shared_response.dot(curve.T))\n", + " q = u.dot(vt)\n", + " aligned_curve = q.T.dot(shared_response)\n", + " return aligned_curve" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca(projection='3d')\n", + "ax.plot(curve[0, :], curve[1, :], curve[2, :], '-g', label='original', lineWidth=5)\n", + "\n", + "proj = find_orthogonal_transform(srm.s_)\n", + "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-b', label='SRM', lineWidth=3)\n", + "proj = find_orthogonal_transform(rsrm.r_)\n", + "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2.0 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file From 7da4e7a42b00e94bef393e55f0a7da7a79324129 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 20:57:22 -0800 Subject: [PATCH 12/32] corrected notebook --- examples/funcalign/rsrm_synthetic_reconstruction.ipynb | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb index 728da32da..f55991a27 100644 --- a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb +++ b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb @@ -1,9 +1,8 @@ { "cells": [ { - "cell_type": "heading", + "cell_type": "markdown", "metadata": {}, - "level": 1, "source": [ "This example replicates the synthetic experiment in the paper\n", "\"Capturing Shared and Individual Information in fMRI Data\", J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke" @@ -156,8 +155,7 @@ "proj = find_orthogonal_transform(srm.s_)\n", "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-b', label='SRM', lineWidth=3)\n", "proj = find_orthogonal_transform(rsrm.r_)\n", - "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)\n", - "\n" + "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)" ] } ], From 2973a67f718053e4c7ad64beda99d660ec402353 Mon Sep 17 00:00:00 2001 From: TuKo Date: Wed, 22 Nov 2017 21:34:50 -0800 Subject: [PATCH 13/32] corrected a few errors in example --- .../rsrm_synthetic_reconstruction.ipynb | 54 +++++++++++-------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb index f55991a27..b82ab4521 100644 --- a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb +++ b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb @@ -4,8 +4,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This example replicates the synthetic experiment in the paper\n", - "\"Capturing Shared and Individual Information in fMRI Data\", J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke" + "## This example replicates the synthetic experiment in the paper\n", + "## \"Capturing Shared and Individual Information in fMRI Data\", J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke" ] }, { @@ -25,8 +25,8 @@ "import scipy.io\n", "from scipy.stats import stats\n", "import numpy as np\n", - "import brainiak.funcalign.srm.SRM as SRM\n", - "import brainiak.funcalign.rsrm.RSRM as RSRM\n", + "import brainiak.funcalign.srm as srm\n", + "import brainiak.funcalign.rsrm as rsrm\n", "from mpl_toolkits.mplot3d import Axes3D\n", "import matplotlib.pyplot as plt" ] @@ -48,7 +48,9 @@ "samples = 200\n", "subjects = 10\n", "features = 3\n", - "snr = 20 # in dB" + "snr = 20 # in dB\n", + "amplitude = 8.0 # approximate amplitude used for the values in Si\n", + "p = 0.05 # probability of being non-zero in Si" ] }, { @@ -71,7 +73,7 @@ "x = r * np.sin(theta)\n", "y = r * np.cos(theta)\n", "curve = np.vstack((x, y, z))\n", - "print('Print curve max, min values:', np.max(curve), np.min(curve))\n", + "print('Curve max, min values:', np.max(curve), np.min(curve))\n", "\n", "# Create the subjects' data\n", "data = [None] * subjects\n", @@ -86,12 +88,15 @@ "# Compute noise_sigma from desired SNR\n", "noise_level = noise_level / (10 ** (snr / 10))\n", "noise_level = np.sqrt(noise_level / subjects / voxels / samples)\n", - "print(noise_level)\n", "\n", "for s in range(subjects):\n", " n = noise_level * np.random.randn(voxels, samples)\n", " data[s] += n\n", - " print(20 * np.log10(np.sum(np.abs(data[s])) / np.sum(np.abs(n)) ))" + " \n", + "S = [None] * subjects\n", + "for s in range(subjects):\n", + " S[s] = (np.random.rand(data[s].shape[0], samples) < p) * ((np.random.rand(data[s].shape[0], samples) * amplitude - amplitude/2) )\n", + " data[s] += S[s]" ] }, { @@ -107,12 +112,14 @@ "metadata": {}, "outputs": [], "source": [ - "srm = SRM(features=3, n_iter=20)\n", - "srm.fit(data)\n", + "algo_srm = srm.SRM(features=3, n_iter=20)\n", + "algo_srm.fit(data)\n", + "\n", "\n", + "algo_rsrm = rsrm.RSRM(features=3, gamma=0.35, n_iter=20)\n", + "algo_rsrm.fit(data)\n", "\n", - "rsrm = RSRM(features=3, gamma=0.35, n_iter=20)\n", - "rsrm.fit(data)" + "print('Done')" ] }, { @@ -128,7 +135,7 @@ "metadata": {}, "outputs": [], "source": [ - "def find_orthogonal_transform(shared_response):\n", + "def find_orthogonal_transform(shared_response, curve):\n", " u,_,vt = np.linalg.svd(shared_response.dot(curve.T))\n", " q = u.dot(vt)\n", " aligned_curve = q.T.dot(shared_response)\n", @@ -152,32 +159,33 @@ "ax = fig.gca(projection='3d')\n", "ax.plot(curve[0, :], curve[1, :], curve[2, :], '-g', label='original', lineWidth=5)\n", "\n", - "proj = find_orthogonal_transform(srm.s_)\n", + "proj = find_orthogonal_transform(algo_srm.s_, curve)\n", "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-b', label='SRM', lineWidth=3)\n", - "proj = find_orthogonal_transform(rsrm.r_)\n", - "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)" + "proj = find_orthogonal_transform(algo_rsrm.r_, curve)\n", + "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)\n", + "plt.legend()" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2.0 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" + "pygments_lexer": "ipython3", + "version": "3.6.3" } }, "nbformat": 4, - "nbformat_minor": 0 -} \ No newline at end of file + "nbformat_minor": 1 +} From 97d6c5089a152666e71b9023b79c417c30889aab Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Wed, 22 Nov 2017 21:54:55 -0800 Subject: [PATCH 14/32] Added news fragment --- docs/newsfragments/302.feature | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/newsfragments/302.feature diff --git a/docs/newsfragments/302.feature b/docs/newsfragments/302.feature new file mode 100644 index 000000000..16c36c4f7 --- /dev/null +++ b/docs/newsfragments/302.feature @@ -0,0 +1 @@ +Added the Robust Shared Response Model method. \ No newline at end of file From 10a3dc0fbc73d0f945d104099d538e9af7031122 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:30:00 -0800 Subject: [PATCH 15/32] Gamma doc on sparsity of solution --- brainiak/funcalign/rsrm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index a0597aa2b..4845080de 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -56,7 +56,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): gamma : float, default: 1.0 Regularization parameter for the sparseness of the individual - components. + components. Higher values yield sparser individual components. rand_seed : int, default: 0 Seed for initializing the random number generator. From 0f722c4131b614da688776b449802c03492294b5 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:34:27 -0800 Subject: [PATCH 16/32] Change Outlier to Individual everywhere --- brainiak/funcalign/rsrm.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 4845080de..bddef3be7 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -40,8 +40,8 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): """Robust Shared Response Model (RSRM) Given multi-subject data, factorize it as a shared response R among all - subjects, an orthogonal transform W per subject, and an outlying - (individual) sparse component S per subject: + subjects, an orthogonal transform W per subject, and an individual + (outlying) sparse component S per subject: .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N @@ -72,7 +72,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): The shared response. s_ : list of array, element i has shape=[voxels_i, samples] - The individual components with outlying data for each subject. + The individual components for each subject. random_state_: `RandomState` Random number generator initialized using rand_seed @@ -159,7 +159,7 @@ def transform(self, X, y=None): Shared responses from input data (X) s : list of 2D arrays, element i has shape=[voxels_i, samples_i] - Outlying data obtained from fitting model to input data (X) + Individual data obtained from fitting model to input data (X) """ # Check if the model exist if hasattr(self, 'w_') is False: @@ -181,7 +181,7 @@ def transform(self, X, y=None): def _transform_new_data(self, X, subject): """Transform new data for a subjects by projecting to the shared subspace and - computing the outlying individual information. + computing the individual information. Parameters ---------- @@ -199,7 +199,7 @@ def _transform_new_data(self, X, subject): Shared response from input data (X) S : array, shape=[voxels, samples] - Outlying data obtained from fitting model to input data (X) + Individual data obtained from fitting model to input data (X) """ S = np.zeros_like(X) R = None @@ -273,7 +273,7 @@ def _rsrm(self, X): # Initialization W = self._init_transforms(subjs, voxels, features, self.random_state_) - S = self._init_outliers(subjs, voxels, TRs) + S = self._init_individual(subjs, voxels, TRs) R = self._update_shared_response(X, S, W, subjs, features, TRs) if logger.isEnabledFor(logging.INFO): @@ -283,7 +283,7 @@ def _rsrm(self, X): # Main loop for i in range(self.n_iter): W = self._update_transforms(X, S, R, subjs) - S = self._update_outliers(X, W, R, subjs, self.lam) + S = self._update_individual(X, W, R, subjs, self.lam) R = self._update_shared_response(X, S, W, subjs, features, TRs) # Print objective function every iteration if logger.isEnabledFor(logging.INFO): @@ -377,8 +377,8 @@ def _objective_function(X, W, R, S, subjs, gamma): return func @staticmethod - def _update_outliers(X, W, R, subjs, gamma): - """Updates the outlier terms `S_i`. + def _update_individual(X, W, R, subjs, gamma): + """Update the individual components `S_i`. Parameters ---------- @@ -412,8 +412,8 @@ def _update_outliers(X, W, R, subjs, gamma): return S @staticmethod - def _init_outliers(subjs, voxels, TRs): - """Initializes the outlier terms `S_i` to empty (all zeros). + def _init_individual(subjs, voxels, TRs): + """Initializes the individual components `S_i` to empty (all zeros). Parameters ---------- From 276c1fd8f6143ad6c296c39965183158548bbd82 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:37:00 -0800 Subject: [PATCH 17/32] Changed samples by timepoints --- brainiak/funcalign/rsrm.py | 72 +++++++++++++++++++------------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index bddef3be7..9f3677054 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -68,10 +68,10 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): w_ : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) for each subject. - r_ : array, shape=[features, samples] + r_ : array, shape=[features, timepoints] The shared response. - s_ : list of array, element i has shape=[voxels_i, samples] + s_ : list of array, element i has shape=[voxels_i, timepoints] The individual components for each subject. random_state_: `RandomState` @@ -81,7 +81,7 @@ class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): ---- The number of voxels may be different between subjects. However, the - number of samples for the alignment data must be the same across + number of timepoints for the alignment data must be the same across subjects. The Robust Shared Response Model is approximated using the @@ -102,7 +102,7 @@ def fit(self, X, y=None): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data of one subject. y : not used @@ -121,7 +121,7 @@ def fit(self, X, y=None): # Check for input data sizes if X[0].shape[1] < self.features: raise ValueError( - "There are not enough samples to train the model with " + "There are not enough timepoints to train the model with " "{0:d} features.".format(self.features)) # Check if all subjects have same number of TRs for alignment @@ -130,7 +130,7 @@ def fit(self, X, y=None): for subject in range(number_subjects): assert_all_finite(X[subject]) if X[subject].shape[1] != number_trs: - raise ValueError("Different number of alignment samples " + raise ValueError("Different number of alignment timepoints " "between subjects.") # Create a new random state @@ -147,7 +147,7 @@ def transform(self, X, y=None): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples_i] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints_i] Each element in the list contains the fMRI data of one subject. y : not used @@ -155,10 +155,10 @@ def transform(self, X, y=None): Returns ------- - r : list of 2D arrays, element i has shape=[features_i, samples_i] + r : list of 2D arrays, element i has shape=[features_i, timepoints_i] Shared responses from input data (X) - s : list of 2D arrays, element i has shape=[voxels_i, samples_i] + s : list of 2D arrays, element i has shape=[voxels_i, timepoints_i] Individual data obtained from fitting model to input data (X) """ # Check if the model exist @@ -186,7 +186,7 @@ def _transform_new_data(self, X, subject): Parameters ---------- - X : array, shape=[voxels, samples] + X : array, shape=[voxels, timepoints] The fMRI data of the subject. subject : int @@ -195,10 +195,10 @@ def _transform_new_data(self, X, subject): Returns ------- - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] Shared response from input data (X) - S : array, shape=[voxels, samples] + S : array, shape=[voxels, timepoints] Individual data obtained from fitting model to input data (X) """ S = np.zeros_like(X) @@ -214,7 +214,7 @@ def transform_subject(self, X, y=None): Parameters ---------- - X : 2D array, shape=[voxels, samples] + X : 2D array, shape=[voxels, timepoints] The fMRI data of the new subject. y : not used @@ -225,7 +225,7 @@ def transform_subject(self, X, y=None): w : 2D array, shape=[voxels, features] Orthogonal mapping `W_{new}` for new subject - s : 2D array, shape=[voxels, samples] + s : 2D array, shape=[voxels, timepoints] Individual term `S_{new}` for new subject """ # Check if the model exist @@ -234,7 +234,7 @@ def transform_subject(self, X, y=None): # Check the number of TRs in the subject if X.shape[1] != self.r_.shape[1]: - raise ValueError("The number of samples(TRs) does not match the" + raise ValueError("The number of timepoints(TRs) does not match the" "one in the model.") s = np.zeros_like(X) @@ -250,7 +250,7 @@ def _rsrm(self, X): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject. @@ -260,10 +260,10 @@ def _rsrm(self, X): W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The shared response. - S : list of array, element i has shape=[voxels_i, samples] + S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. """ subjs = len(X) @@ -343,17 +343,17 @@ def _objective_function(X, W, R, S, subjs, gamma): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject. W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The shared response. - S : list of array, element i has shape=[voxels_i, samples] + S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. subjs : int @@ -383,14 +383,14 @@ def _update_individual(X, W, R, subjs, gamma): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject. W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The shared response. subjs : int @@ -403,7 +403,7 @@ def _update_individual(X, W, R, subjs, gamma): Returns ------- - S : list of array, element i has shape=[voxels_i, samples] + S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. """ S = [] @@ -425,12 +425,12 @@ def _init_individual(subjs, voxels, TRs): A list with the number of voxels per subject. TRs : int - The number of samples in the data. + The number of timepoints in the data. Returns ------- - S : list of 2D array, element i has shape=[voxels_i, samples] + S : list of 2D array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject initialized to zero. """ @@ -443,11 +443,11 @@ def _update_shared_response(X, S, W, subjs, features, TRs): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject. - S : list of array, element i has shape=[voxels_i, samples] + S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. W : list of array, element i has shape=[voxels_i, features] @@ -460,12 +460,12 @@ def _update_shared_response(X, S, W, subjs, features, TRs): The number of features in the model. TRs : int - The number of samples in the data. + The number of timepoints in the data. Returns ------- - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The updated shared response. """ R = np.zeros((features, TRs)) @@ -482,13 +482,13 @@ def _update_transform_subject(X, S, R): Parameters ---------- - X : array, shape=[voxels, samples] + X : array, shape=[voxels, timepoints] The fMRI data for aligning the subject. - S : array, shape=[voxels, samples] + S : array, shape=[voxels, timepoints] The individual component :math:`S_i` for the subject. - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The shared response. Returns @@ -509,14 +509,14 @@ def _update_transforms(X, S, R, subjs): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, samples] + X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject.ß - S : list of array, element i has shape=[voxels_i, samples] + S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. - R : array, shape=[features, samples] + R : array, shape=[features, timepoints] The shared response. subjs : int From 6c77d45ac029838674a21e27d0292b41d26a3f22 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:40:22 -0800 Subject: [PATCH 18/32] Removed unused parameters --- brainiak/funcalign/rsrm.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 9f3677054..9e14128c1 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -96,7 +96,7 @@ def __init__(self, n_iter=10, features=50, gamma=1.0, rand_seed=0): self.lam = gamma self.rand_seed = rand_seed - def fit(self, X, y=None): + def fit(self, X): """Compute the Robust Shared Response Model Parameters @@ -104,8 +104,6 @@ def fit(self, X, y=None): X : list of 2D arrays, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data of one subject. - - y : not used """ logger.info('Starting RSRM') @@ -141,7 +139,7 @@ def fit(self, X, y=None): return self - def transform(self, X, y=None): + def transform(self, X): """Use the model to transform new data to Shared Response space Parameters @@ -150,8 +148,6 @@ def transform(self, X, y=None): X : list of 2D arrays, element i has shape=[voxels_i, timepoints_i] Each element in the list contains the fMRI data of one subject. - y : not used - Returns ------- @@ -208,7 +204,7 @@ def _transform_new_data(self, X, subject): S = self._shrink(X - self.w_[subject].dot(R), self.lam) return R, S - def transform_subject(self, X, y=None): + def transform_subject(self, X): """Transform a new subject using the existing model Parameters @@ -217,8 +213,6 @@ def transform_subject(self, X, y=None): X : 2D array, shape=[voxels, timepoints] The fMRI data of the new subject. - y : not used - Returns ------- From 74b73567bc7566f22c82ad2c69e1d48f63271cfe Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:42:42 -0800 Subject: [PATCH 19/32] orthogonalization comment --- brainiak/funcalign/rsrm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 9e14128c1..d598a3b26 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -322,6 +322,7 @@ def _init_transforms(self, subjs, voxels, features, random_state): # Draw a random W for each subject W = [random_state.random_sample((voxels[i], features)) for i in range(subjs)] + # Make it orthogonal it with QR decomposition for i in range(subjs): W[i], _ = np.linalg.qr(W[i]) From 5974b50e623b27452b6fefef642271292a233699 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 10:47:36 -0800 Subject: [PATCH 20/32] Removed subjs from obj_func() --- brainiak/funcalign/rsrm.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index d598a3b26..5286928df 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -271,7 +271,7 @@ def _rsrm(self, X): R = self._update_shared_response(X, S, W, subjs, features, TRs) if logger.isEnabledFor(logging.INFO): - objective = self._objective_function(X, W, R, S, subjs, self.lam) + objective = self._objective_function(X, W, R, S, self.lam) logger.info('Objective function %f' % objective) # Main loop @@ -281,8 +281,7 @@ def _rsrm(self, X): R = self._update_shared_response(X, S, W, subjs, features, TRs) # Print objective function every iteration if logger.isEnabledFor(logging.INFO): - objective = self._objective_function(X, W, R, S, subjs, - self.lam) + objective = self._objective_function(X, W, R, S, self.lam) logger.info('Objective function %f' % objective) return W, R, S @@ -329,7 +328,7 @@ def _init_transforms(self, subjs, voxels, features, random_state): return W @staticmethod - def _objective_function(X, W, R, S, subjs, gamma): + def _objective_function(X, W, R, S, gamma): """Evaluate the objective function. .. math:: \sum_{i=1}^{N} 1/2 \| X_i - W_i R - S_i \|_F^2 @@ -338,7 +337,7 @@ def _objective_function(X, W, R, S, subjs, gamma): Parameters ---------- - X : list of 2D arrays, element i has shape=[voxels_i, timepoints] + X : list of array, element i has shape=[voxels_i, timepoints] Each element in the list contains the fMRI data for alignment of one subject. @@ -351,9 +350,6 @@ def _objective_function(X, W, R, S, subjs, gamma): S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. - subjs : int - The number of subjects. - gamma : float, default: 1.0 Regularization parameter for the sparseness of the individual components. @@ -365,6 +361,7 @@ def _objective_function(X, W, R, S, subjs, gamma): The RSRM objective function evaluated on the parameters to this function. """ + subjs = len(X) func = .0 for i in range(subjs): func += 0.5 * np.sum((X[i] - W[i].dot(R) - S[i])**2) \ From 2870d4460b16f7da5257179c950f4e298b72604c Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:04:20 -0800 Subject: [PATCH 21/32] Removed subjs and TRs from unnecessary places, and doc on update of R added --- brainiak/funcalign/rsrm.py | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 5286928df..b44b7a784 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -268,7 +268,7 @@ def _rsrm(self, X): # Initialization W = self._init_transforms(subjs, voxels, features, self.random_state_) S = self._init_individual(subjs, voxels, TRs) - R = self._update_shared_response(X, S, W, subjs, features, TRs) + R = self._update_shared_response(X, S, W, features) if logger.isEnabledFor(logging.INFO): objective = self._objective_function(X, W, R, S, self.lam) @@ -276,9 +276,9 @@ def _rsrm(self, X): # Main loop for i in range(self.n_iter): - W = self._update_transforms(X, S, R, subjs) - S = self._update_individual(X, W, R, subjs, self.lam) - R = self._update_shared_response(X, S, W, subjs, features, TRs) + W = self._update_transforms(X, S, R) + S = self._update_individual(X, W, R, self.lam) + R = self._update_shared_response(X, S, W, features) # Print objective function every iteration if logger.isEnabledFor(logging.INFO): objective = self._objective_function(X, W, R, S, self.lam) @@ -429,7 +429,7 @@ def _init_individual(subjs, voxels, TRs): return [np.zeros((voxels[i], TRs)) for i in range(subjs)] @staticmethod - def _update_shared_response(X, S, W, subjs, features, TRs): + def _update_shared_response(X, S, W, features): """Update the shared response `R`. Parameters @@ -445,25 +445,22 @@ def _update_shared_response(X, S, W, subjs, features, TRs): W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. - subjs : int - The number of subjects. - features : int The number of features in the model. - TRs : int - The number of timepoints in the data. - Returns ------- R : array, shape=[features, timepoints] The updated shared response. """ + subjs = len(X) + TRs = X[0].shape[1] R = np.zeros((features, TRs)) + # Project the subject data with the individual component removed into + # the shared subspace and average over all subjects. for i in range(subjs): - R += W[i].T.dot(X[i]) - R -= W[i].T.dot(S[i]) + R += W[i].T.dot(X[i]-S[i]) R /= subjs return R @@ -495,7 +492,7 @@ def _update_transform_subject(X, S, R): return U.dot(V) @staticmethod - def _update_transforms(X, S, R, subjs): + def _update_transforms(X, S, R): """Updates the mappings `W_i` for each subject. Parameters @@ -511,15 +508,13 @@ def _update_transforms(X, S, R, subjs): R : array, shape=[features, timepoints] The shared response. - subjs : int - The number of subjects. - Returns ------- W : list of array, element i has shape=[voxels_i, features] The orthogonal transforms (mappings) :math:`W_i` for each subject. """ + subjs = len(X) W = [] for i in range(subjs): W.append(RSRM._update_transform_subject(X[i], S[i], R)) From 4b2f071671e95cce8eaf76bd56118504feb0c6d2 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:06:01 -0800 Subject: [PATCH 22/32] Changed variables to show specific subject treatment --- brainiak/funcalign/rsrm.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index b44b7a784..528e7a64d 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -465,16 +465,16 @@ def _update_shared_response(X, S, W, features): return R @staticmethod - def _update_transform_subject(X, S, R): + def _update_transform_subject(Xi, Si, R): """Updates the mappings `W_i` for one subject. Parameters ---------- - X : array, shape=[voxels, timepoints] - The fMRI data for aligning the subject. + Xi : array, shape=[voxels, timepoints] + The fMRI data :math:`X_i` for aligning the subject. - S : array, shape=[voxels, timepoints] + Si : array, shape=[voxels, timepoints] The individual component :math:`S_i` for the subject. R : array, shape=[features, timepoints] @@ -483,11 +483,11 @@ def _update_transform_subject(X, S, R): Returns ------- - W : array, shape=[voxels, features] + Wi : array, shape=[voxels, features] The orthogonal transform (mapping) :math:`W_i` for the subject. """ - A = X.dot(R.T) - A -= S.dot(R.T) + A = Xi.dot(R.T) + A -= Si.dot(R.T) U, _, V = np.linalg.svd(A, full_matrices=False) return U.dot(V) From ae5ca6f6f82d30f01c2649e3e8f1b8ea50e4daf6 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:08:34 -0800 Subject: [PATCH 23/32] Documented transform computation --- brainiak/funcalign/rsrm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 528e7a64d..a76fc2c8a 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -488,6 +488,7 @@ def _update_transform_subject(Xi, Si, R): """ A = Xi.dot(R.T) A -= Si.dot(R.T) + # Solve the Procrustes problem U, _, V = np.linalg.svd(A, full_matrices=False) return U.dot(V) From e9c309877cacdd9e90766bcc9d49da73c82d8bce Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:38:35 -0800 Subject: [PATCH 24/32] Removed unnecessary Classifier Mixin --- brainiak/funcalign/rsrm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index a76fc2c8a..4f8ebd74a 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -25,7 +25,7 @@ import logging import numpy as np -from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin +from sklearn.base import BaseEstimator, TransformerMixin from sklearn.utils import assert_all_finite from sklearn.utils.validation import NotFittedError @@ -36,14 +36,14 @@ logger = logging.getLogger(__name__) -class RSRM(BaseEstimator, ClassifierMixin, TransformerMixin): +class RSRM(BaseEstimator, TransformerMixin): """Robust Shared Response Model (RSRM) Given multi-subject data, factorize it as a shared response R among all subjects, an orthogonal transform W per subject, and an individual (outlying) sparse component S per subject: - .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N + .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N Parameters ---------- From 6ff0682e4575a43d483f376d9f9114e12d405e18 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:50:07 -0800 Subject: [PATCH 25/32] Method documentation extended --- brainiak/funcalign/rsrm.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 4f8ebd74a..51e605b12 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -45,6 +45,20 @@ class RSRM(BaseEstimator, TransformerMixin): .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N + This unsupervised model allows to learn idiosyncratic information for + subjects and simultaneously improve the shared response estimation. + The model has similar properties to the Shared Response Model (SRM) with + the addition of the individual components. + + The model is estimated solving the following optimization problem: + + .. math:: \\min_{W_i, S_i, R}\\sum_i 1/2\|X_i - W_i R - S_i\|_F^2 + .. math:: + \\gamma\|S_i\|_1 + .. math:: s.t. \\qquad W_i^TW_i = I \\forall i=1 \\dots N + + The solution to this problem is obtained by applying a Block-Coordinate + Descent procedure. More details can be found in [Turek2017]_. + Parameters ---------- From 3a690efa5882cbbc9bc7febffc8f22d2c1a68c88 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 11:55:53 -0800 Subject: [PATCH 26/32] removed trailing spaces --- brainiak/funcalign/rsrm.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 51e605b12..40eb109b4 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -40,25 +40,25 @@ class RSRM(BaseEstimator, TransformerMixin): """Robust Shared Response Model (RSRM) Given multi-subject data, factorize it as a shared response R among all - subjects, an orthogonal transform W per subject, and an individual + subjects, an orthogonal transform W per subject, and an individual (outlying) sparse component S per subject: - .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N + .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N - This unsupervised model allows to learn idiosyncratic information for - subjects and simultaneously improve the shared response estimation. - The model has similar properties to the Shared Response Model (SRM) with + This unsupervised model allows to learn idiosyncratic information for + subjects and simultaneously improve the shared response estimation. + The model has similar properties to the Shared Response Model (SRM) with the addition of the individual components. - + The model is estimated solving the following optimization problem: - - .. math:: \\min_{W_i, S_i, R}\\sum_i 1/2\|X_i - W_i R - S_i\|_F^2 + + .. math:: \\min_{W_i, S_i, R}\\sum_i 1/2\|X_i - W_i R - S_i\|_F^2 .. math:: + \\gamma\|S_i\|_1 .. math:: s.t. \\qquad W_i^TW_i = I \\forall i=1 \\dots N - - The solution to this problem is obtained by applying a Block-Coordinate + + The solution to this problem is obtained by applying a Block-Coordinate Descent procedure. More details can be found in [Turek2017]_. - + Parameters ---------- From 2563fb108b1ff173bffd6776510639e43e14d810 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 12:01:10 -0800 Subject: [PATCH 27/32] removed unnecessary code --- brainiak/funcalign/rsrm.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 40eb109b4..772aa6190 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -383,7 +383,7 @@ def _objective_function(X, W, R, S, gamma): return func @staticmethod - def _update_individual(X, W, R, subjs, gamma): + def _update_individual(X, W, R, gamma): """Update the individual components `S_i`. Parameters @@ -399,9 +399,6 @@ def _update_individual(X, W, R, subjs, gamma): R : array, shape=[features, timepoints] The shared response. - subjs : int - The number of subjects. - gamma : float, default: 1.0 Regularization parameter for the sparseness of the individual components. @@ -412,6 +409,7 @@ def _update_individual(X, W, R, subjs, gamma): S : list of array, element i has shape=[voxels_i, timepoints] The individual component :math:`S_i` for each subject. """ + subjs = len(X) S = [] for i in range(subjs): S.append(RSRM._shrink(X[i] - W[i].dot(R), gamma)) From 6d6825cf06fa623086d7f7de0a0cc18a7072c7c7 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 12:27:02 -0800 Subject: [PATCH 28/32] corrected latex --- brainiak/funcalign/rsrm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 772aa6190..6ec4d0c6c 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -54,7 +54,7 @@ class RSRM(BaseEstimator, TransformerMixin): .. math:: \\min_{W_i, S_i, R}\\sum_i 1/2\|X_i - W_i R - S_i\|_F^2 .. math:: + \\gamma\|S_i\|_1 - .. math:: s.t. \\qquad W_i^TW_i = I \\forall i=1 \\dots N + .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N The solution to this problem is obtained by applying a Block-Coordinate Descent procedure. More details can be found in [Turek2017]_. From 61b9d149f4d4d51f2d04d1283d612ec39bcb7892 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Mon, 27 Nov 2017 12:31:14 -0800 Subject: [PATCH 29/32] corrected another latex line --- brainiak/funcalign/rsrm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 6ec4d0c6c..34efac464 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -52,7 +52,7 @@ class RSRM(BaseEstimator, TransformerMixin): The model is estimated solving the following optimization problem: - .. math:: \\min_{W_i, S_i, R}\\sum_i 1/2\|X_i - W_i R - S_i\|_F^2 + .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 .. math:: + \\gamma\|S_i\|_1 .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N From 431c13015b475e54fe956b229cbb0f91af283846 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Tue, 28 Nov 2017 09:22:53 -0800 Subject: [PATCH 30/32] corrected math latex problem in multiple lines --- brainiak/funcalign/rsrm.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 34efac464..f903dfdb8 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -52,8 +52,7 @@ class RSRM(BaseEstimator, TransformerMixin): The model is estimated solving the following optimization problem: - .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 - .. math:: + \\gamma\|S_i\|_1 + .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 + \\gamma\|S_i\|_1 # noqa .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N The solution to this problem is obtained by applying a Block-Coordinate From e7df74471e4a2bfbf0caebd037d315d632a47091 Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Tue, 5 Dec 2017 17:51:04 -0800 Subject: [PATCH 31/32] recorrected latex --- brainiak/funcalign/rsrm.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index f903dfdb8..3df9abdc0 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -52,7 +52,8 @@ class RSRM(BaseEstimator, TransformerMixin): The model is estimated solving the following optimization problem: - .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 + \\gamma\|S_i\|_1 # noqa + .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 + .. math:: + \\gamma\|S_i\|_1 .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N The solution to this problem is obtained by applying a Block-Coordinate From 49b96fa03267b150304b9f23496911b874f457af Mon Sep 17 00:00:00 2001 From: Javier Turek Date: Tue, 5 Dec 2017 17:56:15 -0800 Subject: [PATCH 32/32] removed trailing spaces --- brainiak/funcalign/rsrm.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py index 3df9abdc0..34efac464 100644 --- a/brainiak/funcalign/rsrm.py +++ b/brainiak/funcalign/rsrm.py @@ -52,8 +52,8 @@ class RSRM(BaseEstimator, TransformerMixin): The model is estimated solving the following optimization problem: - .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 - .. math:: + \\gamma\|S_i\|_1 + .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2 + .. math:: + \\gamma\|S_i\|_1 .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N The solution to this problem is obtained by applying a Block-Coordinate