Skip to content

srm: subtract the subject-specific mean in srm.transform() and srm.transform_subject() #495

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
qihongl wants to merge 6 commits into
base: master
Choose a base branch
from
Open

srm: subtract the subject-specific mean in srm.transform() and srm.transform_subject() #495

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
cdc53fa
Merge pull request #2 from brainiak/master
qihongl Apr 9, 2021
41ccf04
init: srm - use the intercept term
qihongl Apr 10, 2021
ebb7047
add intercept for detsrm
qihongl Apr 21, 2021
db338e9
actually it seems okay to overwrite 'data' with 'centered_data', whic…
qihongl Apr 21, 2021
8f810da
oops, fix other tests that uses transform_subject, which was expectin…
qihongl Apr 21, 2021
6d899e8
oops, remove stuff that i commented out and forgot to delete
qihongl Apr 30, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions brainiak/funcalign/rsrm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,8 @@ class RSRM(BaseEstimator, TransformerMixin):

The number of voxels may be different between subjects. However, the
number of timepoints for the alignment data must be the same across
subjects.
subjects. Note that unlike SRM, DetSRM does not handle vector shifts
(intercepts) across subjects.

The Robust Shared Response Model is approximated using the
Block-Coordinate Descent (BCD) algorithm proposed in [Turek2017]_.
Expand Down Expand Up @@ -380,7 +381,7 @@ def _objective_function(X, W, R, S, gamma):
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]))
+ gamma * np.sum(np.abs(S[i]))
return func

@staticmethod
Expand Down Expand Up @@ -473,7 +474,7 @@ def _update_shared_response(X, S, W, features):
# 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]-S[i])
R += W[i].T.dot(X[i] - S[i])
R /= subjs
return R

Expand Down
49 changes: 35 additions & 14 deletions brainiak/funcalign/srm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,8 +296,8 @@ def transform(self, X, y=None):
s = [None] * len(X)
for subject in range(len(X)):
if X[subject] is not None:
s[subject] = self.w_[subject].T.dot(X[subject])

s[subject] = self.w_[subject].T.dot(
X[subject] - self.mu_[subject][:, np.newaxis])
return s

def _init_structures(self, data, subjects):
Expand Down Expand Up @@ -419,7 +419,10 @@ def _update_transform_subject(Xi, S):

def transform_subject(self, X):
"""Transform a new subject using the existing model.
The subject is assumed to have recieved equivalent stimulation
The subject is assumed to have recieved equivalent stimulation. In
particular, to transform the new subject X with w and mu, one can do
the following:
shared_X = w.T @ (X - mu[:, np.newaxis])

Parameters
----------
Expand All @@ -432,6 +435,8 @@ def transform_subject(self, X):

w : 2D array, shape=[voxels, features]
Orthogonal mapping `W_{new}` for new subject
mu : 1D array, shape=[voxels]
The voxel means for the new subject

"""
# Check if the model exist
Expand All @@ -442,10 +447,10 @@ def transform_subject(self, X):
if X.shape[1] != self.s_.shape[1]:
raise ValueError("The number of timepoints(TRs) does not match the"
"one in the model.")

w = self._update_transform_subject(X, self.s_)

return w
# get the intercept for mean centering, as procrustes doesn't handle it
mu = np.mean(X, axis=1)
w = self._update_transform_subject(X - mu[:, np.newaxis], self.s_)
return w, mu

def save(self, file):
"""Save fitted SRM to .npz file.
Expand Down Expand Up @@ -562,7 +567,7 @@ def _srm(self, data):
for subject in range(subjects):
if data[subject] is not None:
wt_invpsi_x += (w[subject].T.dot(x[subject])) \
/ rho2[subject]
/ rho2[subject]
trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject]

wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM)
Expand Down Expand Up @@ -650,6 +655,9 @@ class DetSRM(BaseEstimator, TransformerMixin):
s_ : array, shape=[features, samples]
The shared response.

mu_ : list of array, element i has shape=[voxels_i]
The voxel means over the samples for each subject.

random_state_: `RandomState`
Random number generator initialized using rand_seed

Expand Down Expand Up @@ -742,8 +750,8 @@ def transform(self, X, y=None):

s = [None] * len(X)
for subject in range(len(X)):
s[subject] = self.w_[subject].T.dot(X[subject])

s[subject] = self.w_[subject].T.dot(
X[subject] - self.mu_[subject][:, np.newaxis])
return s

def _objective_function(self, data, w, s):
Expand Down Expand Up @@ -818,11 +826,16 @@ def _update_transform_subject(Xi, S):

Wi : array, shape=[voxels, features]
The orthogonal transform (mapping) :math:`W_i` for the subject.

mu : 1D array, shape=[voxels]
The voxel means for the new subject
"""
A = Xi.dot(S.T)
# estimate the intercept and center the data
mu = np.mean(Xi, axis=1)
A = (Xi - mu[:, np.newaxis]).dot(S.T)
# Solve the Procrustes problem
U, _, V = np.linalg.svd(A, full_matrices=False)
return U.dot(V)
return U.dot(V), mu

def transform_subject(self, X):
"""Transform a new subject using the existing model.
Copy link

@TuKo TuKo Oct 25, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As the function _update_transform_subject() was modified to return mu, then we should update the documentation or code of transform_subject().

Copy link
Member Author

@qihongl qihongl Oct 29, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks! what do you mean? I did update the doc string of transform_subject() about mu

Expand Down Expand Up @@ -880,13 +893,20 @@ def _srm(self, data):
np.random.RandomState(self.random_state_.randint(2 ** 32))
for i in range(len(data))]

# compute subject specific intercept
self.mu_ = [np.mean(data[s], axis=1) for s in range(subjects)]
# center the data
data = [data[s] - self.mu_[s][:, np.newaxis]
for s in range(subjects)]

# Initialization step: initialize the outputs with initial values,
# voxels with the number of voxels in each subject.
w, _ = _init_w_transforms(data, self.features, random_states)
shared_response = self._compute_shared_response(data, w)
if logger.isEnabledFor(logging.INFO):
# Calculate the current objective function value
objective = self._objective_function(data, w, shared_response)
objective = self._objective_function(
data, w, shared_response)
logger.info('Objective function %f' % objective)

# Main loop of the algorithm
Expand All @@ -907,7 +927,8 @@ def _srm(self, data):

if logger.isEnabledFor(logging.INFO):
# Calculate the current objective function value
objective = self._objective_function(data, w, shared_response)
objective = self._objective_function(
data, w, shared_response)
logger.info('Objective function %f' % objective)

return w, shared_response
2 changes: 2 additions & 0 deletions brainiak/funcalign/sssrm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,8 @@ class SSSRM(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
subjects. The number of labeled samples per subject can be different.
Note that unlike SRM, DetSRM does not handle vector shifts (intercepts)
across subjects.

The Semi-Supervised Shared Response Model is approximated using the
Block-Coordinate Descent (BCD) algorithm proposed in [Turek2016]_.
Expand Down
Loading