OSCR

Single-neuron network topology governs neural computation and learning in primate cortex.

Code ↔ Paper

2 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 2 matches
  1. [1] § Methods › Quantify the influence of single neuron on the activity dynamics within network › Demixed principal component analysis (dPCA) ↔ python/dPCA/dPCA.py, lines 21–93 · score 0.69 · demixed principal component, population activity, dPCA, dimensionality, variance
  2. [2] § Methods › Quantify the evolution of population activity › Alignment index ↔ python/dPCA/dPCA.py, lines 21–93 · score 0.63 · population activity, variance explained, principal component, PCA, axis, dimensional

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 · 40 KB · MIT · 2 matches

  1. """ demixed Principal Component Analysis
  2. """
  4. # Author: Wieland Brendel <[email hidden]>
  5. #
  6. # License: BSD 3 clause
  8. from __future__ import print_function
  10. import numpy as np
  11. from collections import OrderedDict
  12. from itertools import combinations, chain
  13. from scipy.sparse.linalg import svds
  14. from scipy.linalg import pinv
  16. from sklearn.base import BaseEstimator
  17. from sklearn.utils.extmath import randomized_svd
  18. import numexpr as ne
  19. from .utils import shuffle2D, classification, denoise_mask
  21. class dPCA(BaseEstimator):
  22. """ demixed Principal component analysis (dPCA)
  24. dPCA is a linear dimensionality reduction technique that automatically discovers
  25. and highlights the essential features of complex population activities. The
  26. population activity is decomposed into a few demixed components that capture most
  27. of the variance in the data and that highlight the dynamic tuning of the population
  28. to various task parameters, such as stimuli, decisions, rewards, etc.
  30. Parameters
  31. ----------
  32. labels : int or string
  33. Labels of feature axis.
  35. If int the corresponding number of labels are selected from the alphabet 'abcde...'
  37. join : None or dict
  38. Parameter combinations to join
  40. If a data set has parametrized by time t and stimulus s, then dPCA will split
  41. the data into marginalizations corresponding to 't', 's' and 'ts'. At times,
  42. we want to join different marginalizations (like 's' and 'ts'), e.g. if
  43. we are only interested in the time-modulated stimulus components. In this case,
  44. we would pass {'ts' : ['s','ts']}.
  46. regularizer : None, float, 'auto'
  47. Regularization parameter. If None or 0, then no regularization is applied.
  48. For float, the regularization weight is regularizer*var(data). If 'auto', the
  49. optimal regularization parameter is found during fitting (might take some time).
  51. n_components : None, int or dict
  52. Number of components to keep.
  54. If n_components is int, then the same number of components are kept in every
  55. marginalization. Otherwise, the dict allows to set the number of components
  56. in each marginalization (e.g. {'t' : 10, 'ts' : 5}). Defaults to 10.
  58. copy : bool
  59. If False, data passed to fit are overwritten and running
  60. fit(X).transform(X) will not yield the expected results,
  61. use fit_transform(X) instead.
  63. n_iter : int (default: 0)
  64. Number of iterations for randomized SVD solver (sklearn).
  66. Attributes
  67. ----------
  68. explained_variance_ratio_ : dict with arrays, [n_components]
  69. Dictionary in which each key refers to one marginalization and the \
  70. value is a vector with the percentage of variance explained by each of \
  71. the marginal components.
  73. Notes
  74. -----
  75. Implements the dPCA model from:
  76. D Kobak*, W Brendel*, C Constantinidis, C Feierstein, A Kepecs, Z Mainen, \
  77. R Romo, X-L Qi, N Uchida, C Machens
  78. Demixed principal component analysis of population activity in higher \
  79. cortical areas reveals independent representation of task parameters,
  82. Examples
  83. --------
  85. >>> import numpy as np
  86. >>> from dPCA import dPCA
  87. >>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
  88. >>> dpca = dPCA(n_components=2)
  89. >>> dpca.fit(X)
  90. PCA(copy=True, n_components=2, whiten=False)
  91. >>> print(pca.explained_variance_ratio_)
  92. [ 0.99244... 0.00755...]
  93. """
  94. def __init__(self, labels=None, join=None, n_components=10, regularizer=None, copy=True, n_iter=0):
  95. # create labels from alphabet if not provided
  96. if isinstance(labels,str):
  97. self.labels = labels
  98. elif isinstance(labels,int):
  99. alphabet = 'abcdefghijklmnopqrstuvwxyz'
  100. self.labels = alphabet[:labels]
  101. else:
  102. raise TypeError('Wrong type for labels. Please either set labels to the number of variables or provide the axis labels as a single string of characters (like "ts" for time and stimulus)')
  104. self._join = join
  105. self.join = join
  106. self.regularizer = 0 if regularizer == None else regularizer
  107. self.opt_regularizer_flag = regularizer == 'auto'
  108. self.n_components = n_components
  109. self.copy = copy
  110. self.marginalizations = self._get_parameter_combinations()
  111. self.n_iter = n_iter
  113. # set debug mode, 0 = no reports, 1 = warnings, 2 = warnings & progress, >2 = everything
  114. self.debug = 2
  116. if regularizer == 'auto':
  117. print("""You chose to determine the regularization parameter automatically. This can
  118. take substantial time and grows linearly with the number of crossvalidation
  119. folds. The latter can be set by changing self.n_trials (default = 3). Similarly,
  120. use self.protect to set the list of axes that are not supposed to get to get shuffled
  121. (e.g. upon splitting the data into test- and training, time-points should always
  122. be drawn from the same trial, i.e. self.protect = ['t']). This can significantly
  123. speed up the code.""")
  125. self.n_trials = 3
  126. self.protect = None
  128. def fit(self, X, trialX=None):
  129. """Fit the model with X.
  131. Parameters
  132. ----------
  133. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  134. Training data, where n_samples in the number of samples
  135. and n_features_j is the number of the j-features (where the axis correspond
  136. to different parameters).
  138. Returns
  139. -------
  140. self : object
  141. Returns the instance itself.
  142. """
  143. self._fit(X,trialX=trialX)
  144. return self
  146. def fit_transform(self, X, trialX=None):
  147. """Fit the model with X and apply the dimensionality reduction on X.
  149. Parameters
  150. ----------
  151. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  152. Training data, where n_samples in the number of samples
  153. and n_features_j is the number of the j-features (where the axis correspond
  154. to different parameters).
  156. Returns
  157. -------
  158. X_new : dict with arrays with the same shape as X
  159. Dictionary in which each key refers to one marginalization and the value is the
  160. latent component.
  162. """
  163. self._fit(X,trialX=trialX)
  165. return self.transform(X)
  167. def _get_parameter_combinations(self,join=True):
  168. ''' Returns all parameter combinations, e.g. for labels = 'xyz'
  170. {'x' : (0,), 'y' : (1,), 'z' : (2,), 'xy' : (0,1), 'xz' : (0,2), 'yz' : (1,2), 'xyz' : (0,1,2)}
  172. If join == True, parameter combinations are condensed according to self._join, Otherwise all
  173. combinations are returned.
  174. '''
  175. # subsets = () (0,) (1,) (2,) (0,1) (0,2) (1,2) (0,1,2)"
  176. subsets = list(chain.from_iterable(combinations(list(range(len(self.labels))), r) for r in range(len(self.labels))))
  178. # delete empty set & add (0,1,2)
  179. del subsets[0]
  180. subsets.append(list(range(len(self.labels))))
  182. # create dictionary
  183. pcombs = OrderedDict()
  184. for subset in subsets:
  185. key = ''.join([self.labels[i] for i in subset])
  186. pcombs[key] = set(subset)
  188. # condense dict if not None
  189. if isinstance(self._join,dict) and join:
  190. for key, combs in self._join.items():
  191. tmp = [pcombs[comb] for comb in combs]
  193. for comb in combs:
  194. del pcombs[comb]
  196. pcombs[key] = tmp
  198. return pcombs
  200. def _marginalize(self,X,save_memory=False):
  201. """ Marginalize the data matrix
  203. Parameters
  204. ----------
  205. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  206. Training data, where n_samples in the number of samples
  207. and n_features_j is the number of the j-features (where the axis correspond
  208. to different parameters).
  210. save_memory : bool, set to True if memory really is an issue (though optimization is not perfect yet)
  212. Returns
  213. -------
  214. mXs : dictionary, with values corresponding to the marginalized data (and the key refers to the marginalization)
  215. """
  217. def mmean(X,axes,expand=False):
  218. ''' Takes mean along several axis (given as list). If expand the averaged dimensions will be filled with
  219. new axis to retain the dimension.
  220. '''
  221. Z = X.copy()
  223. for ax in np.sort(axes)[::-1]:
  224. Z = np.mean(Z,ax)
  226. if expand == True:
  227. Z = np.expand_dims(Z,ax)
  229. return Z
  231. def dense_marg(Y,mYs):
  232. ''' The original marginalizations as returned by "get_marginalizations" are sparse in the sense that
  233. marginalized axis are newaxis. This functions blows them up to the original size of the data set
  234. (need for optimization).
  235. '''
  236. tmp = np.zeros_like(Y)
  237. for key in list(mYs.keys()):
  239. mYs[key] = (tmp + mYs[key]).reshape((Y.shape[0],-1))
  241. return mYs
  243. Xres = X.copy() # residual of data
  245. # center data
  246. Xres -= np.mean(Xres.reshape((Xres.shape[0],-1)),-1).reshape((Xres.shape[0],) + (len(Xres.shape)-1)*(1,))
  248. # init dict with marginals
  249. Xmargs = OrderedDict()
  251. # get parameter combinations
  252. pcombs = self._get_parameter_combinations(join=False)
  254. # subtract the mean
  255. S = list(pcombs.values())[-1] # full set of indices
  257. if save_memory:
  258. for key, phi in pcombs.items():
  259. S_without_phi = list(S - phi)
  261. # compute marginalization and save
  262. Xmargs[key] = mmean(Xres,np.array(S_without_phi)+1,expand=True)
  264. # subtract the marginalization from the data
  265. Xres -= Xmargs[key]
  266. else:
  267. # efficient precomputation of means
  268. pre_mean = {}
  270. for key, phi in pcombs.items():
  271. if len(key) == 1:
  272. pre_mean[key] = mmean(Xres,np.array(list(phi))+1,expand=True)
  273. else:
  274. pre_mean[key] = mmean(pre_mean[key[:-1]],np.array([list(phi)[-1]])+1,expand=True)
  276. # compute marginalizations
  277. for key, phi in pcombs.items():
  278. key_without_phi = ''.join(filter(lambda ch: ch not in key, self.labels))
  279. # self.labels.translate(None, key)
  281. # build local dictionary for numexpr
  282. X = pre_mean[key_without_phi] if len(key_without_phi) > 0 else Xres
  284. if len(key) > 1:
  285. subsets = list(chain.from_iterable(combinations(key, r) for r in range(1,len(key))))
  286. subsets = [''.join(subset) for subset in subsets]
  287. local_dict = {subset : Xmargs[subset] for subset in subsets}
  288. local_dict['X'] = X
  290. Xmargs[key] = ne.evaluate('X - ' + ' - '.join(subsets),local_dict=local_dict)
  291. else:
  292. Xmargs[key] = X
  294. # condense dict if not None
  295. if isinstance(self._join,dict):
  296. for key, combs in self._join.items():
  297. Xshape = np.ones(len(self.labels)+1,dtype='int')
  298. for comb in combs:
  299. sh = np.array(Xmargs[comb].shape)
  300. Xshape[(sh-1).nonzero()] = sh[(sh-1).nonzero()]
  302. tmp = np.zeros(Xshape)
  304. for comb in combs:
  305. tmp += Xmargs[comb]
  306. del Xmargs[comb]
  308. Xmargs[key] = tmp
  310. Xmargs = dense_marg(X,Xmargs)
  312. return Xmargs
  314. def _optimize_regularization(self,X,trialX,center=True,lams='auto'):
  315. """ Optimization routine to find optimal regularization parameter.
  317. TO DO: Routine is pretty dumb right now (go through predetermined
  318. list and find minimum). There are several ways to speed it up.
  319. """
  321. # center data
  322. if center:
  323. X = X - np.mean(X.reshape((X.shape[0],-1)),1).reshape((X.shape[0],)\
  324. + len(self.labels)*(1,))
  326. # compute variance of data
  327. varX = np.sum(X**2)
  329. # test different inits and regularization parameters
  330. if lams == 'auto':
  331. N = 45
  332. lams = np.logspace(0,N,num=N, base=1.4, endpoint=False)*1e-7
  334. # compute crossvalidated score over n_trials repetitions
  335. scores = self.crossval_score(lams,X,trialX,mean=False)
  337. # take mean over total scores
  338. totalscore = np.mean(np.sum(np.dstack([scores[key] for key in list(scores.keys())]),-1),0)
  340. # Raise warning if optimal lambda lies at boundaries
  341. if np.argmin(totalscore) == 0 or np.argmin(totalscore) == len(totalscore) - 1:
  342. if self.debug > 0:
  343. print("Warning: Optimal regularization parameter lies at the \
  344. boundary of the search interval. Please provide \
  345. different search list (key: lams).")
  347. # set minimum as new lambda
  348. self.regularizer = lams[np.argmin(totalscore)]
  350. if self.debug > 1:
  351. print('Optimized regularization, optimal lambda = ', self.regularizer)
  352. print('Regularization will be fixed; to compute the optimal \
  353. parameter again on the next fit, please \
  354. set opt_regularizer_flag to True.')
  356. self.opt_regularizer_flag = False
  358. def crossval_score(self,lams,X,trialX,mean=True):
  359. """ Calculates crossvalidation scores for a given set of regularization
  360. parameters. To this end it takes one parameter off the list,
  361. computes the model on a training set and then validates the
  362. reconstruction performance on a validation set.
  364. Parameters
  365. ----------
  366. lams: 1D array of floats
  367. Array of regularization parameters to test.
  369. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  370. Training data, where n_samples in the number of samples
  371. and n_features_j is the number of the j-features (where the
  372. axis correspond to different parameters).
  374. trialX: array-like, shape (n_trials, n_samples, n_features_1, n_features_2, ...)
  375. Trial-by-trial data. Shape is similar to X but with an additional axis at the beginning
  376. with different trials. If different combinations of features have different number
  377. of trials, then set n_samples to the maximum number of trials and fill unoccupied data
  378. points with NaN.
  380. mean: bool (default: True)
  381. Set True if the crossvalidation score should be averaged over
  382. all marginalizations, otherwise False.
  384. Returns
  385. -------
  386. mXs : dictionary, with values corresponding to the marginalized
  387. data (and the key refers to the marginalization)
  388. """
  389. # placeholder for scores
  390. scores = np.zeros((self.n_trials,len(lams))) if mean else {key : np.zeros((self.n_trials,len(lams))) for key in list(self.marginalizations.keys())}
  392. # compute number of samples in each condition
  393. N_samples = self._get_n_samples(trialX,protect=self.protect)
  395. for trial in range(self.n_trials):
  396. print("Starting trial ", trial + 1, "/", self.n_trials)
  398. # perform split into training and test trials
  399. trainX, validX = self.train_test_split(X,trialX,N_samples=N_samples)
  401. # compute marginalization of test and validation data
  402. trainmXs, validmXs = self._marginalize(trainX), self._marginalize(validX)
  404. # compute crossvalidation score for every regularization parameter
  405. for k, lam in enumerate(lams):
  406. # fit dpca model
  407. self.regularizer = lam
  408. self._fit(trainX,mXs=trainmXs,optimize=False)
  410. # compute crossvalidation score
  411. if mean:
  412. scores[trial,k] = self._score(validX,validmXs)
  413. else:
  414. tmp = self._score(validX,validmXs,mean=False)
  415. for key in list(self.marginalizations.keys()):
  416. scores[key][trial,k] = tmp[key]
  418. return scores
  420. def _score(self,X,mXs,mean=True):
  421. """ Scoring for crossvalidation. Predicts one observable (e.g. one neuron) of X at a time, using all other dimensions:
  423. \sum_phi ||X[n] - F_\phi D_phi^{-n} X^{-n}||^2
  425. where phi refers to the marginalization and X^{-n} (D_phi^{-n}) are all rows of X (D) except the n-th row.
  426. """
  427. n_features = X.shape[0]
  428. X = X.reshape((n_features,-1))
  430. error = {key: 0 for key in list(mXs.keys())}
  431. PDY = {key : np.dot(self.P[key],np.dot(self.D[key].T,X)) for key in list(mXs.keys())}
  432. trPD = {key : np.sum(self.P[key]*self.D[key],1) for key in list(mXs.keys())}
  434. for key in list(mXs.keys()):
  435. error[key] = np.sum((mXs[key] - PDY[key] + trPD[key][:,None]*X)**2)
  437. return error if not mean else np.sum(list(error.values()))
  439. def _randomized_dpca(self,X,mXs,pinvX=None):
  440. """ Solves the dPCA minimization problem analytically by using a randomized SVD solver from sklearn.
  442. Returns
  443. -------
  444. P : dict mapping strings to array-like,
  445. Holds encoding matrices for each term in variance decompostions (used in inverse_transform
  446. to map from low-dimensional representation back to original data space).
  448. D : dict mapping strings to array-like,
  449. Holds decoding matrices for each term in variance decompostions (used to transform data
  450. to low-dimensional space).
  452. """
  454. n_features = X.shape[0]
  455. rX = X.reshape((n_features,-1))
  456. pinvX = pinv(rX) if pinvX is None else pinvX
  458. P, D = {}, {}
  460. for key in list(mXs.keys()):
  461. mX = mXs[key].reshape((n_features,-1)) # called X_phi in paper
  462. C = np.dot(mX,pinvX)
  464. if isinstance(self.n_components,dict):
  465. U,s,V = randomized_svd(np.dot(C,rX),n_components=self.n_components[key],n_iter=self.n_iter,random_state=np.random.randint(10e5))
  466. else:
  467. U,s,V = randomized_svd(np.dot(C,rX),n_components=self.n_components,n_iter=self.n_iter,random_state=np.random.randint(10e5))
  469. P[key] = U
  470. D[key] = np.dot(U.T,C).T
  472. return P, D
  474. def _add_regularization(self,Y,mYs,lam,SVD=None,pre_reg=False):
  475. """ Prepares the data matrix and its marginalizations for the randomized_dpca solver (see paper)."""
  476. n_features = Y.shape[0]
  478. if not pre_reg:
  479. regY = np.hstack([Y.reshape((n_features,-1)),lam*np.eye(n_features)])
  480. else:
  481. regY = Y
  482. regY[:,-n_features:] = lam*eye(n_features)
  484. if not pre_reg:
  485. regmYs = OrderedDict()
  487. for key in list(mYs.keys()):
  488. regmYs[key] = np.hstack([mYs[key],np.zeros((n_features,n_features))])
  489. else:
  490. regmYs = mYs
  492. if SVD is not None:
  493. U,s,V = SVD
  495. M = ((s**2 + lam**2)**-1)[:,None]*U.T
  496. pregY = np.dot(np.vstack([V.T*s[None,:],lam*U]),M)
  497. else:
  498. pregY = np.dot(regY.reshape((n_features,-1)).T,np.linalg.inv(np.dot(Y.reshape((n_features,-1)),Y.reshape((n_features,-1)).T) + lam**2*np.eye(n_features)))
  500. return regY, regmYs, pregY
  502. def _fit(self, X, trialX=None, mXs=None, center=True, SVD=None, optimize=True):
  503. """ Fit the model on X
  505. Parameters
  506. ----------
  507. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  508. Training data, where n_samples in the number of samples
  509. and n_features_j is the number of the j-features (where the axis correspond
  510. to different parameters).
  512. trialX: array-like, shape (n_trials, n_samples, n_features_1, n_features_2, ...)
  513. Trial-by-trial data. Shape is similar to X but with an additional axis at the beginning
  514. with different trials. If different combinations of features have different number
  515. of trials, then set n_samples to the maximum number of trials and fill unoccupied data
  516. points with NaN.
  518. mXs: dict with values in the shape of X
  519. Marginalized data, should be the result of dpca._marginalize
  521. center: bool
  522. Centers data if center = True
  524. SVD: list of arrays
  525. Singular-value decomposition of the data. Don't provide!
  527. optimize: bool
  528. Flag to turn automatic optimization of regularization parameter on or off. Needed
  529. internally.
  530. """
  532. def flat2d(A):
  533. ''' Flattens all but the first axis of an ndarray, returns view. '''
  534. return A.reshape((A.shape[0],-1))
  536. # X = check_array(X)
  538. n_features = X.shape[0]
  540. # center data
  541. if center:
  542. X = X - np.mean(flat2d(X),1).reshape((n_features,) + len(self.labels)*(1,))
  544. # marginalize data
  545. if mXs is None:
  546. mXs = self._marginalize(X)
  548. # compute optimal regularization
  549. if self.opt_regularizer_flag and optimize:
  550. if self.debug > 0:
  551. print("Start optimizing regularization.")
  553. if trialX is None:
  554. raise ValueError('To optimize the regularization parameter, the trial-by-trial data trialX needs to be provided.')
  556. self._optimize_regularization(X,trialX)
  558. # add regularization
  559. if self.regularizer > 0:
  560. regX, regmXs, pregX = self._add_regularization(X,mXs,self.regularizer*np.sum(X**2),SVD=SVD)
  561. else:
  562. regX, regmXs, pregX = X, mXs, pinv(X.reshape((n_features,-1)))
  564. # compute closed-form solution
  565. self.P, self.D = self._randomized_dpca(regX,regmXs,pinvX=pregX)
  567. def _zero_mean(self,X):
  568. """ Subtracts the mean from each observable """
  569. return X - np.mean(X.reshape((X.shape[0],-1)),1).reshape((X.shape[0],) + (len(X.shape)-1)*(1,))
  571. def _roll_back(self,X,axes,invert=False):
  572. ''' Rolls all axis in list crossval_protect to the end (or inverts if invert=True) '''
  573. rX = X
  574. axes = np.sort(axes)
  576. if invert:
  577. for ax in reversed(axes):
  578. rX = np.rollaxis(rX,-1,start=ax)
  579. else:
  580. for ax in axes:
  581. rX = np.rollaxis(rX,ax,start=len(X.shape))
  583. return rX
  585. def _get_n_samples(self,trialX,protect=None):
  586. """ Computes the number of samples for each parameter combinations (except along protect) """
  587. n_unprotect = len(trialX.shape) - len(protect) - 1 if protect is not None else len(trialX.shape) - 1
  588. n_protect = len(protect) if protect is not None else 0
  590. return trialX.shape[0] - np.sum(np.isnan(trialX[(np.s_[:],) + (np.s_[:],)*n_unprotect + (0,)*n_protect]),0)
  592. def _check_protected(self,X,protect):
  593. ''' Checks if protect == None or, alternatively, if all protected axis are at the end '''
  594. if protect is None:
  595. protected = True
  596. else:
  597. # convert label in index
  598. protect = [self.labels.index(ax) for ax in protect]
  599. if set(protect) == set(np.arange(len(self.labels)-len(protect),len(self.labels))):
  600. protected = True
  601. else:
  602. protected = False
  603. print('Not all protected axis are at the end! While the algorithm will still work, the performance of the shuffling algorithm will substantially decrease due to unavoidable copies.')
  605. return protected
  607. def train_test_split(self,X,trialX,N_samples=None,sample_ax=0):
  608. """ Splits data in training and validation trial. To this end, we select one data-point in each observable for every
  609. combination of parameters (except along protected axis) for the validation set and average the remaining trial-by-trial
  610. data to get the training set.
  612. Parameters
  613. ----------
  614. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  615. Training data, where n_samples in the number of samples
  616. and n_features_j is the number of the j-features (where the axis correspond
  617. to different parameters).
  619. trialX: array-like, shape (n_trials, n_samples, n_features_1, n_features_2, ...)
  620. Trial-by-trial data. Shape is similar to X but with an additional axis at the beginning
  621. with different trials. If different combinations of features have different number
  622. of trials, then set n_samples to the maximum number of trials and fill unoccupied data
  623. points with NaN.
  625. N_samples: array-like with the same shape as X (except for protected axis).
  626. Number of trials in each condition. If None, computed from trial data.
  629. Returns
  630. -------
  631. trainX: array-like, same shape as X
  632. Training data
  634. blindX: array-like, same shape as X
  635. Validation data
  637. """
  638. def flat2d(A):
  639. ''' Flattens all but the first axis of an ndarray, returns view. '''
  640. return A.reshape((A.shape[0],-1))
  642. protect = self.protect
  644. n_samples = trialX.shape[-1] # number of samples
  645. n_unprotect = len(X.shape) - len(protect) if protect is not None else len(X.shape)
  646. n_protect = len(protect) if protect is not None else 0
  648. if sample_ax != 0:
  649. raise NotImplemented('The sample axis needs to come first.')
  651. # test if all protected axes lie at the end
  652. protected = self._check_protected(trialX,protect)
  654. # reorder matrix to protect certain axis (for speedup)
  655. if not(protected):
  656. # turn crossval_protect into index listX
  657. axes = [self.labels.index(ax) + 2 for ax in protect]
  659. # reorder matrix
  660. trialX = self._roll_back(trialX,axes)
  661. X = np.squeeze(self._roll_back(X[None,...],axes))
  663. # compute number of samples in each condition
  664. if N_samples is None:
  665. N_samples = self._get_n_samples(trialX,protect=self.protect)
  667. # get random indices
  668. idx = (np.random.rand(*N_samples.shape)*N_samples).astype(int)
  670. # select values
  671. blindX = np.empty(trialX.shape[1:])
  673. # iterate over multi_index
  674. it = np.nditer(np.empty(N_samples.shape), flags=['multi_index'])
  676. while not it.finished:
  677. blindX[it.multi_index + (np.s_[:],)*n_protect] = trialX[(idx[it.multi_index],) + it.multi_index + (np.s_[:],)*n_protect]
  678. it.iternext()
  680. # compute trainX
  681. trainX = (X*(N_samples/(N_samples-1))[(np.s_[:],)*n_unprotect + (None,)*n_protect] - blindX/(N_samples-1)[(np.s_[:],)*n_unprotect + (None,)*n_protect])
  683. # inverse rolled axis in blindX
  684. if not(protected):
  685. blindX = self._roll_back(blindX[...,None],axes,invert=True)[...,0]
  686. trainX = self._roll_back(trainX[...,None],axes,invert=True)[...,0]
  688. # remean datasets (both equally)
  689. trainX -= np.mean(flat2d(trainX),1)[(np.s_[:],) + (None,)*(len(X.shape)-1)]
  690. blindX -= np.mean(flat2d(blindX),1)[(np.s_[:],) + (None,)*(len(X.shape)-1)]
  692. return trainX, blindX
  694. def shuffle_labels(self,trialX):
  695. """ Shuffles *inplace* labels between conditions in trial-by-trial data, respecting the number of trials per condition.
  697. Parameters
  698. ----------
  699. trialX: array-like, shape (n_trials, n_samples, n_features_1, n_features_2, ...)
  700. Trial-by-trial data. Shape is similar to X but with an additional axis at the beginning
  701. with different trials. If different combinations of features have different number
  702. of trials, then set n_samples to the maximum number of trials and fill unoccupied data
  703. points with NaN.
  705. """
  707. # import shuffling algorithm from cython source
  708. protect = self.protect
  710. # test if all protected axes lie at the end
  711. protected = self._check_protected(trialX,protect)
  713. # reorder matrix to protect certain axis (for speedup)
  714. if not(protected):
  715. # turn crossval_protect into index list
  716. axes = [self.labels.index(ax) + 2 for ax in protect]
  718. # reorder matrix
  719. trialX = self._roll_back(trialX,axes)
  721. # reshape all non-protect axis into one vector
  722. original_shape = trialX.shape
  723. trialX = trialX.reshape((-1,) + trialX.shape[-len(protect):])
  725. # reshape all protected axis into one
  726. original_shape_protected = trialX.shape
  727. trialX = trialX.reshape((trialX.shape[0],-1))
  729. # shuffle within non-protected axis
  730. shuffle2D(trialX)
  732. # inverse reshaping of protected axis
  733. trialX = trialX.reshape(original_shape_protected)
  735. # inverse reshaping & sample axis
  736. trialX = trialX.reshape(original_shape)
  737. #trialX = np.rollaxis(trialX,0,len(original_shape))
  739. # inverse rolled axis in trialX
  740. if protected:
  741. trialX = self._roll_back(trialX,axes,invert=True)
  743. return trialX
  745. def significance_analysis(self,X,trialX,n_shuffles=100,n_splits=100,n_consecutive=1,axis=None,full=False):
  746. '''
  747. Cross-validated significance analysis of dPCA model. Here the generalization from the training
  748. to test data is tested by a simple classification measure in which one tries to predict the
  749. label of a validation test point from the training data. The performance is tested for n_splits
  750. test and training separations. The classification performance is then compared against
  751. the performance on data with randomly shuffled labels. Only if the performance is higher
  752. then the maximum in the shuffled data we regard the component as significant.
  754. Parameters
  755. ----------
  756. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  757. Training data, where n_samples in the number of samples
  758. and n_features_j is the number of the j-features (where the axis correspond
  759. to different parameters).
  761. trialX: array-like, shape (n_trials, n_samples, n_features_1, n_features_2, ...)
  762. Trial-by-trial data. Shape is similar to X but with an additional axis at the beginning
  763. with different trials. If different combinations of features have different number
  764. of trials, then set n_samples to the maximum number of trials and fill unoccupied data
  765. points with NaN.
  767. n_shuffles: integer
  768. Number of label shuffles over which the maximum is taken (default = 100, which
  769. is equivalent to p > 0.01)
  771. n_splits: integer
  772. Number of train-test splits per shuffle, from which the average performance is
  773. deduced.
  775. n_consecutive: integer
  776. Sometimes individual data points are deemed significant purely by chance. To reduced
  777. such noise one can demand that at least n consecutive data points are rated as significant.
  779. axis: None or True (default = None)
  780. Determines whether the significance is calculated over the last axis. More precisely,
  781. one is often interested in determining the significance of a component over time. In this
  782. case, set axis to True and make sure the last axis is time.
  784. full: Boolean (default = False)
  785. Whether or not all scores are returned. If False, only the significance matrix is returned.
  788. Returns
  789. -------
  790. masks: Dictionary
  791. Dictionary with keys corresponding to the marginalizations and with values that are
  792. binary nparrays that capture the significance of the demixed components.
  794. true_score: Dictionary (only returned when full = True)
  795. Dictionary with the scores of the data.
  797. scores: Dictionary (only returned when full = True)
  798. Dictionary with the scores of the shuffled data.
  800. '''
  801. assert axis in [None, True]
  803. def compute_mean_score(X,trialX,n_splits):
  804. K = 1 if axis is None else X.shape[-1]
  806. if type(self.n_components) == int:
  807. scores = {key : np.empty((self.n_components, n_splits, K)) for key in keys}
  808. else:
  809. scores = {key : np.empty((self.n_components[key], n_splits, K)) for key in keys}
  811. for shuffle in range(n_splits):
  812. print('.', end=' ')
  814. # do train-validation split
  815. trainX, validX = self.train_test_split(X,trialX)
  817. # fit a dPCA model to training data & transform validation data
  818. trainZ = self.fit_transform(trainX)
  819. validZ = self.transform(validX)
  821. # reshape data to match Cython input
  822. for key in keys:
  823. ncomps = self.n_components if type(self.n_components) == int else self.n_components[key]
  825. # mean over all axis not in key
  826. axset = self.marginalizations[key]
  827. axset = axset if type(axset) == set else set.union(*axset)
  828. axes = set(range(len(X.shape)-1)) - axset
  829. for ax in list(axes)[::-1]:
  830. trainZ[key] = np.mean(trainZ[key],axis=ax+1)
  831. validZ[key] = np.mean(validZ[key],axis=ax+1)
  833. # reshape
  834. if len(X.shape)-2 in axset and axis is not None:
  835. trainZ[key] = trainZ[key].reshape((ncomps,-1,K))
  836. validZ[key] = validZ[key].reshape((ncomps,-1,K))
  837. else:
  838. trainZ[key] = trainZ[key].reshape((ncomps,-1,1))
  839. validZ[key] = validZ[key].reshape((ncomps,-1,1))
  841. # compute classification score
  842. for key in keys:
  843. ncomps = self.n_components if type(self.n_components) == int else self.n_components[key]
  844. for comp in range(ncomps):
  845. scores[key][comp, shuffle] = classification(trainZ[key][comp],validZ[key][comp])
  847. for key in keys:
  848. scores[key] = np.nanmean(scores[key], axis=1)
  850. return scores
  852. if self.opt_regularizer_flag:
  853. print("Regularization not optimized yet; start optimization now.")
  854. self._optimize_regularization(X,trialX)
  856. keys = list(self.marginalizations.keys())
  857. keys.remove(self.labels[-1])
  859. # shuffling is in-place, so we need to copy the data
  860. trialX = trialX.copy()
  862. # compute score of original data
  863. print("Compute score of data: ", end=' ')
  864. true_score = compute_mean_score(X,trialX,n_splits)
  865. print("Finished.")
  867. # data collection
  868. scores = {key : [] for key in keys}
  870. # iterate over shuffles
  871. for it in range(n_shuffles):
  872. print("\rCompute score of shuffled data: ", str(it), "/", str(n_shuffles), end=' ')
  874. # shuffle labels
  875. self.shuffle_labels(trialX)
  877. # mean trial-by-trial data
  878. X = np.nanmean(trialX,axis=0)
  880. score = compute_mean_score(X,trialX,n_splits)
  882. for key in keys:
  883. scores[key].append(score[key])
  885. # binary mask, if data score is above maximum shuffled score make true
  886. masks = {}
  887. for key in keys:
  888. maxscore = np.amax(np.dstack(scores[key]),-1)
  889. masks[key] = true_score[key] > maxscore
  891. if n_consecutive > 1:
  892. for key in keys:
  893. mask = masks[key]
  895. for k in range(mask.shape[0]):
  896. masks[key][k,:] = denoise_mask(masks[key][k].astype(np.int32),n_consecutive)
  898. if full:
  899. return masks, true_score, scores
  900. else:
  901. return masks
  903. def transform(self, X, marginalization=None):
  904. """Apply the dimensionality reduction on X.
  906. X is projected on the first principal components previous extracted
  907. from a training set.
  909. Parameters
  910. ----------
  911. X: array-like, shape (n_samples, n_features_1, n_features_2, ...)
  912. Training data, where n_samples in the number of samples
  913. and n_features_j is the number of the j-features (where the axis correspond
  914. to different parameters).
  916. marginalization : str or None
  917. Marginalization subspace upon which to project, if None return dict
  918. with projections on all marginalizations
  920. Returns
  921. -------
  922. X_new : dict with arrays of the same shape as X
  923. Dictionary in which each key refers to one marginalization and the value is the
  924. latent component. If specific marginalization is given, returns only array
  926. """
  927. X = self._zero_mean(X)
  928. total_variance = np.sum(X**2)
  929. Xmargs = self._marginalize(X)
  931. def marginal_variances(marginal):
  932. ''' Computes the relative variance explained of each component
  933. within a marginalization
  934. '''
  935. D, P, Xmarg = self.D[marginal], self.P[marginal], Xmargs[marginal]
  936. return [(np.sum(Xmarg**2)-np.sum((Xmarg-np.outer(P[:, k], np.dot(D[:,k], Xmarg)))**2)) / total_variance for k in range(D.shape[1])]
  938. if marginalization is not None:
  939. D, Xr = self.D[marginalization], X.reshape((X.shape[0],-1))
  940. X_transformed = np.dot(D.T, Xr).reshape((D.shape[1],) + X.shape[1:])
  941. self.explained_variance_ratio_ = {marginalization : marginal_variances(marginalization)}
  942. else:
  943. X_transformed = {}
  944. self.explained_variance_ratio_ = {}
  945. for key in list(self.marginalizations.keys()):
  946. X_transformed[key] = np.dot(self.D[key].T, X.reshape((X.shape[0],-1))).reshape((self.D[key].shape[1],) + X.shape[1:])
  947. self.explained_variance_ratio_[key] = marginal_variances(key)
  949. return X_transformed
  951. def inverse_transform(self, X, marginalization):
  952. """ Transform data back to its original space, i.e.,
  953. return an input X_original whose transform would be X
  955. Parameters
  956. ----------
  957. X : array-like, shape (n_samples, n_components)
  958. New data, where n_samples is the number of samples
  959. and n_components is the number of components.
  961. Returns
  962. -------
  963. X_original array-like, shape (n_samples, n_features)
  964. """
  965. X = self._zero_mean(X)
  966. X_transformed = np.dot(self.P[marginalization],X.reshape((X.shape[0],-1))).reshape((self.P[marginalization].shape[0],) + X.shape[1:])
  968. return X_transformed
  970. def reconstruct(self, X, marginalization):
  971. """ Transform data first into reduced space before projecting
  972. it back into data space. Equivalent to inverse_transform(transform(X)).
  974. Parameters
  975. ----------
  976. X : array-like, shape (n_samples, n_components)
  977. New data, where n_samples is the number of samples
  978. and n_components is the number of components.
  980. Returns
  981. -------
  982. X_original array-like, shape (n_samples, n_features)
  983. """
  984. return self.inverse_transform(self.transform(X,marginalization),marginalization)

