OSCR

Intrinsic cortical geometry is associated with individual differences in local functional organization

Code ↔ Paper

5 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 5 matches
  1. [1] § Methods › Shared geometry embedding ↔ src/variograd_utils/embed_utils.py, lines 78–156 · score 0.68 · Laplacian eigenmaps, Cauchy kernel, reference matrices, affinity, embedding
  2. [2] § Methods › Functional connectivity gradients ↔ scripts/07a.connectivity_JE.py, lines 1–68 · score 0.67 · connectivity matrices, joint embedding, functional connectivity, thresholded, mapping, cortical
  3. [3] § Methods › MRI data ↔ scripts/06.downsample_timeseries.py, lines 1–62 · score 0.63 · fMRI, preprocessing, split, FWHM, smoothed, downsampled
  4. [4] § Methods › Functional connectivity gradients ↔ src/variograd_utils/embed_utils.py, lines 78–156 · score 0.59 · Pearson correlations, joint embedding, affinity, mapping, matrices
  5. [5] § Methods › MRI data ↔ src/variograd_utils/brain_utils.py, lines 238–281 · score 0.54 · fMRI, cortical surfaces, HCP, mesh, hemisphere, vertex

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 · 765 lines · 26 KB · MIT · 2 matches

  1. import warnings
  2. import numpy as np
  3. from scipy.sparse import diags
  4. from scipy.linalg import orthogonal_procrustes
  5. from scipy.optimize import linear_sum_assignment
  6. from sklearn.metrics.pairwise import cosine_similarity, euclidean_distances
  7. from sklearn.decomposition import TruncatedSVD
  8. from variograd_utils.core_utils import vector_wise_corr
  10. class JointEmbedding:
  11. """
  12. JointEmbedding class for computing the joint embedding of two matrices.
  14. Parameters:
  15. ----------
  16. method : str, optional
  17. The embedding method to use. Options are:
  18. - "dme": diffusion map embedding (default)
  19. - "le": Laplacian eigenmap
  20. n_components : int, optional
  21. Number of components to compute (default=2).
  22. alignment : str, optional
  23. The alignment method to use:
  24. - "procrustes": orthogonal Procrustes rotation with scaling.
  25. - "rotation": orthogonal Procrustes rotation (default)
  26. - "sort_flip": reorder the embedding dimensions maximizing absolute
  27. correlation coefficients with th independent refrence. Then multiply
  28. dimensions with a negative coefficient by -1.
  29. - "dot_product": align the embedding using the dot product of the joint reference
  30. and independent reference embeddings.
  31. random_state : int, optional
  32. Random seed of the SVDs.
  33. copy : bool, optional
  34. Whether to copy the input matrices (default=True).
  36. Attributes:
  37. ----------
  38. method : str
  39. The embedding method to use.
  40. n_components : int
  41. Number of components to compute.
  42. alignment : str
  43. The alignment method to use.
  44. random_state : int
  45. Random seed of the SVDs.
  46. copy : bool
  47. Whether to copy the input matrices.
  48. vectors : np.ndarray
  49. The eigenvectors of the embedding.
  50. lambdas : np.ndarray
  51. The eigenvalues of the embedding.
  52. independent_ref : np.ndarray
  53. The independent reference embedding used for alignment.
  55. Methods:
  56. --------
  57. fit_transform(M, R, C=None, affinity="cosine", scale=None, method_kwargs=None)
  58. Compute the joint embedding of M and R using the specified method.
  59. _joint_affinity_matrix(M, R, C=None, affinity="cosine", scale=None)
  60. Computes the joint affinity matrix.
  61. _align_embeddings(embedding, joint_reference, independent_reference, method="rotation")
  62. Align the joint embedding with the independently computed reference embedding.
  63. _affinity_matrix(M, method="cosine", scale=None)
  64. Compute the joint affinity matrix of the input data.
  65. """
  67. def __init__(self, method="dme", n_components=2, alignment=None,
  68. random_state=None, copy=True):
  69. self.method = method
  70. self.n_components = n_components
  71. self.random_state = random_state
  72. self.copy = copy
  73. self.alignment = alignment
  74. self.vectors = self.lambdas = None
  75. self.independent_ref = None
  78. def fit_transform(self, M, R, C=None, affinity="cosine", scale=None, method_kwargs=None):
  79. """
  80. Compute the joint embedding of M and R using the specified method.
  82. Parameters:
  83. ----------
  84. M : np.ndarray
  85. Target matrix to embed.
  86. R : np.ndarray
  87. Reference matrix.
  88. C : np.ndarray, optional
  89. Correspondence matrix.
  90. If affinity is not "precomputed", C is used as is and the
  91. specified affinity method is applied only to while M and R.
  92. affinity : str, optional
  93. The method to compute the affinity matrices. Options are:
  94. - "cosine": cosine similarity (default)
  95. - "correlation": Pearson correlation coefficient
  96. - "linear": linear kernel
  97. - "cauchy": Cauchy kernel
  98. - "gauss": Gaussian kernel
  99. - "precomputed": precomputed affinity matrix.
  100. In this case, M and R are assumed to be affinity matrices
  101. and a correspondence matrix C must be specified.
  102. scale : float, optional
  103. The scaling parameter for the kernel methods.
  104. method_kwargs : dict, optional
  105. Additional keyword arguments for the embedding method.
  107. Returns:
  108. -------
  109. self : JointEmbedding
  110. Fitted instance of JointEmbedding.
  111. embedding_M : np.ndarray
  112. The joint embedding of M.
  113. embedding_R : np.ndarray
  114. The joint embedding of R.
  116. Raises:
  117. -------
  118. ValueError
  119. If the affinity is "precomputed" and C is not specified.
  120. """
  122. if (affinity == "precomputed") & (C is None):
  123. raise ValueError("Precomputed affinity assumes M and R are already affinity matrices,"
  124. + "so a correspondance affinity matrx C must be specified too.")
  125. n = M.shape[0]
  127. R = np.array(R, copy=self.copy)
  128. M = np.array(M, copy=self.copy)
  129. if C is not None:
  130. C = np.array(C, copy=self.copy)
  132. A = self._joint_affinity_matrix(M, R, C=C, affinity=affinity, scale=scale)
  134. method_kwargs = {} if method_kwargs is None else method_kwargs
  135. embedding_function = diffusion_map_embedding if self.method == "dme" else laplacian_eigenmap
  136. embedding, vectors, lambdas = embedding_function(A, n_components=self.n_components,
  137. random_state=self.random_state,
  138. **method_kwargs)
  140. self.vectors = vectors
  141. self.lambdas = lambdas
  142. embedding_R = embedding[:n]
  143. embedding_M = embedding[n:]
  145. if self.alignment is not None:
  146. A = _affinity_matrix(R,method=affinity, scale=scale) if affinity != "precomputed" else R
  147. embedding_R_ind, _, _ = embedding_function(A, n_components=self.n_components,
  148. random_state=self.random_state,
  149. **method_kwargs)
  151. embedding_M, embedding_R = self._align_embeddings(embedding_M, embedding_R,
  152. embedding_R_ind,method=self.alignment)
  153. self.independent_ref = embedding_R_ind
  156. return embedding_M, embedding_R
  159. def _joint_affinity_matrix(self, M, R, C=None, affinity="cosine", scale=None):
  160. """
  161. Computes the joint affinity matrix.
  163. Parameters:
  164. ----------
  165. M : np.ndarray
  166. Target matrix to embed.
  167. R : np.ndarray
  168. Reference matrix.
  169. C : np.ndarray, optional
  170. Correspondence matrix of shape (n_samples_M, n_samples_R).
  171. If affinity is not "precomputed", C is used as is and the
  172. specified affinity method is applied only to while M and R.
  173. affinity : str, optional
  174. The method to compute the affinity matrices. Options are:
  175. - "cosine": cosine similarity (default)
  176. - "correlation": Pearson correlation coefficient
  177. - "linear": linear kernel
  178. - "cauchy": Cauchy kernel
  179. - "gauss": Gaussian kernel
  180. - "precomputed": precomputed affinity matrix.
  181. scale : float, optional
  182. The scaling parameter for the kernel methods.
  184. Returns:
  185. -------
  186. A : np.ndarray
  187. The joint affinity matrix.
  188. """
  190. if C is None:
  191. if affinity in ["cosine", "correlation"]:
  192. A = np.vstack([R, M])
  193. A = _affinity_matrix(A, method=affinity, scale=scale)
  195. elif affinity in ["linear", "cauchy", "gauss"]:
  196. M_aff = _affinity_matrix(M, method=affinity, scale=scale)
  197. R_aff = _affinity_matrix(R, method=affinity, scale=scale)
  198. C = pseudo_sqrt(np.dot(M_aff, R_aff), n_components=100)
  199. A = np.block([[R_aff, C.T],
  200. [C, M_aff]])
  202. else:
  203. if affinity == "precomputed":
  204. A = np.block([[R, C.T],
  205. [C, M]])
  206. else:
  207. A = np.block([[_affinity_matrix(R, method=affinity, scale=scale), C.T],
  208. [C, _affinity_matrix(M, method=affinity, scale=scale)]])
  210. return A
  213. def _align_embeddings(self, embedding, joint_reference, independent_reference,
  214. method="rotation"):
  215. """
  216. Align the joint embedding with the independently computed reference embedding.
  218. Parameters:
  219. ----------
  220. embedding : np.ndarray
  221. The joint embedding.
  222. reference : np.ndarray
  223. The reference embedding.
  224. method : str, optional
  225. The method used to align the joint embedding to the independently
  226. computed reference embedding:
  227. - "procrustes": orthogonal Procrustes rotation with scaling.
  228. - "rotation": orthogonal Procrustes rotation (default)
  229. - "sort_flip": reorder the embedding dimensions maximizing absolute
  230. correlation coefficients with th independent refrence. Then multiply
  231. dimensions with a negative coefficient by -1.
  232. - "dot_product": align the embedding using the dot product of the joint reference
  233. and independent reference embeddings.
  235. Returns:
  236. -------
  237. embedding : np.ndarray
  238. The aligned joint embedding.
  239. reference : np.ndarray
  240. The aligned reference embedding.
  242. """
  244. s = 1
  245. if method == "sort_flip":
  246. idx = argsort_axes(joint_reference.copy(), independent_reference.copy())
  247. joint_reference = joint_reference[:, idx]
  248. embedding = embedding[:, idx]
  249. to_flip = vector_wise_corr(joint_reference.copy(), independent_reference.copy()) < 0
  250. R = np.diag([-1 if i else 1 for i in to_flip])
  252. elif method == "dot_product":
  253. R = dot_product_rotation(joint_reference.copy(), independent_reference.copy())
  255. elif method in ["rotation", "procrustes"]:
  256. R, s = procrustes_rotation(joint_reference.copy(), independent_reference.copy())
  257. s = 1 if method == "rotation" else s
  259. else:
  260. raise ValueError(f"Unknown alignment method: {self.alignment}")
  262. joint_reference = np.dot(joint_reference, R) * s
  263. embedding = np.dot(embedding, R) * s
  266. return embedding, joint_reference
  269. def _affinity_matrix(M, method="cosine", scale=None):
  270. """
  271. Convert a distance measure to affinity.
  273. Parameters
  274. ----------
  275. M : array-like
  276. The input matrix of shape (samples, features)
  277. method : str
  278. The method to compute the affinity matrix. Options are:
  279. - "cosine": cosine similarity
  280. - "correlation": Pearson correlation coefficient
  281. - "linear": linear kernel
  282. - "cauchy": Cauchy kernel
  283. - "gauss": Gaussian kernel
  285. scale : float
  286. The scaling parameter for the kernel
  288. Returns
  289. -------
  290. A : array-like
  291. The affinity matrix
  292. """
  295. if method == "cosine":
  296. A = cosine_similarity(M)
  298. elif method == "correlation":
  299. A = np.corrcoef(M)
  301. elif method in {"linear", "cauchy", "gauss"}:
  302. # A = M if _is_square(M) else euclidean_distances(M)
  303. A = kernel_affinity(A, kernel=method, scale=scale)
  305. else:
  306. raise ValueError("Unknown affinity method")
  308. return A
  311. def kernel_affinity(A, kernel="linear", scale=None):
  312. """
  313. Apply kernel to a matrix A
  315. Parameters
  316. ----------
  317. A : array-like
  318. The input matrix of shape smaples x features
  319. kernel : str
  320. The kernel to apply. Options are "cauchy", "gauss", "linear"
  321. scale : float
  322. The scaling parameter for the kernel
  324. Returns
  325. -------
  326. A : array-like
  327. The kernelized matrix
  328. """
  330. if scale is None:
  331. scale = 1 / A.shape[1]
  333. if kernel == "cauchy":
  334. A = 1.0 / (1.0 + (A ** 2) / (scale ** 2))
  336. elif kernel == "gauss":
  337. A = np.exp(-0.5 * (A ** 2) / (scale ** 2))
  339. elif kernel == "linear":
  340. A = 1 / (1 + A / scale)
  342. else:
  343. raise ValueError("Unknown kernel type")
  345. return A
  348. def spectral_affinity(M, R, n_components=2, random_state=None):
  349. """
  350. Compute the spectral similarity between two matrices using Laplacian Eigenmaps and Procrustes alignment.
  352. This function calculates the cosine similarity between Laplacian Eigenmaps of two matrices,
  353. `M` and `R` of shape observations x features. The embeddings are aligned using
  354. the Orthogonal Procrustes transformation before computing similarity.
  356. Parameters
  357. ----------
  358. M : numpy.ndarray
  359. The first input matrix.
  360. R : numpy.ndarray
  361. The second input matrix (also target of the procrustes alignment).
  362. n_components : int, optional
  363. Number of dimensions for the Laplacian Eigenmap embeddings. Default is 2.
  364. random_state : int, optional
  365. Determines random number generation for eigenmap embedding. Default is None.
  367. Returns
  368. -------
  369. numpy.ndarray
  370. A cosine similarity matrix representing the similarity between aligned
  371. embeddings of `M` and `R`.
  372. """
  374. embedding_M, _, _ = laplacian_eigenmap(M, n_components=n_components, normalized=True, random_state=random_state)
  375. embedding_R, _, _ = laplacian_eigenmap(R, n_components=n_components, normalized=True, random_state=random_state)
  377. embedding_M -= embedding_M.mean(axis=0)
  378. embedding_M /= np.linalg.norm(embedding_M)
  379. embedding_R -= embedding_R.mean(axis=0)
  380. embedding_R /= np.linalg.norm(embedding_R)
  382. rotation, scaling = orthogonal_procrustes(embedding_M, embedding_R)
  383. embedding_M = np.dot(embedding_M, rotation) * scaling
  385. A = cosine_similarity(embedding_M, embedding_R)
  387. return A
  390. def diffusion_map_embedding(A, n_components=2, alpha=0.5, diffusion_time=1, random_state=None):
  391. """
  392. Computes the joint diffusion map embedding of an affinity matrix.
  394. Parameters:
  395. ----------
  396. A : np.ndarray
  397. Target matrix to embed.
  398. n_components: int, optional
  399. Number of components to keep (default=2).
  400. alpha: float, optional
  401. Controls Laplacian normalization, balancing local and global structures in embeddings.
  402. Alpha <= 0.5 emphasize local structures, alpha >= 1 focus on global organization.
  403. diffusion_time: float, optional
  404. Determines the scale of the random walk process (default=1).
  405. random_state: int, optional
  406. Random seed of the SVD.
  408. Returns
  409. -------
  410. embedding : numpy.ndarray of shape (n_samples, n_components)
  411. The diffusion map embedding computed from the random walk Laplacian.
  412. vectors : numpy.ndarray of shape (n_samples, n_components)
  413. The corresponding singular vectors (or eigenvectors) used for the embedding.
  414. lambdas : numpy.ndarray of shape (n_components,)
  415. The singular values (or eigenvalues) associated with the embedding.
  417. """
  419. if np.any(A < 0):
  420. warnings.warn(f"{np.sum(A < 0)} negative values in the affinity matrix set to 0. "
  421. + f"Mean negative value: {np.mean(A[A < 0])}({np.std(A[A < 0])})",
  422. RuntimeWarning)
  423. A[A<0] = 0
  425. L = _random_walk_laplacian(A, alpha=alpha)
  427. embedding, vectors, lambdas = _diffusion_map(L, n_components=n_components,
  428. diffusion_time=diffusion_time,
  429. random_state=random_state)
  431. return embedding, vectors, lambdas
  434. def _diffusion_map(L, n_components=2, diffusion_time=1, random_state=None):
  435. """
  436. Computes the diffusion map of the random walk Laplacian L.
  438. Parameters:
  439. ----------
  440. L: np.ndarray
  441. Random walk Laplacian.
  442. n_components: int, optional
  443. Number of components to compute (default=2).
  444. diffusion_time: float, optional
  445. Determines the scale of the random walk process (default=1).
  446. random_state: int, optional
  447. Random seed of the SVD.
  449. Returns
  450. -------
  451. embedding : numpy.ndarray of shape (n_samples, n_components)
  452. The diffusion map embedding computed from the random walk Laplacian.
  453. vectors : numpy.ndarray of shape (n_samples, n_components)
  454. The corresponding singular vectors (or eigenvectors) used for the embedding.
  455. lambdas : numpy.ndarray of shape (n_components,)
  456. The singular values (or eigenvalues) associated with the embedding.
  457. """
  459. n_components = n_components + 1
  460. embedding = TruncatedSVD(n_components=n_components,
  461. random_state=random_state
  462. ).fit(L)
  464. vectors = embedding.components_.T
  465. lambdas = embedding.singular_values_
  466. if np.any(vectors[:, 0] == 0):
  467. warnings.warn("0 values found in the first eigenvector; 1e-15 was added to all vectors.",
  468. RuntimeWarning)
  469. vectors += 1e-16
  471. psi = vectors / np.tile(vectors[:, 0], (vectors.shape[1], 1)).T
  472. lambdas[1:] = np.power(lambdas[1:], diffusion_time)
  473. embedding = psi[:, 1:n_components] @ np.diag(lambdas[1:n_components], 0)
  474. lambdas = lambdas[1:]
  475. vectors = vectors[:, 1:]
  477. return embedding, vectors, lambdas
  480. def _random_walk_laplacian(A, alpha=0.5):
  481. """
  482. Computes the random walk Laplacian of an affinity matrix A.
  484. Parameters:
  485. ----------
  486. A: np.ndarray
  487. Affinity matrix.
  488. alpha: float, optional
  489. Controls Laplacian normalization, balancing local and global structures in embeddings.
  490. Alpha <= 0.5 emphasize local structures, alpha >= 1 focus on global organization.
  492. Returns:
  493. -------
  494. L : np.ndarray
  495. The random walk Laplacian of A.
  496. """
  498. if (A.shape[0] != A.shape[1]) or not np.allclose(A, A.T):
  499. raise ValueError("A must be a squared, symmetrical affinity matrix.")
  501. # Compute the normalized Laplacian
  502. degree = np.sum(A, axis=1)
  503. d_alpha = diags(np.power(degree, -alpha))
  504. L_alpha = d_alpha @ A @ d_alpha
  506. # Compute the random walk Laplacian
  507. degree = np.sum(L_alpha, axis=1)
  508. d_alpha = diags(np.power(degree, -1))
  509. L = d_alpha @ L_alpha
  511. return L
  514. def laplacian_eigenmap(A, n_components=2, normalized=True, random_state=None):
  515. """
  516. Computes the spectral embedding of an affinity matrix.
  518. Parameters:
  519. ----------
  520. A : np.ndarray
  521. Target matrix to embed.
  522. n_components: int, optional
  523. Number of components to compute (default=2).
  524. normalized: bool, optional
  525. Whether to normalize the Laplacian matrix (default=True).
  526. random_state: int, optional
  527. Random seed of the SVD.
  529. Returns
  530. -------
  531. embedding : numpy.ndarray of shape (n_samples, n_components)
  532. The laplaician eigenmap embedding of L.
  533. vectors : numpy.ndarray of shape (n_samples, n_components)
  534. The corresponding singular vectors (or eigenvectors) used for the embedding.
  535. lambdas : numpy.ndarray of shape (n_components,)
  536. The singular values (or eigenvalues) associated with the embedding.
  538. """
  540. if np.any(A < 0):
  541. warnings.warn(f"{np.sum(A < 0)} negative values in the affinity matrix set to 0. "
  542. + f"Mean negative value: {np.mean(A[A < 0])}({np.std(A[A < 0])})",
  543. RuntimeWarning)
  544. A[A<0] = 0
  546. L = _laplacian(A, normalized=normalized)
  548. n_components = n_components + 1
  549. svd = TruncatedSVD(n_components=n_components, random_state=random_state)
  550. embedding = svd.fit_transform(L)[:, 1:]
  551. vectors = svd.components_.T[:, 1:]
  552. lambdas = svd.singular_values_[1:]
  554. return embedding, vectors, lambdas
  557. def _laplacian(A, normalized=True):
  558. """
  559. Computes the Laplacian of an affinity matrix A.
  561. Parameters:
  562. ----------
  563. A : np.ndarray
  564. Affinity matrix.
  565. normalized : bool, optional
  566. Compute the normalized Laplacian (default=True).
  568. Returns:
  569. -------
  570. L : np.ndarray
  571. Laplacian matrix of M.
  572. """
  574. if (A.shape[0] != A.shape[1]) or not np.allclose(A, A.T):
  575. raise ValueError("A must be a squared, symmetrical affinity matrix.")
  577. # Calculate Laplacian
  578. D = np.sum(A, axis=1)
  579. L = D - A
  581. # Normalize Laplacian
  582. if normalized:
  583. D_inv_sqrt = 1.0 / np.sqrt(D)
  584. D_inv_sqrt[np.isinf(D_inv_sqrt)] = 0
  585. D_inv_sqrt = np.diag(D_inv_sqrt)
  586. L = D_inv_sqrt @ (L @ D_inv_sqrt)
  588. return L
  591. def pseudo_sqrt(X, n_components=100):
  592. """
  593. Compute the pseudo-square root of a symmetric matrix.
  595. This function computes a low-rank approximation of the pseudo-square root of
  596. a symmetric matrix using Singular Value Decomposition (SVD).
  598. Parameters
  599. ----------
  600. X : numpy.ndarray
  601. The input square matrix to compute the pseudo-square root for.
  602. n_components : int, optional
  603. The number of components to retain in the low-rank approximation.
  604. Default is 100.
  606. Returns
  607. -------
  608. numpy.ndarray
  609. The approximated symmetric pseudo-square root of the input matrix.
  611. Notes
  612. -----
  613. - If `X` is not symmetric, it is symmetrized as `(X + X.T) / 2`.
  614. - The pseudo-square root is computed as `U * sqrt(S) * U.T`, where `U` and `S`
  615. - The method assumes `X` has a valid SVD decomposition.
  617. Raises
  618. ------
  619. ValueError
  620. If the input matrix `X` is not square.
  622. """
  624. if (X.shape[0] != X.shape[1]) or (X.ndim != 2):
  625. raise ValueError("Input matrix X must be square.")
  627. if not np.allclose(X, X.T, atol=1e-10):
  628. X = (X + X.T) / 2
  630. svd = TruncatedSVD(n_components=n_components)
  631. svd.fit(X)
  632. U = svd.components_.T
  633. S = np.diag(svd.singular_values_)
  634. Ssqrt = S ** 0.5
  636. return np.dot(np.dot(U, Ssqrt), U.T)
  640. def dot_product_rotation(M, R):
  641. """
  642. Compute the rotated dot product matrix between two datasets.
  644. This function takes two matrices, `M` and `R`, and computes the roation matrix
  645. necessary to align them using their dot product.
  647. Parameters
  648. ----------
  649. M : numpy.ndarray
  650. A 2D array of shape `(n_samples, n_features)`. The matrix to be rotated.
  651. R : numpy.ndarray
  652. A 2D array of shape `(n_samples, n_features)`. The target reference matrix.
  654. Returns
  655. -------
  656. numpy.ndarray
  657. A 2D array of shape `(n_features, n_features)`. The rotation matrix from M to R.
  659. Notes
  660. -----
  661. - The input matrices `M` and `R` must have the same shape.
  662. - This function normalizes each row (sample) to have unit norm, ensuring that the dot
  663. product represents cosine similarity between the corresponding features.
  665. Raises
  666. ------
  667. ValueError
  668. If the shapes of `M` and `R` do not match, or if any row contains only zeros.
  669. """
  671. if M.shape != R.shape:
  672. raise ValueError("Input matrices M and R must have the same shape.")
  673. if np.any(np.linalg.norm(M, axis=1) == 0) or np.any(np.linalg.norm(R, axis=1) == 0):
  674. raise ValueError("Rows of M and R must not be all zeros.")
  676. R -= R.mean(axis=1, keepdims=True)
  677. M -= M.mean(axis=1, keepdims=True)
  678. R /= np.linalg.norm(R, axis=1, keepdims=True)
  679. M /= np.linalg.norm(M, axis=1, keepdims=True)
  680. rotation_mat = np.dot(M.T, R)
  681. return rotation_mat
  685. def procrustes_rotation(M, R):
  686. """
  687. Compute the rotation matrix and scaling factor to align a matric M to a reference R
  688. solving the orthogonal Procrustes problem. This is essentially a wrapper over
  689. scipy.linalg.orthogonal_procrustes() that includes the normalization step (see
  690. scipy.spatial.procrustes()).
  692. Parameters
  693. ----------
  694. M : numpy.ndarray
  695. A 2D array of shape `(n_samples, n_features)`. The matrix to be rotated.
  696. R : numpy.ndarray
  697. A 2D array of shape `(n_samples, n_features)`. The target reference matrix.
  699. Returns
  700. -------
  701. R : (N, N) ndarray
  702. A 2D array of shape `(n_features, n_features)`. The rotation matrix from M to R.
  703. scale : float
  704. Sum of the singular values of ``A.conj().T @ B``.
  706. Raises
  707. ------
  708. ValueError
  709. If the shapes of `M` and `R` do not match, or if any row contains only zeros.
  710. """
  712. if M.shape != R.shape:
  713. raise ValueError("Input matrices M and R must have the same shape.")
  714. if np.any(np.linalg.norm(M, axis=1) == 0) or np.any(np.linalg.norm(R, axis=1) == 0):
  715. raise ValueError("Rows of M and R must not be all zeros.")
  717. R -= R.mean(axis=0)
  718. R /= np.linalg.norm(R)
  719. M -= M.mean(axis=0)
  720. M /= np.linalg.norm(M)
  722. return orthogonal_procrustes(M, R)
  725. def argsort_axes(M, R):
  726. """
  727. Determine optimal axis reordering of a subject embedding to match a reference embedding.
  729. This function computes the absolute Pearson correlation between each pair of axes
  730. in the subject embedding `M` and reference embedding `R`, and solves the optimal
  731. one-to-one axis assignment using the Hungarian algorithm to maximize total correlation.
  733. Parameters
  734. ----------
  735. M : numpy.ndarray
  736. A 2D array of shape `(n_samples, n_features)`. The embedding
  737. dimensions to reorder.
  738. R : numpy.ndarray
  739. A 2D array of shape `(n_samples, n_features)`. The target
  740. reference embedding.
  742. Returns
  743. -------
  744. reorder_idx : ndarray of shape (n_dims,)
  745. Indices that reorder the axes of `M` to best match the axes of `R`
  746. in correlation magnitude. Use as: `M[:, reorder_idx]`.
  748. Notes
  749. -----
  750. - Sign alignment (i.e., flipping axes) is not handled here.
  751. - Axis correspondence is determined using absolute Pearson correlation.
  753. See Also
  754. --------
  755. vector_wise_corr : Computes vector-wise Pearson correlation.
  756. """
  758. ndims = M.shape[1]
  760. if R.shape[1] != ndims:
  761. raise ValueError("Both embeddings must have the same number of dimensions")
  763. S = np.array([vector_wise_corr(R, m.reshape(-1,1)) for m in M.T])
  765. return linear_sum_assignment(abs(S), maximize=True)[1]

