Skip to content

Commit

Permalink
GPFA renamed xsm -> latent_variable
Browse files Browse the repository at this point in the history
  • Loading branch information
dizcza committed Nov 11, 2020
1 parent 6f81f6c commit 295c6bd
Show file tree
Hide file tree
Showing 4 changed files with 65 additions and 48 deletions.
36 changes: 23 additions & 13 deletions elephant/gpfa/gpfa.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""
Gaussian-process factor analysis (GPFA) is a dimensionality reduction method
[#f1]_ for neural trajectory visualization of parallel spike trains. GPFA applies
factor analysis (FA) to time-binned spike count data to reduce the
[#f1]_ for neural trajectory visualization of parallel spike trains. GPFA
applies factor analysis (FA) to time-binned spike count data to reduce the
dimensionality and at the same time smoothes the resulting low-dimensional
trajectories by fitting a Gaussian process (GP) model to them.
Expand Down Expand Up @@ -198,13 +198,16 @@ class GPFA(sklearn.base.BaseEstimator):
...
>>> gpfa = GPFA(bin_size=20*pq.ms, x_dim=8)
>>> gpfa.fit(data)
>>> results = gpfa.transform(data, returned_data=['xorth', 'xsm'])
>>> xorth = results['xorth']; xsm = results['xsm']
>>> results = gpfa.transform(data, returned_data=['latent_variable_orth',
... 'latent_variable'])
>>> latent_variable_orth = results['latent_variable_orth']
>>> latent_variable = results['latent_variable']
or simply
>>> results = GPFA(bin_size=20*pq.ms, x_dim=8).fit_transform(data,
... returned_data=['xorth', 'xsm'])
... returned_data=['latent_variable_orth',
... 'latent_variable'])
"""

@deprecated_alias(binsize='bin_size')
Expand All @@ -219,7 +222,12 @@ def __init__(self, bin_size=20 * pq.ms, x_dim=3, min_var_frac=0.01,
self.em_tol = em_tol
self.em_max_iters = em_max_iters
self.freq_ll = freq_ll
self.valid_data_names = ('xorth', 'xsm', 'Vsm', 'VsmGP', 'y')
self.valid_data_names = (
'latent_variable_orth',
'latent_variable',
'Vsm',
'VsmGP',
'y')
self.verbose = verbose

if not isinstance(self.bin_size, pq.Quantity):
Expand Down Expand Up @@ -317,7 +325,7 @@ def _format_training_data(self, spiketrains):
seq['y'] = seq['y'][self.has_spikes_bool, :]
return seqs

def transform(self, spiketrains, returned_data=['xorth']):
def transform(self, spiketrains, returned_data=['latent_variable_orth']):
"""
Obtain trajectories of neural activity in a low-dimensional latent
variable space by inferring the posterior mean of the obtained GPFA
Expand All @@ -338,9 +346,10 @@ def transform(self, spiketrains, returned_data=['xorth']):
The dimensionality reduction transform generates the following
resultant data:
'xorth': orthonormalized posterior mean of latent variable
'latent_variable_orth': orthonormalized posterior mean of latent
variable
'xsm': posterior mean of latent variable before
'latent_variable': posterior mean of latent variable before
orthonormalization
'Vsm': posterior covariance between latent variables
Expand All @@ -352,7 +361,7 @@ def transform(self, spiketrains, returned_data=['xorth']):
`returned_data` specifies the keys by which the data dict is
returned.
Default is ['xorth'].
Default is ['latent_variable_orth'].
Returns
-------
Expand All @@ -367,9 +376,9 @@ def transform(self, spiketrains, returned_data=['xorth']):
shape, specific to each data type, containing the corresponding
data for the n-th trial:
`xorth`: (#latent_vars, #bins) np.ndarray
`latent_variable_orth`: (#latent_vars, #bins) np.ndarray
`xsm`: (#latent_vars, #bins) np.ndarray
`latent_variable`: (#latent_vars, #bins) np.ndarray
`y`: (#units, #bins) np.ndarray
Expand Down Expand Up @@ -410,7 +419,8 @@ def transform(self, spiketrains, returned_data=['xorth']):
return seqs[returned_data[0]]
return {x: seqs[x] for x in returned_data}

def fit_transform(self, spiketrains, returned_data=['xorth']):
def fit_transform(self, spiketrains, returned_data=[
'latent_variable_orth']):
"""
Fit the model with `spiketrains` data and apply the dimensionality
reduction on `spiketrains`.
Expand Down
41 changes: 22 additions & 19 deletions elephant/gpfa/gpfa_core.py
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ def em(params_init, seqs_train, max_iters=500, tol=1.0E-8, min_var_frac=0.01,
seqs_latent : np.recarray
a copy of the training data structure, augmented with the new
fields:
xsm : np.ndarray of shape (#latent_vars x #bins)
latent_variable : np.ndarray of shape (#latent_vars x #bins)
posterior mean of latent variables at each time bin
Vsm : np.ndarray of shape (#latent_vars, #latent_vars, #bins)
posterior covariance between latent variables at each
Expand Down Expand Up @@ -244,11 +244,12 @@ def em(params_init, seqs_train, max_iters=500, tol=1.0E-8, min_var_frac=0.01,
sum_p_auto = np.zeros((x_dim, x_dim))
for seq_latent in seqs_latent:
sum_p_auto += seq_latent['Vsm'].sum(axis=2) \
+ seq_latent['xsm'].dot(seq_latent['xsm'].T)
+ seq_latent['latent_variable'].dot(
seq_latent['latent_variable'].T)
y = np.hstack(seqs_train['y'])
xsm = np.hstack(seqs_latent['xsm'])
sum_yxtrans = y.dot(xsm.T)
sum_xall = xsm.sum(axis=1)[:, np.newaxis]
latent_variable = np.hstack(seqs_latent['latent_variable'])
sum_yxtrans = y.dot(latent_variable.T)
sum_xall = latent_variable.sum(axis=1)[:, np.newaxis]
sum_yall = y.sum(axis=1)[:, np.newaxis]

# term is (xDim+1) x (xDim+1)
Expand All @@ -262,8 +263,8 @@ def em(params_init, seqs_train, max_iters=500, tol=1.0E-8, min_var_frac=0.01,

# yCent must be based on the new d
# yCent = bsxfun(@minus, [seq.y], currentParams.d);
# R = (yCent * yCent' - (yCent * [seq.xsm]') * currentParams.C')
# / sum(T);
# R = (yCent * yCent' - (yCent * [seq.latent_variable]') * \
# currentParams.C') / sum(T);
c = params['C']
d = params['d'][:, np.newaxis]
if params['notes']['RforceDiagonal']:
Expand Down Expand Up @@ -344,7 +345,7 @@ def exact_inference_with_ll(seqs, params, get_ll=True):
seqs_latent : np.recarray
a copy of the input data structure, augmented with the new
fields:
xsm : (#latent_vars, #bins) np.ndarray
latent_variable : (#latent_vars, #bins) np.ndarray
posterior mean of latent variables at each time bin
Vsm : (#latent_vars, #latent_vars, #bins) np.ndarray
posterior covariance between latent variables at each
Expand All @@ -359,7 +360,7 @@ def exact_inference_with_ll(seqs, params, get_ll=True):

# copy the contents of the input data structure to output structure
dtype_out = [(x, seqs[x].dtype) for x in seqs.dtype.names]
dtype_out.extend([('xsm', np.object), ('Vsm', np.object),
dtype_out.extend([('latent_variable', np.object), ('Vsm', np.object),
('VsmGP', np.object)])
seqs_latent = np.empty(len(seqs), dtype=dtype_out)
for dtype_name in seqs.dtype.names:
Expand Down Expand Up @@ -424,12 +425,13 @@ def exact_inference_with_ll(seqs, params, get_ll=True):
blk_prod = k_big[:x_dim * t_half, :].dot(
gpfa_util.fill_persymm(np.eye(x_dim * t_half, x_dim * t) -
blk_prod, x_dim, t))
# xsmMat is (xDim*T) x length(nList)
xsm_mat = gpfa_util.fill_persymm(blk_prod, x_dim, t).dot(term1_mat)
# latent_variableMat is (xDim*T) x length(nList)
latent_variable_mat = gpfa_util.fill_persymm(
blk_prod, x_dim, t).dot(term1_mat)

for i, n in enumerate(n_list):
seqs_latent[n]['xsm'] = xsm_mat[:, i].reshape((x_dim, t),
order='F')
seqs_latent[n]['latent_variable'] = \
latent_variable_mat[:, i].reshape((x_dim, t), order='F')
seqs_latent[n]['Vsm'] = vsm
seqs_latent[n]['VsmGP'] = vsm_gp

Expand Down Expand Up @@ -536,7 +538,7 @@ def orthonormalize(params_est, seqs):
number of timesteps
y : np.ndarray of shape (#units, #bins)
neural data
xsm : np.ndarray of shape (#latent_vars, #bins)
latent_variable : np.ndarray of shape (#latent_vars, #bins)
posterior mean of latent variables at each time bin
Vsm : np.ndarray of shape (#latent_vars, #latent_vars, #bins)
posterior covariance between latent variables at each
Expand All @@ -550,13 +552,14 @@ def orthonormalize(params_est, seqs):
Estimated model parameters, including `Corth`, obtained by
orthonormalizing the columns of C.
seqs : np.recarray
Training data structure that contains the new field `xorth`,
the orthonormalized neural trajectories.
Training data structure that contains the new field
`latent_variable_orth`, the orthonormalized neural trajectories.
"""
C = params_est['C']
X = np.hstack(seqs['xsm'])
Xorth, Corth, _ = gpfa_util.orthonormalize(X, C)
seqs = gpfa_util.segment_by_trial(seqs, Xorth, 'xorth')
X = np.hstack(seqs['latent_variable'])
latent_variable_orth, Corth, _ = gpfa_util.orthonormalize(X, C)
seqs = gpfa_util.segment_by_trial(
seqs, latent_variable_orth, 'latent_variable_orth')

params_est['Corth'] = Corth

Expand Down
12 changes: 6 additions & 6 deletions elephant/gpfa/gpfa_util.py
Original file line number Diff line number Diff line change
Expand Up @@ -423,8 +423,8 @@ def make_precomp(seqs, xDim):
# Loop once for each trial (each of nList)
for n in precomp[i]['Tu'][j]['nList']:
precomp[i]['Tu'][j]['PautoSUM'] += seqs[n]['VsmGP'][:, :, i] \
+ np.outer(
seqs[n]['xsm'][i, :], seqs[n]['xsm'][i, :])
+ np.outer(seqs[n]['latent_variable'][i, :],
seqs[n]['latent_variable'][i, :])
return precomp


Expand Down Expand Up @@ -512,7 +512,7 @@ def orthonormalize(x, l):
Returns
-------
Xorth : (xDim, T) np.ndarray
latent_variable_orth : (xDim, T) np.ndarray
Orthonormalized latent variables
Lorth : (yDim, xDim) np.ndarray
Orthonormalized loading matrix
Expand All @@ -523,15 +523,15 @@ def orthonormalize(x, l):
if xDim == 1:
TT = np.sqrt(np.dot(l.T, l))
Lorth = rdiv(l, TT)
Xorth = np.dot(TT, x)
latent_variable_orth = np.dot(TT, x)
else:
UU, DD, VV = sp.linalg.svd(l, full_matrices=False)
# TT is transform matrix
TT = np.dot(np.diag(DD), VV)

Lorth = UU
Xorth = np.dot(TT, x)
return Xorth, Lorth, TT
latent_variable_orth = np.dot(TT, x)
return latent_variable_orth, Lorth, TT


def segment_by_trial(seqs, x, fn):
Expand Down
24 changes: 14 additions & 10 deletions elephant/test/test_gpfa.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,17 +91,19 @@ def gen_test_data(rates, durs, shapes=(1, 1, 1, 1)):
def test_data1(self):
gpfa = GPFA(x_dim=self.x_dim, em_max_iters=self.n_iters)
gpfa.fit(self.data1)
xorth = gpfa.transform(self.data1)
latent_variable_orth = gpfa.transform(self.data1)
self.assertAlmostEqual(gpfa.transform_info['log_likelihood'],
-8172.004695554373, places=5)
# Since data1 is inherently 2 dimensional, only the first two
# dimensions of xorth should have finite power.
# dimensions of latent_variable_orth should have finite power.
for i in [0, 1]:
self.assertNotEqual(xorth[0][i].mean(), 0)
self.assertNotEqual(xorth[0][i].var(), 0)
self.assertNotEqual(latent_variable_orth[0][i].mean(), 0)
self.assertNotEqual(latent_variable_orth[0][i].var(), 0)
for i in [2, 3]:
self.assertAlmostEqual(xorth[0][i].mean(), 0, places=2)
self.assertAlmostEqual(xorth[0][i].var(), 0, places=2)
self.assertAlmostEqual(latent_variable_orth[0][i].mean(), 0,
places=2)
self.assertAlmostEqual(latent_variable_orth[0][i].var(), 0,
places=2)

def test_transform_testing_data(self):
# check if the num. of neurons in the test data matches the
Expand Down Expand Up @@ -169,12 +171,14 @@ def test_fit_transform(self):
gpfa1 = GPFA(bin_size=self.bin_size, x_dim=self.x_dim,
em_max_iters=self.n_iters)
gpfa1.fit(self.data1)
xorth1 = gpfa1.transform(self.data1)
xorth2 = GPFA(bin_size=self.bin_size, x_dim=self.x_dim,
em_max_iters=self.n_iters).fit_transform(self.data1)
latent_variable_orth1 = gpfa1.transform(self.data1)
latent_variable_orth2 = GPFA(
bin_size=self.bin_size, x_dim=self.x_dim,
em_max_iters=self.n_iters).fit_transform(self.data1)
for i in range(len(self.data1)):
for j in range(self.x_dim):
assert_array_almost_equal(xorth1[i][j], xorth2[i][j])
assert_array_almost_equal(latent_variable_orth1[i][j],
latent_variable_orth2[i][j])

def test_get_seq_sqrt(self):
data = [self.data2[0]]
Expand Down

0 comments on commit 295c6bd

Please sign in to comment.