OSCR

Interactions across hemispheres in prefrontal cortex reflect global cognitive processing.

Code ↔ Paper

20 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 20 matches · 5 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Behavioral task ↔ ex/ex_SimpleTaskDemo.m, lines 1–25 · score 0.80 · memory guided saccade, go cue, target flash, saccade initiation, radius, locations
  2. [2] § Methods › Behavioral task ↔ ex/ex_SimpleJoystickTaskDemo.m, lines 112–259 · score 0.74 · target flash, saccade initiation, target location, go, reward, cue
  3. [3] § Methods › Fitting pCCA-FA and pCCA to simulated data ↔ pcca_fa_mdl.py, lines 11–54 · score 0.74 · EM algorithm, cross validated, area dimensionalities, pCCA FA model, latent variables, dshared
  4. [4] § Methods › Fitting pCCA-FA ↔ pcca_fa_mdl.py, lines 519–601 · score 0.74 · fold cross validation, area dimensionalities, pCCA FA model, likelihood, integer, d1
  5. [5] § Methods › Fitting pCCA-FA ↔ pcca_fa_mdl.py, lines 11–54 · score 0.70 · expectation maximization, EM algorithm, pCCA FA model, orthogonal, fit, matrix
  6. [6] § Methods › Preprocessing of neural activity ↔ ex_bci/utils/getChannelsKeepWithDat.m, the whole file · a weak match · score 0.67 · low firing rates, Fano factor, coincident, electrodes, binned, spikes
  7. [7] § Methods › Computing signal correlation ↔ main_analyses/compute_rsc.py, lines 29–78 · score 0.66 · signal correlation, area rsc, permutation, opposite, tuning, raw
  8. [8] § Methods › Separation of slow and fast components ↔ plot_figureS1.ipynb, lines 18–67 · score 0.66 · timescale component, remove slow timescale, slow component, fast, correlation, activity
  9. [9] § Methods › Preprocessing of neural activity ↔ ex_bci/utils/preprocessDatWithNoStimTrial.m, the whole file · a weak match · score 0.66 · low firing rates, Fano factor, coincident, Preprocessing, binned, spikes
  10. [10] § Methods › Simulated data generation ↔ sim_pcca_fa.py, lines 25–167 · score 0.65 · area co fluctuation, model parameters, orthogonalized, simulated, vector, matrix
  11. [11] § Methods › Fitting pCCA-FA and pCCA to simulated data ↔ pcca_fa_mdl.py, lines 519–601 · score 0.64 · EM algorithm, cross validated, area dimensionality, pCCA FA, fit, model
  12. [12] § Methods › Preprocessing of pupil diameter ↔ main_analyses/compile_pupil_data.m, the whole file · a weak match · score 0.62 · delay period, subtracted, smoothed, pupil, noise, window
  13. [13] § Methods › Probabilistic canonical correlation analysis - factor analysis (pCCA-FA) ↔ sim_pcca_fa.py, lines 25–167 · score 0.61 · independent variance, area latent variables, area component, n1, d1, d2
  14. [14] § Methods › Preprocessing of pupil diameter ↔ ex_control/samp.m, the whole file · a weak match · score 0.59 · EyeLink, Eye position, smoothed, tracking, pupil, window
  15. [15] § Results › Across-area latent variables were related to pupil diameter ↔ main_analyses/compile_pupil_data.m, the whole file · a weak match · score 0.56 · event related, target onset, delay period, pupil
  16. [16] § Methods › Pupil diameter regression ↔ main_analyses/compute_pupil_pred_1d.py, lines 47–98 · score 0.54 · linear regression, trial pupil, r2, pCCA FA, latent, model
  17. [17] § Methods › Pupil diameter regression ↔ main_analyses/compute_pupil_pred.py, lines 47–96 · score 0.54 · linear regression, trial pupil, r2, pCCA FA, latent, model
  18. [18] § Results › pCCA-FA partitions across- and within-area shared variance ↔ plot_figure3.ipynb, lines 173–255 · score 0.52 · ground truth sv, model validation, dshared, pCCA FA, dimensionality, variance
  19. [19] § Methods › Probabilistic canonical correlation analysis - factor analysis (pCCA-FA) ↔ ex_bci/oldBci/FA_distanceBCIsystemBlock/calibrateDistanceBCIFA.m, lines 101–207 · score 0.52 · low rank, covariance matrices, correlation, FA, neurons
  20. [20] § Methods › Probabilistic canonical correlation analysis - factor analysis (pCCA-FA) ↔ ex_bci/oldBci/FA_distanceBCIsystemBlock2BCI/calibrateDistanceBCIFA.m, lines 79–195 · score 0.52 · low rank, covariance matrices, correlation, FA, neurons

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 984 lines · 45 KB · Apache-2.0 · 4 matches

  1. import numpy as np
  2. import cca.prob_cca as pcca
  3. import fa.factor_analysis as fa_mdl
  4. import scipy.linalg as slin
  5. import sklearn.model_selection as ms
  6. from joblib import Parallel,delayed
  7. from functools import partial
  8. from psutil import cpu_count
  9. from tqdm import tqdm
  11. class pcca_fa:
  12. '''
  13. pCCA-FA is a dimensionality reduction framework that combines probabilistic canonical correlation analysis (pCCA)
  14. and factor analysis (FA) to model across- and within- dataset interactions.
  16. This class implements the pCCA-FA model, stores parameters, and contains methods for fitting the model to data and computing model metrics.
  18. Methods
  19. -------
  20. train()
  21. Fit a pCCA-FA model to data using expectation-maximization (EM) algorithm.
  22. get_loading_matrices()
  23. Get across- and within-area loading matrices of the fit model.
  24. get_canonical_directions()
  25. Get canonical directions from the parameters of the fit model, as in canonical correlation analysis (CCA).
  26. get_correlative_modes()
  27. Transforms across-area loading matrices to their correlative modes.
  28. get_params()
  29. Get parameters of the fit model.
  30. set_params()
  31. Set parameters of the model.
  32. estep()
  33. Compute expectation of the posterior, according to the E-step of the EM algorithm.
  34. orthogonalize()
  35. Orthogonalize across- and within-area loading matrices using singular value decomposition.
  36. orthogonalize_latents()
  37. Orthogonalize latent variables (posterior means) using singular value decomposition.
  38. crossvalidate()
  39. Perform k-fold cross-validation to select hyperparameters (optimal across- and within-area dimensionality),
  40. then fit a pCCA-FA model with the selected hyperparameters.
  42. Model metric methods
  43. -------
  44. compute_load_sim()
  45. Compute loading similarity in each across- and within-area loading matrix.
  46. compute_dshared()
  47. Compute shared dimensionality (d_shared) in each across- and within-area loading matrix.
  48. compute_part_ratio()
  49. Compute part ratio in each across- and within-area loading matrix.
  50. compute_psv()
  51. Compute percentage of shared variance (%sv) in each across- and within-area loading matrix.
  52. compute_metrics()
  53. Wrapper to compute loading similarity, d_shared, part ratio, %sv, and canonical correlations.
  54. '''
  56. def __init__(self,min_var=0.01):
  57. '''
  58. Initialize pCCA-FA model class.
  60. Parameters:
  61. min_var (float): Used to set the variance floor, to prevent numerical underflow.
  62. '''
  63. self.params = []
  64. self.min_var = min_var
  66. def train(self,X_1,X_2,d,d1,d2,tol=1e-6,max_iter=int(1e6),verbose=False,rand_seed=None,warmstart=True,X_1_early_stop=None,X_2_early_stop=None,start_params=None,parallelize=True):
  67. '''
  68. Fit a pCCA-FA model to data using expectation-maximization (EM) algorithm.
  70. Parameters:
  71. X_1 (array): Array of size N (trials) x n1 (neurons), spike counts in area 1
  72. X_2 (array): Array of size N (trials) x n2 (neurons), spike counts in area 2
  73. d (int): Across-area dimensionality
  74. d1 (int): Within-area dimensionality for area 1
  75. d2 (int): Within-area dimensionality for area 2
  76. tol (float): Tolerance for convergence of the EM algorithm
  77. max_iter (int): Maximum number of iterations of the EM algorithm
  78. verbose (bool): Flag to print out updates during training
  79. rand_seed (int): Seed for random number generator, provide to ensure reproducibility
  80. warmstart (bool): Whether to initialize starting parameters of EM algorithm using pCCA and FA
  81. X_1_early_stop (array): Array of size N (trials) x n1 (neurons), test spike counts in area 1
  82. X_2_early_stop (array): Array of size N (trials) x n2 (neurons), test spike counts in area 2
  83. start_params (dict): Dictionary containing pCCA-FA model parameters to initialize EM algorithm
  84. use_process (bool): Whether to run training in a separate process for better performance
  86. Returns:
  87. LL (array): Training data log likelihood at each iteration of EM algorithm
  88. testLL (array): If using early_stop test data, contains test data log likelihood at each iteration of EM algorithm. Empty array otherwise.
  89. '''
  91. if parallelize:
  92. # Define function to run in parallel
  93. def _train_wrapper(X_1, X_2, d, d1, d2, min_var, tol, max_iter, verbose, rand_seed,
  94. warmstart, X_1_early_stop, X_2_early_stop, start_params):
  95. model = pcca_fa(min_var=min_var)
  96. LL, testLL = model.train(X_1, X_2, d, d1, d2,
  97. tol=tol, max_iter=max_iter,
  98. verbose=verbose, rand_seed=rand_seed,
  99. warmstart=warmstart,
  100. X_1_early_stop=X_1_early_stop,
  101. X_2_early_stop=X_2_early_stop,
  102. start_params=start_params,
  103. parallelize=False)
  104. return LL, testLL, model.get_params()
  106. # Run in parallel with all arguments explicitly passed
  107. result = Parallel(n_jobs=cpu_count(logical=False), backend='loky')([
  108. delayed(_train_wrapper)(X_1, X_2, d, d1, d2, self.min_var, tol, max_iter,
  109. verbose, rand_seed, warmstart, X_1_early_stop,
  110. X_2_early_stop, start_params)]
  111. )[0]
  113. LL, testLL, self.params = result
  114. return LL, testLL
  116. # Regular in-process training
  117. # set random seed
  118. if not(rand_seed is None):
  119. np.random.seed(rand_seed)
  121. early_stop = not(X_1_early_stop is None) and not(X_2_early_stop is None)
  123. # set some useful parameters
  124. N,n1 = X_1.shape
  125. _,n2 = X_2.shape
  126. mu_x1,mu_x2 = X_1.mean(axis=0),X_2.mean(axis=0)
  127. X_1c,X_2c = (X_1-mu_x1), (X_2-mu_x2)
  128. X_total = np.concatenate((X_1c,X_2c),axis=1)
  129. covX_1 = 1/N * (X_1c.T).dot(X_1c)
  130. covX_2 = 1/N * (X_2c.T).dot(X_2c)
  131. sampleCov = 1/N * (X_total.T).dot(X_total)
  132. var_floor = self.min_var*np.diag(sampleCov)
  133. Iz = np.identity(d+d1+d2)
  134. const = (n1+n2)*np.log(2*np.pi)
  136. if early_stop:
  137. X_1c_test,X_2c_test = X_1_early_stop - mu_x1, X_2_early_stop - mu_x2
  138. X_total_test = np.concatenate((X_1c_test,X_2c_test),axis=1)
  139. cov_test = (1/N)*(X_total_test.T.dot(X_total_test))
  141. # check that covariance is full rank
  142. if np.linalg.matrix_rank(sampleCov)==(n1+n2):
  143. x1_scale = np.exp(2/n1*np.sum(np.log(np.diag(slin.cholesky(covX_1)))))
  144. x2_scale = np.exp(2/n2*np.sum(np.log(np.diag(slin.cholesky(covX_2)))))
  145. else:
  146. raise np.linalg.LinAlgError(f'Covariance matrix is low rank ({np.linalg.matrix_rank(sampleCov):d}, should be {n1+n2:d})')
  148. # initialize parameters
  149. if warmstart:
  150. # get across-area loading matrices from pCCA and within-area loading matrices from FA
  151. tmp = pcca.prob_cca()
  152. tmp.train_maxLL(X_1,X_2,d)
  153. W_1 = tmp.get_params()['W_x']
  154. W_2 = tmp.get_params()['W_y']
  155. tmp = fa_mdl.factor_analysis()
  156. tmp.train(X_1,d1,rand_seed=rand_seed)
  157. L_1 = tmp.get_params()['L']
  158. tmp = fa_mdl.factor_analysis()
  159. tmp.train(X_2,d2,rand_seed=rand_seed)
  160. L_2 = tmp.get_params()['L']
  161. Psi = np.diag(sampleCov)
  162. elif not(start_params is None):
  163. # allow for specifying parameter initialization
  164. W_1 = start_params['W_1']
  165. W_2 = start_params['W_2']
  166. L_1 = start_params['L_1']
  167. L_2 = start_params['L_2']
  168. Psi = np.abs(np.append(start_params['psi_1'], start_params['psi_2']))
  169. else:
  170. if d > 0:
  171. W_1 = np.random.randn(n1,d) * np.sqrt(x1_scale/d)
  172. W_2 = np.random.randn(n2,d) * np.sqrt(x2_scale/d)
  173. else:
  174. W_1 = np.random.randn(n1,d)
  175. W_2 = np.random.randn(n2,d)
  176. if d1 > 0:
  177. L_1 = np.random.randn(n1,d1) * np.sqrt(x1_scale/d1)
  178. else:
  179. L_1 = np.random.randn(n1,d1)
  180. if d2 > 0:
  181. L_2 = np.random.randn(n2,d2) * np.sqrt(x2_scale/d2)
  182. else:
  183. L_2 = np.random.randn(n2,d2)
  184. Psi = np.diag(sampleCov)
  186. # define L_total - joint loading matrix
  187. L_top = np.concatenate((W_1,L_1,np.zeros((n1,d2))),axis=1)
  188. L_bottom = np.concatenate((W_2,np.zeros((n2,d1)),L_2),axis=1)
  189. L_total = np.concatenate((L_top,L_bottom),axis=0)
  191. L_mask = np.ones(L_total.shape)
  192. L_mask[:n1,(d+d1):] = np.zeros((n1,d2))
  193. L_mask[n1:,d:(d+d1)] = np.zeros((n2,d1))
  195. # EM algorithm
  196. LL = []
  197. testLL = []
  198. for i in range(max_iter):
  199. # E-step: set q(z) = p(z,zx,zy|x,y)
  200. iPsi = np.diag(1/Psi)
  201. iPsiL = iPsi.dot(L_total)
  202. if d==0 and d1==0 and d2==0:
  203. iSig = iPsi
  204. else:
  205. iSig = iPsi - iPsiL.dot(slin.inv(Iz+(L_total.T).dot(iPsiL))).dot(iPsiL.T)
  206. iSigL = iSig.dot(L_total)
  207. cov_iSigL = sampleCov.dot(iSigL)
  208. E_zz = Iz - (L_total.T).dot(iSigL) + (iSigL.T).dot(cov_iSigL)
  210. # compute log likelihood
  211. logDet = 2*np.sum(np.log(np.diag(slin.cholesky(iSig))))
  212. curr_LL = -N/2 * (const - logDet + np.trace(iSig.dot(sampleCov)))
  213. LL.append(curr_LL)
  214. if early_stop:
  215. curr_testLL = -N/2 * (const - logDet + np.trace(iSig.dot(cov_test)))
  216. testLL.append(curr_testLL)
  217. if verbose:
  218. print('EM iteration ',i,', LL={:.2f}'.format(curr_LL))
  220. # check for convergence (training LL increases by less than tol, or testLL decreases)
  221. if i>1:
  222. if (LL[-1]-LL[-2])<tol or (early_stop and testLL[-1]<testLL[-2]):
  223. break
  225. # M-step: compute new L and Psi
  226. if not(d==0 and d1==0 and d2==0):
  227. L_total = cov_iSigL.dot(slin.inv(E_zz))
  228. L_total = L_total * L_mask
  229. Psi = np.diag(sampleCov) - np.diag(cov_iSigL.dot(L_total.T))
  230. Psi = np.maximum(Psi,var_floor)
  232. # get final parameters after convergence or max_iter
  233. W_1, W_2 = L_total[:n1,:d], L_total[n1:,:d]
  234. L_1, L_2 = L_total[:n1,d:(d+d1)], L_total[n1:,(d+d1):]
  235. psi_1, psi_2 = Psi[:n1], Psi[n1:]
  237. # create parameter dict
  238. self.params = {
  239. 'mu_x1':mu_x1,'mu_x2':mu_x2, # estimated mean per neuron
  240. 'L_total':L_total, # maximum likelihood estimated matrix
  241. 'W_1':W_1,'W_2':W_2, # across-area loading matrices
  242. 'L_1':L_1,'L_2':L_2, # within-area loading matrices
  243. 'psi_1':psi_1,'psi_2':psi_2, # private variance per neuron
  244. 'd':d,'d1':d1,'d2':d2, # selected dimensionalities
  245. }
  247. return np.array(LL), np.array(testLL)
  249. def get_loading_matrices(self):
  250. '''
  251. Get across- and within-area loading matrices of the fit model.
  253. Returns:
  254. W_1 (array): Array of size n1 (neurons) x d (latents) containing the loadings for across-area latent variables onto neurons in area 1
  255. W_2 (array): Array of size n2 (neurons) x d (latents) containing the loadings for across-area latent variables onto neurons in area 2
  256. L_1 (array): Array of size n1 (neurons) x d1 (latents) containing the loadings for within-area latent variables onto neurons in area 1
  257. L_2 (array): Array of size n2 (neurons) x d2 (latents) containing the loadings for within-area latent variables onto neurons in area 2
  258. '''
  260. n1 = len(self.params['mu_x1'])
  261. d, d1 = self.params['d'], self.params['d1']
  262. L_total = self.params['L_total']
  264. # get final parameters
  265. W_1, W_2 = L_total[:n1,:d], L_total[n1:,:d]
  266. L_1, L_2 = L_total[:n1,d:(d+d1)], L_total[n1:,(d+d1):]
  268. return W_1, W_2, L_1, L_2
  270. def get_canonical_directions(self):
  271. '''
  272. Get canonical directions from the parameters of the fit model, as in canonical correlation analysis (CCA).
  274. Returns:
  275. canonical_dirs_x (array): Array of size n1 (neurons) x d (latents) whose columns contain the canonical directions for area 1
  276. canonical_dirs_y (array): Array of size n2 (neurons) x d (latents) whose columns contain the canonical directions for area 2
  277. rho (array): Array of size d (latents) x 1 containing the corresponding canonical correlations
  278. '''
  280. W_1, W_2, L_1, L_2 = self.get_loading_matrices()
  281. psi_1, psi_2 = self.params['psi_1'], self.params['psi_2']
  282. d = self.params['d']
  284. # compute canonical correlations
  285. est_covX_1 = W_1.dot(W_1.T) + L_1.dot(L_1.T) + np.diag(psi_1)
  286. est_covX_2 = W_2.dot(W_2.T) + L_2.dot(L_2.T) + np.diag(psi_2)
  287. est_covX_1X_2 = W_1.dot(W_2.T)
  288. inv_sqrt_covX_1 = slin.inv(slin.sqrtm(est_covX_1))
  289. inv_sqrt_covX_2 = slin.inv(slin.sqrtm(est_covX_2))
  290. K = inv_sqrt_covX_1.dot(est_covX_1X_2).dot(inv_sqrt_covX_2)
  291. u,s,vt = slin.svd(K)
  292. rho = s[0:d]
  294. canonical_dirs_x = slin.inv(slin.sqrtm(est_covX_1)) @ u[:,:d]
  295. canonical_dirs_y = slin.inv(slin.sqrtm(est_covX_2)) @ vt[:d,:].T
  297. return (canonical_dirs_x, canonical_dirs_y), rho
  299. def get_correlative_modes(self):
  300. '''
  301. Transforms across-area loading matrices to their correlative modes.
  303. Follows equations in Bach & Jordan, 2005.
  305. Returns:
  306. CorrModes_x (array): Array of size n1 (neurons) x d (latents) whose columns contain the correlative modes for area 1
  307. CorrModes_y (array): Array of size n2 (neurons) x d (latents) whose columns contain the correlative modes for area 2
  308. '''
  310. W_1, W_2, L_1, L_2 = self.get_loading_matrices()
  311. psi_1, psi_2 = self.params['psi_1'], self.params['psi_2']
  312. d = self.params['d']
  314. # compute canonical correlations
  315. est_covX_1 = W_1.dot(W_1.T) + L_1.dot(L_1.T) + np.diag(psi_1)
  316. est_covX_2 = W_2.dot(W_2.T) + L_2.dot(L_2.T) + np.diag(psi_2)
  317. est_covX_1X_2 = W_1.dot(W_2.T)
  318. inv_sqrt_covX_1 = slin.inv(slin.sqrtm(est_covX_1))
  319. inv_sqrt_covX_2 = slin.inv(slin.sqrtm(est_covX_2))
  320. K = inv_sqrt_covX_1.dot(est_covX_1X_2).dot(inv_sqrt_covX_2)
  321. u,s,vt = slin.svd(K)
  322. rho = s[0:d]
  324. # order W_1, W_2 by canon corrs
  325. pd = np.diag(np.sqrt(rho))
  326. CorrModes_x = slin.sqrtm(est_covX_1).dot(u[:,0:d]).dot(pd)
  327. CorrModes_y = slin.sqrtm(est_covX_2).dot(vt[0:d,:].T).dot(pd)
  329. return CorrModes_x, CorrModes_y
  331. def get_params(self):
  332. '''
  333. Get parameters of the fit model.
  335. Returns:
  336. params (dict): Dictionary containing each parameter of the pCCA-FA model
  337. '''
  338. return self.params
  340. def set_params(self,params):
  341. '''
  342. Set parameters of the model.
  344. Parameters:
  345. params (dict): Dictionary containing each parameter of the pCCA-FA model
  346. '''
  347. self.params = params
  349. def estep(self,X_1,X_2):
  350. '''
  351. Compute expectation of the posterior, according to the E-step of the EM algorithm.
  353. Parameters:
  354. X_1 (array): Array of size N (trials) x n1 (neurons), spike counts in area 1
  355. X_2 (array): Array of size N (trials) x n2 (neurons), spike counts in area 2
  357. Returns:
  358. z (dict): Dictionary containing the mean and covariance of the posterior
  359. LL (float): Log likelihood of the provided spike counts X_1 and X_2 under the fit model parameters
  360. '''
  362. N,n1 = X_1.shape
  363. _,n2 = X_2.shape
  365. d,d1,d2 = self.params['d'],self.params['d1'],self.params['d2']
  367. # get model parameters
  368. mu_x1,mu_x2 = self.params['mu_x1'],self.params['mu_x2']
  369. L_total = self.params['L_total']
  370. psi_1 = self.params['psi_1']
  371. psi_2 = self.params['psi_2']
  372. Psi = np.diag(np.concatenate((psi_1,psi_2)))
  374. # center data and compute covariances
  375. X_1c = X_1-mu_x1
  376. X_2c = X_2-mu_x2
  377. X_total = np.concatenate((X_1c,X_2c),axis=1)
  378. sampleCov = 1/N * (X_total.T).dot(X_total)
  380. # compute z
  381. Iz = np.identity(d+d1+d2)
  382. C = L_total.dot(L_total.T) + Psi
  383. invC = slin.inv(C)
  384. z_mu = X_total.dot(invC).dot(L_total)
  385. z_cov = np.diag(np.diag(Iz - (L_total.T).dot(invC).dot(L_total)))
  387. # compute LL
  388. const = (n1+n2)*np.log(2*np.pi)
  389. logDet = 2*np.sum(np.log(np.diag(slin.cholesky(C))))
  390. LL = -N/2 * (const + logDet + np.trace(invC.dot(sampleCov)))
  392. # return posterior and LL
  393. z = {
  394. 'z_mu':z_mu[:,:d],
  395. 'z_cov':z_cov[:d,:d],
  396. 'zx1_mu':z_mu[:,d:(d+d1)],
  397. 'zx1_cov':z_cov[d:(d+d1),d:(d+d1)],
  398. 'zx2_mu':z_mu[:,(d+d1):],
  399. 'zx2_cov':z_cov[(d+d1):,(d+d1):],
  400. }
  401. return z, LL
  403. def orthogonalize(self,across_mode='paired'):
  404. '''
  405. Orthogonalize across- and within-area loading matrices using singular value decomposition.
  407. Note: this also transforms loading matrices to be in covariant modes (as opposed to correlative modes)
  409. Parameters:
  410. across_mode (str): Parameter to indicate whether to orthogonalize the across-area loading matrices jointly ('paired') or individually in each area ('unpaired')
  412. Returns:
  413. W_1_norm (array): Array of size n1 (neurons) x d (latents) containing orthogonal columns with the loadings for across-area latent variables onto neurons in area 1
  414. W_2_norm (array): Array of size n2 (neurons) x d (latents) containing orthogonal columns the loadings for across-area latent variables onto neurons in area 2
  415. L_1_norm (array): Array of size n1 (neurons) x d1 (latents) containing orthogonal columns the loadings for within-area latent variables onto neurons in area 1
  416. L_2_norm (array): Array of size n2 (neurons) x d2 (latents) containing orthogonal columns the loadings for within-area latent variables onto neurons in area 2
  417. '''
  419. n1 = len(self.params['mu_x1'])
  420. W_1, W_2, L_1, L_2 = self.get_loading_matrices() # output from maximum likelihood estimation
  422. # within-area loading matrices
  423. u,s,_ = slin.svd(L_1,full_matrices=False)
  424. L_1_norm = u @ np.diag(s)
  425. u,s,_ = slin.svd(L_2,full_matrices=False)
  426. L_2_norm = u @ np.diag(s)
  428. # across-area loading matrices
  429. if across_mode == 'paired':
  430. W_total = np.concatenate((W_1,W_2),axis=0)
  431. u,s,_ = slin.svd(W_total,full_matrices=False)
  432. W_1_norm = u[:n1,:] @ np.diag(s)
  433. W_2_norm = u[n1:,:] @ np.diag(s)
  434. elif across_mode == 'unpaired':
  435. u,s,_ = slin.svd(W_1,full_matrices=False)
  436. W_1_norm = u @ np.diag(s)
  437. u,s,_ = slin.svd(W_2,full_matrices=False)
  438. W_2_norm = u @ np.diag(s)
  439. else:
  440. raise ValueError('across-mode must be "paired" or "unpaired"')
  442. return W_1_norm, W_2_norm, L_1_norm, L_2_norm
  444. def orthogonalize_latents(self,zx_mu,zy_mu,do_across=False,z_mu=None,across_mode='paired'):
  445. '''
  446. Orthogonalize latent variables (posterior means) using singular value decomposition.
  448. Parameters:
  449. zx_mu (array): Array of size N (trials) x d1 (latents) containing the within-area latent variables or posterior mean in area 1
  450. zy_mu (array): Array of size N (trials) x d2 (latents) containing the within-area latent variables or posterior mean in area 2
  451. do_across (bool): Whether to orthogonalize the across-area latent variables (True) or not (False)
  452. z_mu (array): Array of size N (trials) x d (latents) containing the across-area latent variables or posterior mean. Only used if do_across is True
  453. across_mode (str): Parameter to indicate whether to orthogonalize the across-area latent variables jointly ('paired') or individually in each area ('unpaired'). Only used if do_across is True
  455. Returns:
  456. z_orth (dict): Dictionary containing the orthogonalized latent variables
  457. W_orth (dict): Dictionary containing the orthogonalized loading matrices
  458. '''
  460. W_1, W_2, L_1, L_2 = self.get_loading_matrices() # output from maximum likelihood estimation
  461. n1 = L_1.shape[0]
  463. # orthogonalize zx
  464. u,s,vt = slin.svd(L_1,full_matrices=False)
  465. L_1_orth = u
  466. TT = np.diag(s).dot(vt)
  467. z_x1 = (TT.dot(zx_mu.T)).T
  469. # orthogonalize zy
  470. u,s,vt = slin.svd(L_2,full_matrices=False)
  471. L_2_orth = u
  472. TT = np.diag(s).dot(vt)
  473. z_x2 = (TT.dot(zy_mu.T)).T
  475. # orthogonalize across-area
  476. across_z_orth = {}
  477. if do_across:
  478. if across_mode == 'paired':
  479. # orthogonalize across-area latents using both area's loading matrix
  480. W_total = np.concatenate((W_1,W_2),axis=0)
  481. u,s,vt = slin.svd(W_total,full_matrices=False)
  482. W_1_orth = u[:n1,:]
  483. W_2_orth = u[n1:,:]
  484. TT = np.diag(s).dot(vt)
  485. z = (TT.dot(z_mu.T)).T
  486. across_z_orth['area1'] = z
  487. across_z_orth['area2'] = z
  488. elif across_mode == 'unpaired':
  489. # orthogonalize across-area latents using each area's loading matrix
  490. u,s,vt = slin.svd(W_1,full_matrices=False)
  491. W_1_orth = u
  492. TT = np.diag(s).dot(vt)
  493. z = (TT.dot(z_mu.T)).T
  494. across_z_orth['area1'] = z
  496. u,s,vt = slin.svd(W_2,full_matrices=False)
  497. W_2_orth = u
  498. TT = np.diag(s).dot(vt)
  499. z = (TT.dot(z_mu.T)).T
  500. across_z_orth['area2'] = z
  501. else:
  502. raise ValueError('across-mode must be "paired" or "unpaired"')
  504. # return z_orth, W_orth
  505. z_orth = {
  506. 'z':across_z_orth, # across area latent variables, empty if do_across is False
  507. 'z1':z_x1, # within-area latent variables for area 1
  508. 'z2':z_x2 # within-area latent variables for area 2
  509. }
  510. W_orth = {
  511. 'W_1':W_1_orth, # across-area loading matrix for area 1
  512. 'W_2':W_2_orth, # across-area loading matrix for area 2
  513. 'L_1':L_1_orth, # within-area loading matrix for area 1
  514. 'L_2':L_2_orth # within-area loading matrix for area 2
  515. }
  517. return z_orth, W_orth
  519. def crossvalidate(self,X_1,X_2,d_list=np.linspace(0,8,9),d1_list=np.linspace(0,8,9),d2_list=np.linspace(0,8,9),n_folds=10,verbose=True,max_iter=int(1e6),tol=1e-6,warmstart=True,rand_seed=None,parallelize=False,early_stop=False):
  520. '''
  521. Perform k-fold cross-validation to select hyperparameters (optimal across- and within-area dimensionality), then fit a pCCA-FA model with the selected hyperparameters.
  523. Parameters:
  524. X_1 (array): Array of size N (trials) x n1 (neurons), spike counts in area 1
  525. X_2 (array): Array of size N (trials) x n2 (neurons), spike counts in area 2
  526. d_list (array): 1-dimensional array containing the across-area dimensionalities to test
  527. d1_list (array): 1-dimensional array containing the within-area dimensionalities to test for area 1
  528. d2_list (array): 1-dimensional array containing the within-area dimensionalities to test for area 2
  529. n_folds (int): The number of folds (k) for cross-validation
  530. verbose (bool): Flag to print out updates during training
  531. max_iter (int): Maximum number of iterations of the EM algorithm
  532. tol (float): Tolerance for convergence of the EM algorithm
  533. warmstart (bool): Whether to initialize starting parameters of EM algorithm using pCCA and FA
  534. rand_seed (int): Seed for random number generator, provide to ensure reproducibility
  535. parallelize (bool): Whether to parallelize cross-validation folds (True) or not (False), to reduce run time
  536. early_stop (bool): Whether to use early_stop (True) or not (False) on the testing data of each cross-validation fold
  538. Returns:
  539. results (dict): Dictionary containing the lists of tested dimensionalities and their corresponding cross-validated data log likelihood and prediction errors,
  540. as well as the selected dimensionalities and its corresponding cross-validated log likelihood
  541. '''
  543. # set random seed
  544. if not(rand_seed is None):
  545. np.random.seed(rand_seed)
  547. # make sure z dims are integers
  548. d_list,d1_list,d2_list = np.meshgrid(d_list.astype(int),d1_list.astype(int),d2_list.astype(int))
  549. d_list = np.matrix.flatten(d_list)
  550. d1_list = np.matrix.flatten(d1_list)
  551. d2_list = np.matrix.flatten(d2_list)
  552. results = {'d_list':d_list,'d1_list':d1_list,'d2_list':d2_list}
  554. # create k-fold iterator
  555. if verbose:
  556. print('Crossvalidating pCCA-FA model to choose # of dims...')
  557. cv_kfold = ms.KFold(n_splits=n_folds,shuffle=True,random_state=rand_seed)
  559. # iterate through train/test splits
  560. i = 0
  561. LLs,PEs = np.zeros([n_folds,len(d_list)]),np.zeros([n_folds,len(d_list)])
  562. for train_idx,test_idx in cv_kfold.split(X_1):
  563. if verbose:
  564. print(' Fold ',i+1,' of ',n_folds,'...')
  565. X_1_train,X_1_test = X_1[train_idx], X_1[test_idx]
  566. X_2_train,X_2_test = X_2[train_idx], X_2[test_idx]
  568. # iterate through each d, provide training and testing trials to the helper function
  569. func = partial(self._cv_helper,X_1train=X_1_train,X_2train=X_2_train,X_1test=X_1_test,X_2test=X_2_test,\
  570. rand_seed=rand_seed,max_iter=max_iter,tol=tol,warmstart=warmstart,early_stop=early_stop)
  571. if parallelize:
  572. tmp = Parallel(n_jobs=cpu_count(logical=False),backend='loky')\
  573. (delayed(func)(d_list[j],d1_list[j],d2_list[j]) for j in range(len(d_list)))
  574. LLs[i,:] = [val[0] for val in tmp]
  575. PEs[i,:] = [val[1] for val in tmp]
  576. else:
  577. for j in tqdm(range(len(d_list))):
  578. tmp = func(d_list[j],d1_list[j],d2_list[j])
  579. LLs[i,j],PEs[i,j] = tmp[0],tmp[1]
  581. i = i+1
  583. sum_LLs = LLs.sum(axis=0)
  584. sum_SEs = PEs.sum(axis=0)
  585. results['LLs'] = sum_LLs
  586. results['PEs'] = sum_SEs
  588. # find the best # of z dimensions and train final pCCA-FA model
  589. max_idx = np.argmax(sum_LLs)
  590. d,d1,d2 = d_list[max_idx],d1_list[max_idx],d2_list[max_idx]
  591. results['d']=d
  592. results['d1']=d1
  593. results['d2']=d2
  594. results['final_LL'] = sum_LLs[max_idx]
  595. self.train(X_1,X_2,d,d1,d2) # sets params of the final model
  597. self.compute_cv_canonical_corrs(X_1,X_2,n_folds=n_folds,verbose=verbose,max_iter=max_iter,tol=tol,warmstart=warmstart,rand_seed=rand_seed)
  599. self.cv_results = results
  601. return results
  603. def _cv_helper(self,d,d1,d2,X_1train,X_2train,X_1test,X_2test,rand_seed=None,max_iter=int(1e5),tol=1e-6,warmstart=True,early_stop=False):
  604. '''
  605. Helper function for crossvalidate().
  607. Runs one train-test split and computes the log-likelihood and prediction error on the testing data.
  609. Parameters:
  610. d (int): Across-area dimensionality
  611. d1 (int): Within-area dimensionality for area 1
  612. d2 (int): Across-area dimensionality for area 2
  613. X_1train (array): Array of size Ntrain (trials) x n1 (neurons), training spike counts in area 1
  614. X_2train (array): Array of size Ntrain (trials) x n2 (neurons), training spike counts in area 2
  615. X_1test (array): Array of size Ntest (trials) x n1 (neurons), testing spike counts in area 1
  616. X_2test (array): Array of size Ntest (trials) x n2 (neurons), testing spike counts in area 2
  617. rand_seed (int): Seed for random number generator, provide to ensure reproducibility
  618. max_iter (int): Maximum number of iterations of the EM algorithm
  619. tol (float): Tolerance for convergence of the EM algorithm
  620. warmstart (bool): Whether to initialize starting parameters of EM algorithm using pCCA and FA
  621. early_stop (bool): Whether to use early_stop (True) or not (False) on the testing data
  623. Returns:
  624. LL (float): Cross-validated data log likelihood of the testing data
  625. PE (float): Prediction error of the testing data using leave-one-out prediction
  626. '''
  628. tmp = pcca_fa()
  629. if early_stop:
  630. tmp.train(X_1train,X_2train,d,d1,d2,rand_seed=rand_seed,max_iter=max_iter,tol=tol,warmstart=warmstart,X_1_early_stop=X_1test,X_2_early_stop=X_2test,parallelize=False)
  631. else:
  632. tmp.train(X_1train,X_2train,d,d1,d2,rand_seed=rand_seed,max_iter=max_iter,tol=tol,warmstart=warmstart,parallelize=False)
  633. # log-likelihood
  634. _,LL = tmp.estep(X_1test,X_2test)
  635. # prediction error
  636. X_1test_pred,X_2test_pred = tmp._leaveoneout_pred(X_1test,X_2test)
  637. PE = np.sum(np.square(X_1test_pred - X_1test)) + np.sum(np.square(X_2test_pred - X_2test))
  639. return (LL,PE)
  641. def _leaveoneout_pred(self,X_1,X_2):
  642. '''
  643. Helper function for crossvalidate().
  645. Runs leave-one-out prediction on provided data.
  647. Parameters:
  648. X_1 (array): Array of size N (trials) x n1 (neurons), spike counts in area 1
  649. X_2 (array): Array of size N (trials) x n2 (neurons), spike counts in area 2
  651. Returns:
  652. pred_x (array): Array of size N (trials) x n1 (neurons) containing the prediction errors for area 1
  653. pred_y (array): Array of size N (trials) x n2 (neurons) containing the prediction errors for area 2
  654. '''
  656. N,n1 = X_1.shape # trials x neurons
  657. n2 = X_2.shape[1]
  658. X_total = np.concatenate((X_1,X_2),axis=1)
  660. # extract model parameters
  661. Psi = np.concatenate((self.params['psi_1'],self.params['psi_2']),axis=0)
  662. mu = np.concatenate((self.params['mu_x1'],self.params['mu_x2']),axis=0)
  663. L_total = self.params['L_total']
  665. # compute covariances
  666. mdl_cov = L_total.dot(L_total.T) + np.diag(Psi)
  667. inv_cov = slin.inv(mdl_cov)
  669. # compute conditional expectations (predictions)
  670. n_total = n1+n2
  671. pred_total = np.zeros((N,n_total))
  672. for i in range(n_total):
  673. E = np.delete(np.delete(inv_cov,i,axis=0),i,axis=1)
  674. f = np.delete(inv_cov[:,i],i,axis=0)
  675. h = inv_cov[i,i]
  676. inv_term = E - (1/h)*np.outer(f,f)
  677. proj_term = np.delete(mdl_cov[i,:],i) # 1 x n_total-1
  678. mean_term = np.delete(X_total,i,axis=1) - np.delete(mu,i,axis=0).T # N x n_total-1
  679. pred = mu[i] + proj_term.dot(inv_term.dot(mean_term.T)) # 1 x N
  680. pred_total[:,i] = pred.T
  682. pred_x = pred_total[:,:n1] # predictions for neurons in area 1
  683. pred_y = pred_total[:,n1:] # predictions for neurons in area 2
  685. return pred_x,pred_y
  687. def compute_cv_canonical_corrs(self,X_1,X_2,n_folds=10,verbose=False,max_iter=int(1e5),tol=1e-6,warmstart=True,rand_seed=None):
  688. '''
  689. Get cross-validated canonical correlations from the fit model.
  691. Parameters:
  692. Parameters:
  693. X_1 (array): Array of size N (trials) x n1 (neurons), spike counts in area 1
  694. X_2 (array): Array of size N (trials) x n2 (neurons), spike counts in area 2
  695. n_folds (int): The number of folds (k) for cross-validation
  696. verbose (bool): Flag to print out updates during training
  697. max_iter (int): Maximum number of iterations of the EM algorithm
  698. tol (float): Tolerance for convergence of the EM algorithm
  699. warmstart (bool): Whether to initialize starting parameters of EM algorithm using pCCA and FA
  700. rand_seed (int): Seed for random number generator, provide to ensure reproducibility
  702. Returns:
  703. cv_rho (array): Array of size d (latents) x 1 containing the cross-validated canonical correlations
  704. '''
  706. # set random seed
  707. if not(rand_seed is None):
  708. np.random.seed(rand_seed)
  710. # check if model has been fit
  711. if not self.params:
  712. raise ValueError('Model must be fit before computing cross-validated canonical correlations. Run train() or crossvalidate() first.')
  714. # cross-validate to get cross-validated canonical correlations
  715. if verbose:
  716. print('Crossvalidating pCCA-FA model to compute canon corrs...')
  718. # set up needed parameters
  719. d,d1,d2 = self.params['d'],self.params['d1'],self.params['d2']
  720. N = X_1.shape[0]
  722. cv_kfold = ms.KFold(n_splits=n_folds,shuffle=True,random_state=rand_seed)
  723. zx1,zx2 = np.zeros((2,N,d))
  724. i=0
  725. for train_idx,test_idx in cv_kfold.split(X_1):
  726. if verbose:
  727. print(' Fold ',i+1,' of ',n_folds,'...')
  729. X_1_train,X_1_test = X_1[train_idx], X_1[test_idx]
  730. X_2_train,X_2_test = X_2[train_idx], X_2[test_idx]
  732. tmp = pcca_fa()
  733. tmp.train(X_1_train,X_2_train,d,d1,d2,rand_seed=rand_seed,max_iter=max_iter,tol=tol,warmstart=warmstart)
  734. W_1,W_2,L_1,L_2 = tmp.get_loading_matrices() # take direct EM outputs to compute E-step
  735. tmp_params = tmp.get_params()
  737. # compute pCCA E-step: E[z|x] and E[z|y]
  738. X_1c = X_1_test - tmp_params['mu_x1']
  739. Cx1 = W_1 @ W_1.T + (L_1 @ L_1.T + np.diag(tmp_params['psi_1']))
  740. invCx1 = slin.inv(Cx1)
  741. zx1_mu = X_1c.dot(invCx1).dot(W_1)
  743. X_2c = X_2_test - tmp_params['mu_x2']
  744. Cx2 = W_2 @ W_2.T + (L_2 @ L_2.T + np.diag(tmp_params['psi_2']))
  745. invCx2 = slin.inv(Cx2)
  746. zx2_mu = X_2c.dot(invCx2).dot(W_2)
  748. zx1[test_idx,:] = zx1_mu
  749. zx2[test_idx,:] = zx2_mu
  751. i+=1
  753. cv_rho = np.zeros(d)
  754. for i in range(d):
  755. tmp = np.corrcoef(zx1[:,i],zx2[:,i])
  756. cv_rho[i] = tmp[0,1]
  758. self.params['cv_rho'] = cv_rho
  759. return cv_rho
  761. def compute_load_sim(self):
  762. '''
  763. Compute loading similarity in each across- and within-area loading matrix.
  765. Returns:
  766. ls (dict): Dictionary containing the loading similarity for each across- and within-area loading matrix
  768. '''
  769. n1 = self.params['W_1'].shape[0]
  770. n2 = self.params['W_2'].shape[0]
  772. # first, orthonormalize each loading matrix
  773. W_1,_,_ = slin.svd(self.params['W_1'],full_matrices=False)
  774. W_2,_,_ = slin.svd(self.params['W_2'],full_matrices=False)
  775. L_1,_,_ = slin.svd(self.params['L_1'],full_matrices=False)
  776. L_2,_,_ = slin.svd(self.params['L_2'],full_matrices=False)
  778. # calculate loading similarity - following equation in Umakantha, Morina, Cowley, et al., 2021.
  779. ls_W_1 = 1 - n1*W_1.var(axis=0,ddof=0)
  780. ls_W_2 = 1 - n2*W_2.var(axis=0,ddof=0)
  781. ls_L_1 = 1 - n1*L_1.var(axis=0,ddof=0)
  782. ls_L_2 = 1 - n2*L_2.var(axis=0,ddof=0)
  784. ls = {
  785. 'ls_W_1':ls_W_1, # across-area loading similarity for area 1, involves W_1
  786. 'ls_W_2':ls_W_2, # across-area loading similarity for area 2, involves W_2
  787. 'ls_L_1':ls_L_1, # within-area loading similarity for area 1, involves L_1
  788. 'ls_L_2':ls_L_2 # within-area loading similarity for area 2, involves L_2
  789. }
  790. return ls
  792. def compute_dshared(self,cutoff_thresh=0.95):
  793. '''
  794. Compute shared dimensionality (d_shared) in each across- and within-area loading matrix.
  796. Parameters:
  797. cutoff_thresh (float): Cutoff percentage (0-1) of across- or within-area shared variance to explain for selecting d_shared
  799. Returns:
  800. dshared (dict): Dictionary containing the across- and within-area d_shared for each area
  801. '''
  803. W_1,W_2,L_1,L_2 = self.get_loading_matrices()
  805. # for across-area
  806. if self.params['d'] > 0:
  807. # area 1
  808. shared = W_1.dot(W_1.T)
  809. s = slin.svdvals(shared) # eigenvalues of WWT
  810. var_exp = np.cumsum(s)/np.sum(s)
  811. dims = np.where(var_exp >= (cutoff_thresh - 1e-9))[0]
  812. dshared_W_1 = dims[0]+1
  814. # area 2
  815. shared = W_2.dot(W_2.T)
  816. s = slin.svdvals(shared) # eigenvalues of WWT
  817. var_exp = np.cumsum(s)/np.sum(s)
  818. dims = np.where(var_exp >= (cutoff_thresh - 1e-9))[0]
  819. dshared_W_2 = dims[0]+1
  821. # overall
  822. W_total = np.concatenate((W_1,W_2),axis=0)
  823. shared = W_total.dot(W_total.T)
  824. s = slin.svdvals(shared) # eigenvalues of WWT
  825. var_exp = np.cumsum(s)/np.sum(s)
  826. dims = np.where(var_exp >= (cutoff_thresh - 1e-9))[0]
  827. dshared_W_total = dims[0]+1
  828. else:
  829. dshared_W_1 = 0
  830. dshared_W_2 = 0
  831. dshared_W_total = 0
  833. # for within area 1
  834. if self.params['d1'] > 0:
  835. shared = L_1.dot(L_1.T)
  836. s = slin.svdvals(shared) # eigenvalues of LLT
  837. var_exp = np.cumsum(s)/np.sum(s)
  838. dims = np.where(var_exp >= (cutoff_thresh - 1e-9))[0]
  839. dshared_L_1 = dims[0]+1
  840. else:
  841. dshared_L_1 = 0
  843. # for within area 2
  844. if self.params['d2'] > 0:
  845. shared = L_2.dot(L_2.T)
  846. s = slin.svdvals(shared) # eigenvalues of LLT
  847. var_exp = np.cumsum(s)/np.sum(s)
  848. dims = np.where(var_exp >= (cutoff_thresh - 1e-9))[0]
  849. dshared_L_2 = dims[0]+1
  850. else:
  851. dshared_L_2 = 0
  853. dshared = {
  854. 'dshared_W_1':dshared_W_1, # d_shared for across-area shared variance in area 1, involves W_1
  855. 'dshared_W_2':dshared_W_2, # d_shared for across-area shared variance in area 2, involves W_2
  856. 'dshared_L_1':dshared_L_1, # d_shared for within-area shared variance in area 1, involves L_1
  857. 'dshared_L_2':dshared_L_2, # d_shared for within-area shared variance in area 2, involves L_2
  858. 'dshared_W_total':dshared_W_total # d_shared for across-area shared variance jointly for area 1 and 2, involves W_1 and W_2
  859. }
  861. return dshared
  863. def compute_part_ratio(self):
  864. '''
  865. Compute part ratio in each across- and within-area loading matrix.
  867. Returns:
  868. pr (dict): Dictionary containing the part ratio for each across- and within-area loading matrix
  869. '''
  871. W_1,W_2,L_1,L_2 = self.get_loading_matrices()
  873. # for across-area
  874. shared = W_1.dot(W_1.T)
  875. s = slin.svdvals(shared)
  876. pr_W_1 = np.square(s.sum()) / np.square(s).sum()
  878. shared = W_2.dot(W_2.T)
  879. s = slin.svdvals(shared)
  880. pr_W_2 = np.square(s.sum()) / np.square(s).sum()
  882. # overall
  883. W_total = np.concatenate((W_1,W_2),axis=0)
  884. shared = W_total.dot(W_total.T)
  885. s = slin.svdvals(shared)
  886. pr_W_total = np.square(s.sum()) / np.square(s).sum()
  888. # for within area 1
  889. shared = L_1.dot(L_1.T)
  890. s = slin.svdvals(shared)
  891. pr_L_1 = np.square(s.sum()) / np.square(s).sum()
  893. # for within area 2
  894. shared = L_2.dot(L_2.T)
  895. s = slin.svdvals(shared)
  896. pr_L_2 = np.square(s.sum()) / np.square(s).sum()
  898. pr = {
  899. 'pr_W_1':pr_W_1, # part ratio for across-area loading matrix in area 1, involves W_1
  900. 'pr_W_2':pr_W_2, # part ratio for across-area loading matrix in area 2, involves W_2
  901. 'pr_L_1':pr_L_1, # part ratio for within-area loading matrix in area 1, involves L_1
  902. 'pr_L_2':pr_L_2, # part ratio for within-area loading matrix in area 2, involves L_2
  903. 'pr_W_total':pr_W_total # part ratio for across-area loading matrix jointly for area 1 and 2, involves W_1 and W_2
  904. }
  906. return pr
  908. def compute_psv(self):
  909. '''
  910. Compute percentage of shared variance (%sv) in each across- and within-area loading matrix.
  912. Returns:
  913. psv (dict): Dictionary containing the across- and within-area %sv for neurons in each area
  914. '''
  916. W_1,W_2,L_1,L_2 = self.get_loading_matrices()
  917. psi_1,psi_2 = self.params['psi_1'],self.params['psi_2']
  919. shared_across_x1 = np.diag(W_1.dot(W_1.T))
  920. shared_across_x2 = np.diag(W_2.dot(W_2.T))
  921. shared_within_x1 = np.diag(L_1.dot(L_1.T))
  922. shared_within_x2 = np.diag(L_2.dot(L_2.T))
  923. total_x = shared_across_x1 + shared_within_x1 + psi_1
  924. total_y = shared_across_x2 + shared_within_x2 + psi_2
  926. # for area 1
  927. psv_W_1 = (shared_across_x1 / total_x).flatten() * 100
  928. psv_L_1 = (shared_within_x1 / total_x).flatten() * 100
  929. ind_var_x1 = (psi_1 / total_x).flatten() * 100
  930. avg_psv_W_1 = np.mean(psv_W_1)
  931. avg_psv_L_1 = np.mean(psv_L_1)
  933. # for area 2
  934. psv_W_2 = (shared_across_x2 / total_y).flatten() * 100
  935. psv_L_2 = (shared_within_x2 / total_y).flatten() * 100
  936. ind_var_x2 = (psi_2 / total_y).flatten() * 100
  937. avg_psv_W_2 = np.mean(psv_W_2)
  938. avg_psv_L_2 = np.mean(psv_L_2)
  940. # overall
  941. avg_psv_W_total = np.mean(np.concatenate((psv_W_1,psv_W_2)))
  942. avg_psv_L_total = np.mean(np.concatenate((psv_L_1,psv_L_2)))
  944. psv = {
  945. 'psv_W_1':psv_W_1, # percent of across-area shared variance for each neuron in area 1
  946. 'psv_W_2':psv_W_2, # percent of across-area shared variance for each neuron in area 2
  947. 'psv_L_1':psv_L_1, # percent of within-area shared variance for each neuron in area 1
  948. 'psv_L_2':psv_L_2, # percent of within-area shared variance for each neuron in area 2
  949. 'avg_psv_W_1':avg_psv_W_1, # percent of across-area shared variance, averaged across neurons in area 1
  950. 'avg_psv_W_2':avg_psv_W_2, # percent of across-area shared variance, averaged across neurons in area 1
  951. 'avg_psv_L_1':avg_psv_L_1, # percent of across-area shared variance, averaged across neurons in area 1
  952. 'avg_psv_L_2':avg_psv_L_2, # percent of across-area shared variance, averaged across neurons in area 1
  953. 'ind_var_x1':ind_var_x1, # percent of independent variance for each neuron in area 1
  954. 'ind_var_x2':ind_var_x2, # percent of independent variance for each neuron in area 2
  955. 'avg_psv_W_total':avg_psv_W_total, # percent of across-area shared variance, averaged across all neurons
  956. 'avg_psv_L_total':avg_psv_L_total, # percent of within-area shared variance, averaged across all neurons
  957. }
  958. return psv
  960. def compute_metrics(self,cutoff_thresh=0.95):
  961. '''
  962. Wrapper to compute loading similarity, d_shared, part ratio, %sv, and canonical correlations.
  964. Returns:
  965. metrics (dict): Dictionary containing the computed metrics (loading similarity, d_shared, part ratio, %sv, and canonical correlations)
  966. '''
  968. dshared = self.compute_dshared(cutoff_thresh=cutoff_thresh)
  969. psv = self.compute_psv()
  970. pr = self.compute_part_ratio()
  971. ls = self.compute_load_sim()
  972. _,rho = self.get_canonical_directions()
  974. metrics = {
  975. 'dshared':dshared, # dictionary of d_shared metric
  976. 'psv':psv, # dictionary of %sv metric
  977. 'part_ratio':pr, # dictionary of part ratio metric
  978. 'load_sim':ls, # dictionary of loading similarity metric
  979. 'rho':rho # array of canonical correlations
  980. }
  981. if 'cv_rho' in self.params:
  982. metrics['cv_rho'] = self.params['cv_rho'] # array of cross-validated canonical correlations (if crossvalidate() was called)
  984. return metrics