embed_utils.py at commit c1a6f0b, under MIT · at the source

Overview

Authors: Francesco Alberti1,2, Pierre-Louis Bazin3, R. Austin Benn1,2, Robert Scholz1,2,4,5, Wei Wei1,2, Alexander Holmes1,2, Victoria Shevchenko1,2, Ulysse Klatzmann1,2, Carla Pallavicini1,2,6,7, Robert Leech8, Daniel S. Margulies1,2
ORCID iDs: Robert Scholz, Victoria Shevchenko, Robert Leech
  1. Université Paris Cité, INCC UMR 8002, CNRS, Paris, France
  2. Centre for Integrative Neuroimaging (OxCIN), FMRIB, Nuffield Department of Clinical Neurosciences, University of Oxford
  3. Full Brain Picture Analytics, Leiden, Netherlands
  4. Wilhelm Wundt Institute for Psychology, Leipzig University, Leipzig, Germany
  5. Max Planck School of Cognition, Leipzig, Germany
  6. National Scientific and Technical Research Council (CONICET), Argentina
  7. Cognitive Neuroscience Center, University of San Andres, Buenos Aires, Argentina
  8. Department of Neuroimaging, King’s College London, UK
Institutions: Université Paris Cité (France); Centre National de la Recherche Scientifique (France); Wellcome Centre for Integrative Neuroimaging (United Kingdom); University of Oxford (United Kingdom); Leipzig University (Germany); Consejo Nacional de Investigaciones Científicas y Técnicas (Argentina); University of San Andrés (Argentina); King's College London (United Kingdom)
Journal: Research Square (preprint)
Dates: published online 8 April 2026
Type: Preprint · Language: English
License: CC BY
Identifiers: DOI 10.21203/rs.3.rs-9200088/v1 · OpenAlex W7151357601
Open access: green, a free copy (OpenAlex)
Status: code verified
Methods: Connectivity, Statistics, fMRI & imaging
Keywords: xx
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institute for Health Research (NIHR) (NIHR203316); Wellcome Trust (203141/Z/16/Z)
Citations: not cited yet (Europe PMC); 39 references in the paper