dPCA.py at commit 1def5b1, under MIT · at the source

Overview

Authors: Zhuangyi Jiang1, Ziang Liu1, Li Shi1, Fang Fang1,2,3,4, Shiming Tang3,5, Yang Zhou1,2,3,4
ORCID iDs: Zhuangyi Jiang, Ziang Liu, Fang Fang, Shiming Tang, Yang Zhou
  1. Peking-Tsinghua Center for Life Sciences, Peking University,Beijing, China
  2. School of Psychological and Cognitive Sciences, Peking University,Beijing, China
  3. PKU-IDG/McGovern Institute for Brain Research, Peking University,Beijing, China
  4. School of Life Sciences, Peking University,Beijing, China
  5. Beijing Key Laboratory of Behavior and Mental Health, Beijing, China
Institutions: Peking University (China)
Journal: Nature communications, volume 17, issue 1, article 5909
Dates: received 21 September 2025; accepted 15 April 2026; published online 30 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-72510-9 · PMID 42062292 · PMCID PMC13338231 · OpenAlex W7159689184
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: non-human primate (organism), cognitive (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Graphs, fMRI & imaging, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: Cortex, Cognitive neuroscience
MeSH: Learning*, Nerve Net*, Neurons*, Parietal Lobe*, Animals, Macaca mulatta, Male, Models, Neurological (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Ministry of Science and Technology of the People's Republic of China (STI2030-Major Projects (2021ZD0203800)); National Natural Science Foundation of China (National Science Foundation of China) (NSFC32171036); Ministry of Education of the People's Republic of China (Ministry of Education of China) (JYB2025XDXM504)
Citations: not cited yet (Europe PMC); 95 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 2 matches between paragraphs and lines of code.

machenslab/dPCA

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 1def5b15854638811a25257dc0b68074ab5a1be0, 29 July 2026
Languages: MATLAB (15), Python (4), Jupyter (1)
Size: 28 files, 20 scripts
Software Heritage: not archived
Found in: the text, “Demixed principal component analysis (dPCA)”
Holds: README, license file, environment (python/pyproject.toml, python/requirements.txt, python/setup.py), continuous integration, 1 notebook
Not found: CITATION.cff, tests, documentation
Tools: NumPy (4 files), Statistics and Machine Learning Toolbox (3 files), Numba (2 files), scikit-learn (2 files), SciPy (2 files)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
22 files
At the source: github.com/machenslab/dPCA

codeocean:8493005

License: none: the authors keep all their rights
State: cannot be verified, verified on 30 September 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 30 September 2026: cannot be verified
  • 30 September 2026: cannot be verified
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:

  • it points to the authors' code: codeocean:8493005

Read it in the paper: doi.org/10.1038/s41467-026-72510-9.

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 20 scripts, each with its path and the digest of its content;
  • 2 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:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1038/s41467-026-72510-9.

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, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 2 keywords, 8 MeSH terms, 3 funders, 94 references.

Cite

This paper

Jiang, Z., Liu, Z., Shi, L., Fang, F., Tang, S., & Zhou, Y. (2026). Single-neuron network topology governs neural computation and learning in primate cortex. Nature communications, 17(1), 5909. https://doi.org/10.1038/s41467-026-72510-9

BibTeX

@article{jiang2026single,
author = {Jiang, Zhuangyi and Liu, Ziang and Shi, Li and Fang, Fang and Tang, Shiming and Zhou, Yang},
title = {{Single-neuron network topology governs neural computation and learning in primate cortex}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {5909},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-72510-9},
url = {https://doi.org/10.1038/s41467-026-72510-9},
pmid = {42062292},
pmcid = {PMC13338231}
}

RIS

TY - JOUR
AU - Jiang, Zhuangyi
AU - Liu, Ziang
AU - Shi, Li
AU - Fang, Fang
AU - Tang, Shiming
AU - Zhou, Yang
TI - Single-neuron network topology governs neural computation and learning in primate cortex
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/30
VL - 17
IS - 1
SP - 5909
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-72510-9
UR - https://doi.org/10.1038/s41467-026-72510-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-72510-9",
"type": "article-journal",
"title": "Single-neuron network topology governs neural computation and learning in primate cortex",
"container-title": "Nature communications",
"author": [
{
"family": "Jiang",
"given": "Zhuangyi"
},
{
"family": "Liu",
"given": "Ziang"
},
{
"family": "Shi",
"given": "Li"
},
{
"family": "Fang",
"given": "Fang"
},
{
"family": "Tang",
"given": "Shiming"
},
{
"family": "Zhou",
"given": "Yang"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "5909",
"DOI": "10.1038/s41467-026-72510-9",
"PMID": "42062292",
"PMCID": "PMC13338231",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-72510-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
30
]
]
}
}

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.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: Numba, Statistics and Machine Learning Toolbox, scikit-learn, 2 other tools, 7 references
[2] doi:10.1038/s41467-026-74466-2 [code]
Neuromorphic hierarchical modular reservoirs.
Journal: Nature communications
In common: Numba, scikit-learn, SciPy, 1 other tool, 7 references
[3] doi:10.1038/s41467-026-71725-0 [code]
Interactions across hemispheres in prefrontal cortex reflect global cognitive processing.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, scikit-learn, SciPy, 1 other tool, non-human primate, cognitive, 5 references
[4] doi:10.1007/s12311-026-02042-x
The Cerebellar Connectome.
Journal: Cerebellum (London, England)
In common: 9 references
[5] doi:10.7554/elife.107518 [code]
Continuous flash suppression of neural responses and population orientation coding in macaque V1.
Journal: eLife
In common: Statistics and Machine Learning Toolbox, non-human primate, 2 references, author Shi-Ming Tang
[6] doi:10.1038/s41467-026-75585-6 [code]
Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease.
Journal: Nature communications
In common: scikit-learn, SciPy, NumPy, 6 references
[7] doi:10.1038/s41467-026-71458-0 [code]
Early differential impact of MeCP2 mutations on functional networks in Rett syndrome patient-derived human cortical organoids.
Journal: Nature communications
In common: Numba, Statistics and Machine Learning Toolbox, scikit-learn, 2 other tools, 4 references
[8] doi:10.1002/hbm.70485 [code]
Exploring the Role of the Rich Club in Network Control of Neurocognitive States.
Journal: Human brain mapping
In common: Statistics and Machine Learning Toolbox, scikit-learn, SciPy, 1 other tool, 5 references
[9] doi:10.1038/s42003-026-09831-4 [code]
Comprehensive large-scale analyses reveal association between brain structure and cognitive ability during adolescence.
Journal: Communications biology
In common: Statistics and Machine Learning Toolbox, scikit-learn, SciPy, 1 other tool, 5 references
[10] doi:10.1371/journal.pbio.3003915 [code]
Noise-invariant representations of sound emerge along the canonical cortical hierarchy.
Journal: PLoS biology
In common: Numba, Statistics and Machine Learning Toolbox, scikit-learn, 2 other tools, 3 references

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.