pcca_fa_mdl.py at commit d13882e, under Apache-2.0 · at the source

Overview

Authors: Megan E McDonnell1,2, Akash Umakantha1,2,3, Ryan C Williamson1,2,3,4, Matthew A Smith1,2,5, Byron M Yu1,2,5,6
ORCID iDs: Megan E McDonnell, Matthew A Smith, Byron M Yu
  1. Neuroscience Institute, Carnegie Mellon University, Pittsburgh, PA USA
  2. Center for the Neural Basis of Cognition, Carnegie Mellon University & University of Pittsburgh, Pittsburgh, PA USA
  3. Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA USA
  4. School of Medicine, University of Pittsburgh, Pittsburgh, PA USA
  5. Department of Biomedical Engineering, Carnegie Mellon University, Pittsburgh, PA USA
  6. Department of Electrical and Computer Engineering, Carnegie Mellon University, Pittsburgh, PA USA
Institutions: University of Pittsburgh (United States); Center for the Neural Basis of Cognition (United States); Carnegie Mellon University (United States)
Journal: Nature communications, volume 17, issue 1, article 5088
Dates: received 3 July 2025; accepted 25 March 2026; published online 11 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-71725-0 · PMID 41965893 · PMCID PMC13247044 · OpenAlex W4411283442
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.06.12.659406
Status: code verified
Categories: non-human primate (organism), cognitive (subfield)
Methods: Spectral & time-frequency, Statistics, Machine learning, Preprocessing, Connectivity, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: Neural decoding, Cognitive neuroscience
MeSH: Cognition*, Prefrontal Cortex*, Animals, Macaca mulatta, Male, Neurons, Spatial Memory (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Simons Foundation (NC-GB-CULM-00003241-05, 543065); National Science Foundation (NSF) (NCS BCS 1734916/1954107, NCS DRL 2124066/2123911); U.S. Department of Health &amp; Human Services | NIH | National Institute of Mental Health (R01 MH118929, R01 MH128393); NIBIB NIH HHS (R01 EB026953, T32 EB029365); NEI NIH HHS (R01 EY035896, R01 EY029250); NIMH NIH HHS (R01 MH118929, R01 MH128393); U.S. Department of Health &amp; Human Services | NIH | National Institute of Neurological Disorders and Stroke (RF1 NS127107); U.S. Department of Health &amp; Human Services | NIH | National Institute of Biomedical Imaging and Bioengineering (T32 EB029365, R01 EB026953); NINDS NIH HHS (RF1 NS127107); U.S. Department of Health &amp; Human Services | NIH | National Eye Institute (R01 EY029250, R01 EY035896)
Citations: cited by 4 papers (Europe PMC); 150 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 20 matches between paragraphs and lines of code.

meganmcd13/figures-dual-pfc

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 3b59632fb5b84c377d82a16ae7dd491e70eba4b5, 9 March 2026
Languages: Python (19), Jupyter (16), MATLAB (9)
Size: 52 files, 44 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 16 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (33 files), SciPy (27 files), Matplotlib (16 files), pandas (15 files), scikit-learn (5 files), Parallel Computing Toolbox (1 file), Signal Processing Toolbox (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
45 files
At the source: github.com/meganmcd13/figures-dual-pfc

SmithLabNeuro/pcca_fa

License: Apache-2.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: d13882eca2eee8cdb5e6ee4b0047b32b01f08754, 8 March 2026
Languages: Python (2), Jupyter (1)
Size: 11 files, 3 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (environment.yml), 1 notebook
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (3 files), SciPy (2 files), Matplotlib (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
5 files
At the source: github.com/SmithLabNeuro/pcca_fa

SmithLabNeuro/Ex

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 1a97dc0873abf4a2450f42ed8575305feaff7a16, 28 August 2026
Languages: MATLAB (415), C (9), Java (1), C/C++ (1)
Size: 495 files, 426 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Psychtoolbox (76 files), Statistics and Machine Learning Toolbox (16 files), Image Processing Toolbox (2 files), Signal Processing Toolbox (2 files), export_fig (1 file), Optimization Toolbox (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
427 files
At the source: github.com/SmithLabNeuro/Ex

Zenodo 18913386

License: Apache-2.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (3 files), SciPy (2 files), Matplotlib (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
5 files
At the source:

Zenodo 18929016

License: Apache-2.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (33 files), SciPy (27 files), Matplotlib (16 files), pandas (15 files), scikit-learn (5 files), Parallel Computing Toolbox (1 file), Signal Processing Toolbox (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
45 files
At the source:

Zenodo 18913385

License: Apache-2.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (3 files), SciPy (2 files), Matplotlib (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
5 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-71725-0.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 6 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 523 scripts, each with its path and the digest of its content;
  • 20 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41467-026-71725-0.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 2 keywords, 7 MeSH terms, 10 funders, 137 references.

Cite

This paper

McDonnell, M. E., Umakantha, A., Williamson, R. C., Smith, M. A., & Yu, B. M. (2026). Interactions across hemispheres in prefrontal cortex reflect global cognitive processing. Nature communications, 17(1), 5088. https://doi.org/10.1038/s41467-026-71725-0

BibTeX

@article{mcdonnell2026interactions,
author = {McDonnell, Megan E and Umakantha, Akash and Williamson, Ryan C and Smith, Matthew A and Yu, Byron M},
title = {{Interactions across hemispheres in prefrontal cortex reflect global cognitive processing}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {5088},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71725-0},
url = {https://doi.org/10.1038/s41467-026-71725-0},
pmid = {41965893},
pmcid = {PMC13247044}
}

RIS

TY - JOUR
AU - McDonnell, Megan E
AU - Umakantha, Akash
AU - Williamson, Ryan C
AU - Smith, Matthew A
AU - Yu, Byron M
TI - Interactions across hemispheres in prefrontal cortex reflect global cognitive processing
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/11
VL - 17
IS - 1
SP - 5088
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71725-0
UR - https://doi.org/10.1038/s41467-026-71725-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71725-0",
"type": "article-journal",
"title": "Interactions across hemispheres in prefrontal cortex reflect global cognitive processing",
"container-title": "Nature communications",
"author": [
{
"family": "McDonnell",
"given": "Megan E"
},
{
"family": "Umakantha",
"given": "Akash"
},
{
"family": "Williamson",
"given": "Ryan C"
},
{
"family": "Smith",
"given": "Matthew A"
},
{
"family": "Yu",
"given": "Byron M"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "5088",
"DOI": "10.1038/s41467-026-71725-0",
"PMID": "41965893",
"PMCID": "PMC13247044",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71725-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
11
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.neuron.2026.07.016 [code]
Inferring brain-wide interactions using data-constrained recurrent neural network models.
Journal: Neuron
In common: Parallel Computing Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 2 other tools, 14 references
[2] doi:10.7554/elife.99278 [code]
Brain-wide arousal signals are segregated from movement planning in the superior colliculus of the macaque.
Journal: eLife
In common: NumPy, non-human primate, 11 references, author Matthew A Smith
[3] doi:10.1038/s41467-026-75347-4 [code]
Sleep reveals dynamics integrating and segregating movement and stimulus representations in V1.
Journal: Nature communications
In common: export_fig, Optimization Toolbox, Parallel Computing Toolbox, 8 other tools, 7 references
[4] doi:10.1371/journal.pbio.3003915 [code]
Noise-invariant representations of sound emerge along the canonical cortical hierarchy.
Journal: PLoS biology
In common: Optimization Toolbox, Statistics and Machine Learning Toolbox, scikit-learn, 4 other tools, 9 references
[5] doi:10.1038/s41467-026-75705-2 [code]
Redundant prefrontal hemispheres adapt storage strategy to working memory demands.
Journal: Nature communications
In common: scikit-learn, pandas, SciPy, 2 other tools, non-human primate, cognitive, 6 references, author Matthew A Smith
[6] doi:10.1126/sciadv.adz6495
Pupil-linked arousal heterogeneously modulates cell-type-specific sensory processing.
Journal: Science advances
In common: 13 references
[7] doi:10.1038/s41598-026-55225-1 [code]
Benchmarking criteria to determine latent linear dimensionality in neural data.
Journal: Scientific reports
In common: Statistics and Machine Learning Toolbox, pandas, 11 references
[8] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: export_fig, Psychtoolbox, Optimization Toolbox, 9 other tools
[9] doi:10.1126/sciadv.adv5652 [code]
The anterior cingulate cortex modulates pupil-linked arousal.
Journal: Science advances
In common: pandas, SciPy, Matplotlib, 1 other tool, 9 references
[10] doi:10.1002/hbm.70602 [code]
Neuroimaging Correlates of Post-Stroke Pain After Ischemic Stroke: Secondary Analysis of the INSPiRE-TMS Trial.
Journal: Human brain mapping
In common: export_fig, Psychtoolbox, Optimization Toolbox, 8 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.