Abstract

It is widely accepted that the geometry of the cerebral cortex constrains its functional topography. However, how geometric properties contribute to individual differences in cortical organization has not been fully characterized. Here, we investigate whether cortical function varies with cortical geometry across individuals at the local or global scale. We characterize functional organization using the first three gradients of functional connectivity, and project individual surfaces into a shared embedding that captures their intrinsic geometry. Fitting localized spatial models within this embedding via lattice Kriging, we test whether interindividual variation of the gradients is linked to differences in spatial location. These models capture a common spatial structure underlying local gradient transitions across individuals, but do not capture differences in global gradient layout. This suggests that universal geometric properties shape the functional transitions between otherwise stable functional systems, meaningfully contributing to subject-specific functional topography.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

alberti-f/VarioGrad

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: c1a6f0b1ba34caeafac778cb97a6d671fbf9c7c9, 23 June 2026
Languages: Python (22), R (11), Jupyter (3)
Size: 43 files, 36 scripts
Software Heritage: not archived
Found in: “Code Availability”
Holds: README, license file, environment (pyproject.toml, requirements.txt), 3 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (21 files), NiBabel (8 files), scikit-learn (8 files), pandas (7 files), SciPy (7 files), Matplotlib (4 files), Connectome Workbench (4 files), reticulate (2 files), seaborn (2 files), statsmodels (2 files), caret (1 file), h5py (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
38 files
At the source: github.com/alberti-f/VarioGrad

Code Availability

The HCP preprocessing pipelines can be found at https://github.com/Washington-University/HCPpipelines, while the code used to run all subsequent analyses is made available as an installable package on F.A.'s GitHub repository https://github.com/alberti-f/VarioGrad.

Reproduced under the paper's license (CC BY), from the paper cited above.

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:

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

Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. All data are obtainable from the HCP website (https://db.humanconnectome.org/).

Reproduced under the paper's license (CC BY), from the paper cited above.

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, dates, 11 authors, 1 keyword, 2 funders, 37 references.

Cite

This paper

Alberti, F., Bazin, P.-L., Benn, R. A., Scholz, R., Wei, W., Holmes, A., Shevchenko, V., Klatzmann, U., Pallavicini, C., Leech, R., & Margulies, D. S. (2026). Intrinsic cortical geometry is associated with individual differences in local functional organization. Research Square (preprint). https://doi.org/10.21203/rs.3.rs-9200088/v1

BibTeX

@article{alberti2026intrinsic,
author = {Alberti, Francesco and Bazin, Pierre-Louis and Benn, R. Austin and Scholz, Robert and Wei, Wei and Holmes, Alexander and Shevchenko, Victoria and Klatzmann, Ulysse and Pallavicini, Carla and Leech, Robert and Margulies, Daniel S.},
title = {{Intrinsic cortical geometry is associated with individual differences in local functional organization}},
journal = {Research Square (preprint)},
year = {2026},
month = apr,
publisher = {Research Square},
issn = {2693-5015},
doi = {10.21203/rs.3.rs-9200088/v1},
url = {https://doi.org/10.21203/rs.3.rs-9200088/v1}
}

RIS

TY - JOUR
AU - Alberti, Francesco
AU - Bazin, Pierre-Louis
AU - Benn, R. Austin
AU - Scholz, Robert
AU - Wei, Wei
AU - Holmes, Alexander
AU - Shevchenko, Victoria
AU - Klatzmann, Ulysse
AU - Pallavicini, Carla
AU - Leech, Robert
AU - Margulies, Daniel S.
TI - Intrinsic cortical geometry is associated with individual differences in local functional organization
T2 - Research Square (preprint)
J2 - Res Sq
PY - 2026
DA - 2026/04/08
SN - 2693-5015
PB - Research Square
DO - 10.21203/rs.3.rs-9200088/v1
UR - https://doi.org/10.21203/rs.3.rs-9200088/v1
LA - en
ER -

CSL-JSON

{
"id": "10.21203/rs.3.rs-9200088/v1",
"type": "article",
"title": "Intrinsic cortical geometry is associated with individual differences in local functional organization",
"container-title": "Research Square (preprint)",
"author": [
{
"family": "Alberti",
"given": "Francesco"
},
{
"family": "Bazin",
"given": "Pierre-Louis"
},
{
"family": "Benn",
"given": "R. Austin"
},
{
"family": "Scholz",
"given": "Robert"
},
{
"family": "Wei",
"given": "Wei"
},
{
"family": "Holmes",
"given": "Alexander"
},
{
"family": "Shevchenko",
"given": "Victoria"
},
{
"family": "Klatzmann",
"given": "Ulysse"
},
{
"family": "Pallavicini",
"given": "Carla"
},
{
"family": "Leech",
"given": "Robert"
},
{
"family": "Margulies",
"given": "Daniel S."
}
],
"container-title-short": "Res Sq",
"DOI": "10.21203/rs.3.rs-9200088/v1",
"ISSN": "2693-5015",
"publisher": "Research Square",
"URL": "https://doi.org/10.21203/rs.3.rs-9200088/v1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
8
]
]
}
}

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-71270-w [code]
Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan.
Journal: Nature communications
In common: Connectome Workbench, h5py, NiBabel, 6 other tools, 5 references
[2] doi:10.1038/s41467-026-76011-7 [code]
Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: Connectome Workbench, NiBabel, SciPy, 2 other tools, 7 references
[3] doi:10.1038/s41467-026-72931-6 [code]
Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.
Journal: Nature communications
In common: Connectome Workbench, NiBabel, statsmodels, 6 other tools, 4 references
[4] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: reticulate, h5py, NiBabel, 5 other tools, 4 references
[5] doi:10.21203/rs.3.rs-9326213/v1 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: Connectome Workbench, h5py, NiBabel, 6 other tools, 4 references
[6] doi:10.64898/2026.03.09.710558 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: Connectome Workbench, h5py, NiBabel, 6 other tools, 4 references
[7] doi:10.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: h5py, NiBabel, statsmodels, 6 other tools, 4 references
[8] doi:10.1038/s41467-026-76452-0 [code]
Music evokes shared neural representations of imagined narratives across sensory modalities.
Journal: Nature communications
In common: Connectome Workbench, h5py, NiBabel, 7 other tools, 1 reference
[9] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Connectome Workbench, NiBabel, scikit-learn, 3 other tools, 4 references
[10] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: Connectome Workbench, NiBabel, seaborn, 5 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.