OSCR

A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks.

Code ↔ Paper

10 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 10 matches
  1. [1] § Results ↔ langlocfc/plot_template_matching.py, lines 61–106 · score 0.85 · template matching, PM PPr, CG OP, dATN, LanA, PMN
  2. [2] § Results ↔ parcellate/model.py, lines 235–308 · score 0.82 · linear sum assignment, CG OP, template matching, reference atlases, PMN, SAL
  3. [3] § Results › LangFC shows expected topography, distinct from other brain networks ↔ langlocfc/plot_oracle.py, lines 30–107 · score 0.74 · PM PPr, CG OP, dATN, LanA, PMN, SAL
  4. [4] § Results › LangFC shows expected topography, distinct from other brain networks ↔ langlocfc/plot_rest_v_task.py, lines 61–106 · score 0.73 · PM PPr, CG OP, dATN, LanA, PMN, SAL
  5. [5] § Methods › Additional analyses › LangFC periphery (Fig. 2C) ↔ parcellate/model.py, lines 548–606 · score 0.59 · Poisson binomial distribution, space, subnetworks, probabilistic, atlas, parcellation
  6. [6] § Results › LangFC differs in function from nearby networks ↔ langlocfc/plot_performance.py, lines 22–81 · score 0.58 · working memory, MD, auditory, scrambled, intact, listening
  7. [7] § Methods › Parcellation and evaluation procedure › Alignment ↔ parcellate/data.py, lines 242–325 · score 0.56 · linear sum assignment, greedy, alignment, error, networks, parcellations
  8. [8] § Methods › Parcellation and evaluation procedure › Sampling ↔ parcellate/model.py, lines 29–111 · score 0.54 · connectivity matrix, voxelwise, downsampling, detrended, components, binarized
  9. [9] § Methods › fMRI data acquisition and preprocessing ↔ parcellate/data.py, lines 123–181 · score 0.53 · MNI space, FWHM, NIFTI, smoothed, template, resampled
  10. [10] § Methods › Parcellation and evaluation procedure › Alignment ↔ parcellate/data.py, lines 189–222 · score 0.50 · linear sum assignment, greedy, alignment, parcellations

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 · 1,728 lines · 72 KB · no license · 2 matches

  1. import shutil
  2. import time
  3. import datetime
  4. import inspect
  5. import yaml
  6. import numpy as np
  7. import pandas as pd
  8. from scipy import optimize
  9. from sklearn.cluster import MiniBatchKMeans
  10. from sklearn.decomposition import PCA, FastICA
  11. from sklearn.metrics import mutual_info_score, normalized_mutual_info_score, adjusted_rand_score
  12. # from sica.base import StabilizedICA
  13. from nilearn import image, masking, maskers, datasets
  14. import nitransforms as nt
  15. from fast_poibin import PoiBin
  17. from parcellate.cfg import *
  18. from parcellate.data import *
  19. from parcellate.util import *
  22. ######################################
  23. #
  24. # CORE METHODS
  25. #
  26. ######################################
  29. def sample(
  30. output_dir,
  31. functional_paths,
  32. n_networks=100,
  33. fwhm=None,
  34. sample_id=None,
  35. xfm_path=None,
  36. mask_path=None,
  37. detrend=True,
  38. standardize=True,
  39. envelope=False,
  40. independent_runs=False,
  41. data_fraction=1,
  42. tr=2,
  43. low_pass=0.1,
  44. high_pass=0.01,
  45. n_samples=None,
  46. n_components_pca='auto',
  47. n_components_ica='auto',
  48. cluster=True,
  49. target_affine=None,
  50. downsample_to=None,
  51. use_connectivity_profile=False,
  52. use_connectivity_to_regions=True,
  53. binarize_connectivity=True,
  54. transform_connectivity=False,
  55. clustering_kwargs=None,
  56. compress_outputs=True,
  57. dump_kwargs=True,
  58. indent=0,
  59. **kwargs
  60. ):
  61. """
  62. Sample parcellations by clustering the voxel timecourses
  64. :param output_dir: ``str``; Output directory
  65. :param functional_paths: ``list``; Paths to functional images
  66. :param n_networks: ``int``; Number of networks to sample
  67. :param fwhm: ``float`` or ``None``; Full-width at half-maximum for spatial smoothing. If ``None``, no smoothing is
  68. applied.
  69. :param sample_id: ``str``; Sample ID
  70. :param xfm_path: ``str`` or ``None``; if the parcellation is not in MNI space, path to
  71. transformation from MNI to the parcellation space (e.g., native), which will be applied to atlases if
  72. applicable (e.g., Schaeffer et al 18 atlases if `use_connectivity_to_regions` is `True`).
  73. If ``None``, parcellation is assumed to be in MNI space.
  74. :param mask_path: ``str`` or ``None``; Path to mask image. If ``None``, a mask will be computed from the first
  75. functional image.
  76. :param detrend: ``bool``; Whether to detrend the timecourses
  77. :param standardize: ``bool``; Whether to standardize the timecourses
  78. :param envelope: ``bool``; Whether to use the envelope of the timecourses
  79. :param independent_runs: ``bool``; Whether to treat each run as an independent sample
  80. :param data_fraction: ``float``; Fraction of data to use. If < 1, randomly subsample the data.
  81. :param tr: ``float``; Repetition time
  82. :param low_pass: ``float`` or ``None``; Low-pass filter cutoff frequency. If ``None``, no low-pass filtering is
  83. applied.
  84. :param high_pass: ``float`` or ``None``; High-pass filter cutoff frequency. If ``None``, no high-pass filtering is
  85. applied.
  86. :param n_samples: ``int``; Number of samples to draw
  87. :param n_components_pca: ``int`` or ``str``; Number of components to retain after PCA. If ``'auto'``, no PCA is
  88. applied.
  89. :param n_components_ica: ``int`` or ``None``; Number of components to retain after ICA. If ``None``, no ICA is
  90. applied.
  91. :param cluster: ``bool``; Whether to cluster the timecourses
  92. :param target_affine: ``list`` or ``None``; Target affine for spatial resampling. If ``None``, no resampling is
  93. applied.
  94. :param downsample_to: ``int`` or ``None``; If given, downsample the columns of the connectivity matrix to this
  95. voxel size (in mm). If ``None``, no downsampling is applied. Ignored if ``use_connectivity_to_regions`` is
  96. ``True``.
  97. :param use_connectivity_profile: ``bool``; Whether to use the connectivity profile. If ``False``, the raw
  98. timecourses are used.
  99. :param use_connectivity_to_regions: ``bool``; Whether to use the connectivity to regions. If ``False``, the
  100. voxelwise connectivity is used.
  101. :param binarize_connectivity: ``bool``; Whether to binarize the connectivity matrix
  102. :param transform_connectivity: ``bool``; Whether to transform the connectivity matrix
  103. :param clustering_kwargs: ``dict`` or ``None``; Keyword arguments to pass to the clustering algorithm
  104. :param compress_outputs: ``bool``; Whether to compress the output files
  105. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  106. :param indent: ``int``; Indentation level for progress reporting
  107. :param kwargs: ``dict``; Unused keyword arguments
  108. :return: ``None``
  109. """
  110. if len(kwargs):
  111. stderr('WARNING: Unused keyword arguments to `sample()`: %s\n' % ', '.join(kwargs.keys()))
  113. assert isinstance(sample_id, str), 'sample_id must be given as a str'
  114. assert 0 <= data_fraction <= 1, 'data_fraction must be a proportion between 0 and 1'
  116. t0 = time.time()
  118. stderr('%sSampling (sample_id=%s, n_networks=%d)\n' % (' ' * (indent * 2), sample_id, n_networks))
  119. indent += 1
  121. assert isinstance(output_dir, str), 'output_dir must be provided'
  123. sample_dir = get_path(output_dir, 'subdir', 'sample', sample_id)
  124. if not os.path.exists(sample_dir):
  125. os.makedirs(sample_dir)
  127. if clustering_kwargs is None:
  128. clustering_kwargs = dict(
  129. n_init=N_INIT,
  130. init_size=INIT_SIZE
  131. )
  133. if n_samples is None:
  134. n_samples = 0
  136. if dump_kwargs:
  137. kwargs = dict(
  138. output_dir=output_dir,
  139. functional_paths=functional_paths,
  140. n_networks=n_networks,
  141. fwhm=fwhm,
  142. sample_id=sample_id,
  143. xfm_path=xfm_path,
  144. mask_path=mask_path,
  145. detrend=detrend,
  146. standardize=standardize,
  147. envelope=envelope,
  148. independent_runs=independent_runs,
  149. data_fraction=data_fraction,
  150. tr=tr,
  151. low_pass=low_pass,
  152. high_pass=high_pass,
  153. n_samples=n_samples,
  154. n_components_pca=n_components_pca,
  155. n_components_ica=n_components_ica,
  156. cluster=cluster,
  157. target_affine=target_affine,
  158. use_connectivity_profile=use_connectivity_profile,
  159. use_connectivity_to_regions=use_connectivity_to_regions,
  160. binarize_connectivity=binarize_connectivity,
  161. transform_connectivity=transform_connectivity,
  162. clustering_kwargs=clustering_kwargs,
  163. compress_outputs=compress_outputs
  164. )
  165. kwargs_path = get_path(output_dir, 'kwargs', 'sample', sample_id)
  166. with open(kwargs_path, 'w') as f:
  167. yaml.safe_dump(kwargs, f, sort_keys=False)
  168. output_path = get_path(output_dir, 'output', 'sample', sample_id, compressed=compress_outputs)
  170. t1 = time.time()
  171. stderr('%sLoading timecourses' % (' ' * (indent * 2)))
  172. input_data = InputData(
  173. functional_paths=functional_paths,
  174. fwhm=fwhm,
  175. mask_path=mask_path,
  176. standardize=standardize,
  177. detrend=detrend,
  178. envelope=envelope,
  179. tr=tr,
  180. low_pass=low_pass,
  181. high_pass=high_pass
  182. )
  183. v = input_data.v
  184. df = pd.DataFrame([dict(n_trs=input_data.n_trs, n_runs=input_data.n_runs)])
  185. metadata_path = get_path(output_dir, 'metadata', 'sample', sample_id)
  186. df.to_csv(metadata_path, index=False)
  187. stderr(' (%0.2fs)\n' % (time.time() - t1))
  189. stderr('%sN voxels: %d\n' % (' ' * ((indent + 1) * 2), v))
  191. n_runs = input_data.n_runs
  192. samples_all = []
  193. scores_all = []
  194. if independent_runs:
  195. timecourses = input_data.functionals
  196. else:
  197. timecourses = [input_data.timecourses]
  198. for i, timecourse in enumerate(timecourses):
  199. # Sample parcellations by clustering the voxel timecourses
  200. if n_networks > 256:
  201. dtype=np.uint16
  202. else:
  203. dtype=np.uint8
  204. scores = np.zeros(n_samples) # Shape: <n_samples>
  206. X = timecourse
  207. t = X.shape[-1]
  208. X_img = input_data.nii_ref
  209. X_mask = input_data.mask
  210. if target_affine is not None:
  211. stderr('%sSpatial resampling' % (' ' * (indent * 2)))
  212. t1 = time.time()
  213. X_img = input_data.unflatten(X * (1 + 1e-6)) # Hack to force conversion to float
  214. X_img = image.resample_img(
  215. X_img, target_affine=np.diag(np.array(target_affine)), copy_header=True, force_resample=True
  216. )
  217. X_mask = image.new_img_like(input_data.nii_ref, input_data.mask * (1 + 1e-6))
  218. X_mask = image.resample_img(
  219. X_mask, target_affine=np.diag(target_affine), copy_header=True, force_resample=True
  220. )
  221. X_mask = image.get_data(X_mask) > 0.5
  222. X = image.get_data(X_img)[X_mask]
  223. v = X.shape[0]
  224. stderr(' (%0.2fs)\n' % (time.time() - t1))
  225. if n_components_pca:
  226. n_components = n_components_pca
  227. if n_components == 'auto':
  228. n_components = n_networks
  229. stderr('%sPCA transforming (n components = %s)' % (' ' * (indent * 2), n_components))
  230. t1 = time.time()
  231. n_components = min(n_components, t)
  232. m = PCA(n_components=n_components, svd_solver='auto', whiten=True)
  233. X = m.fit_transform(X)
  234. stderr(' (%0.2fs)\n' % (time.time() - t1))
  235. if n_components_ica:
  236. n_components = n_components_ica
  237. if n_components == 'auto':
  238. n_components = n_networks
  239. n_components = min(n_components, X.shape[-1])
  240. stderr('%sICA transforming (n components = %s)' % (' ' * (indent * 2), n_components))
  241. t1 = time.time()
  242. m = FastICA(n_components=n_components, whiten='unit-variance')
  243. X = m.fit_transform(X)
  244. stderr(' (%0.2fs)\n' % (time.time() - t1))
  245. if use_connectivity_profile:
  246. A = standardize_array(X)
  247. if use_connectivity_to_regions:
  248. stderr('%sRetrieving connectivity atlas' % (' ' * (indent * 2)))
  249. t1 = time.time()
  250. B_img = input_data.unflatten(X, mask=X_mask, nii_ref=X_img)
  251. X_mask_img = image.new_img_like(X_img, X_mask > 0.5)
  252. anat_atlas = datasets.fetch_atlas_schaefer_2018(n_rois=1000)
  253. atlas_filename = anat_atlas.maps
  254. atlas_nii = image.smooth_img(atlas_filename, None)
  255. atlas_nii = image.new_img_like(atlas_nii, image.get_data(atlas_nii).astype(np.uint16))
  256. if xfm_path is not None:
  257. xfm = nt.manip.load(xfm_path)
  258. labels = np.unique(image.get_data(atlas_nii))
  259. atlas_nii = image.new_img_like(
  260. atlas_nii,
  261. (image.get_data(atlas_nii)[..., None] == labels[None, None, None, ...]).astype(np.uint16)
  262. )
  263. atlas_nii = nt.resampling.apply(
  264. xfm,
  265. atlas_nii,
  266. input_data.nii_ref,
  267. )
  268. atlas_nii = image.new_img_like(
  269. atlas_nii,
  270. np.argmax(image.get_data(atlas_nii), axis=-1).astype(np.uint16)
  271. )
  272. masker = maskers.NiftiLabelsMasker(labels_img=atlas_nii, mask_img=X_mask_img)
  273. B = standardize_array(masker.fit_transform(B_img).T)
  274. stderr(' (%0.2fs)\n' % (time.time() - t1))
  275. elif downsample_to is not None:
  276. stderr('%sDownsampling columns of connectivity matrix to %d mm' % (' ' * (indent * 2), downsample_to))
  277. t1 = time.time()
  278. B = input_data.unflatten(X * (1 + 1e-6)) # Hack to force conversion to float
  279. B = image.resample_img(
  280. B, target_affine=np.diag(np.array([downsample_to] * 3)), copy_header=True, force_resample=True
  281. )
  282. mask_ = image.new_img_like(input_data.nii_ref, input_data.mask * (1 + 1e-6))
  283. mask_ = image.resample_img(
  284. mask_, target_affine=np.diag(np.array([downsample_to] * 3)), copy_header=True, force_resample=True
  285. )
  286. mask_ = image.get_data(mask_) > 0.5
  287. B = standardize_array(image.get_data(B)[mask_])
  288. stderr(' (%0.2fs)\n' % (time.time() - t1))
  289. else:
  290. B = A
  291. stderr('%sComputing connectivity matrix' % (' ' * (indent * 2)))
  292. t1 = time.time()
  293. X = np.dot(
  294. A,
  295. B.T
  296. )
  297. if binarize_connectivity:
  298. X = (X > np.quantile(X, 0.9)).astype(int)
  299. stderr(' (%0.2fs)\n' % (time.time() - t1))
  300. if transform_connectivity:
  301. stderr('%sTransforming connectivity matrix\n' % (' ' * (indent * 2)))
  302. if n_components_pca:
  303. n_components = n_components_pca
  304. if n_components == 'auto':
  305. n_components = n_networks
  306. stderr('%sPCA transforming (n components = %s)' % (' ' * (indent * 2), n_components))
  307. t1 = time.time()
  308. n_components = min(n_components, t)
  309. m = PCA(n_components=n_components, svd_solver='auto', whiten=True)
  310. X = m.fit_transform(X)
  311. stderr(' (%0.2fs)\n' % (time.time() - t1))
  312. if n_components_ica:
  313. n_components = n_components_ica
  314. if n_components == 'auto':
  315. n_components = n_networks
  316. n_components = min(n_components, X.shape[-1])
  317. stderr('%sICA transforming (n components = %s)\n' % (' ' * (indent * 2), n_components))
  318. t1 = time.time()
  319. m = FastICA(n_components=n_components, whiten='unit-variance')
  320. X = m.fit_transform(X)
  321. stderr(' (%0.2fs)\n' % (time.time() - t1))
  322. if n_samples:
  323. stderr('%sDrawing samples\n' % (' ' * (indent * 2)))
  324. indent += 1
  325. samples = np.zeros((v, n_samples), dtype=dtype) # Shape: <n_voxels, n_samples>
  326. for j in range(n_samples):
  327. if len(timecourses) > 1:
  328. suffix = ' for run %d/%d' % (i + 1, n_runs)
  329. else:
  330. suffix = ''
  331. if n_samples > 1:
  332. stderr('\r%sSample %d/%d%s' % (' ' * (indent * 2), j + 1, n_samples, suffix))
  333. if cluster:
  334. m = MiniBatchKMeans(n_clusters=n_networks, **clustering_kwargs)
  335. _sample = m.fit_predict(X)
  336. _score = m.inertia_
  337. samples[:, j] = _sample
  338. else:
  339. X_ = X
  340. n_components = n_networks
  341. m = FastICA(n_components=n_components, whiten='unit-variance')
  342. X = m.fit_transform(X_)
  343. # Minmax normalize
  344. _sample = X[..., :n_networks]
  345. _sample = np.clip(_sample, 0, np.inf)
  346. _sample = _sample / _sample.max(axis=0, keepdims=True)
  347. _score = 0
  348. if j == 0:
  349. samples = _sample
  350. else:
  351. R = np.dot(standardize_array(samples, axis=0).T, standardize_array(_sample, axis=0))
  352. ix_r, ix_c = optimize.linear_sum_assignment(R, maximize=True)
  353. _sample = _sample[:, ix_c]
  354. samples = (samples * j + _sample) / (j + 1)
  356. scores[j] = _score
  358. if n_samples > 1:
  359. stderr('\n')
  360. indent -= 1
  361. else:
  362. assert n_components_pca or n_components_ica, 'Must use PCA or ICA if not sampling'
  363. assert X.shape == (v, n_networks), 'X must have shape (%s, %s) if not sampling. ' \
  364. 'Got shape %s. Check to make sure that there is ' \
  365. 'at least one matrix decomposition (PCA or ICA) and that the final ' \
  366. 'decomposition has n_components equal to n_networks.' % \
  367. (v, n_networks, str(X.shape))
  369. # Assume a network covers < half the mask volume, flip sign accordingly
  370. X = np.where(np.median(X, axis=0, keepdims=True) > 0, -X, X)
  371. # Clip and normalize (scale is arbitrary)
  372. uq = np.quantile(X, 0.99, axis=0, keepdims=True)
  373. X = np.clip(X, 0, uq) / uq
  374. samples = X.astype(np.float32)
  375. scores = np.zeros((1,))
  377. if target_affine is not None:
  378. samples = input_data.unflatten(samples, mask=X_mask, nii_ref=X_img)
  379. samples = image.resample_to_img(
  380. samples, input_data.nii_ref, interpolation='nearest', copy_header=True, force_resample=True
  381. )
  382. samples = input_data.flatten(samples)
  383. samples_all.append(samples)
  384. scores_all.append(scores)
  386. samples = np.concatenate(samples_all, axis=-1)
  387. samples = input_data.unflatten(samples)
  388. samples.to_filename(output_path)
  389. scores = pd.DataFrame({'sample_score': np.concatenate(scores_all, axis=0)})
  390. evaluation_path = get_path(output_dir, 'evaluation', 'sample', sample_id, compressed=compress_outputs)
  391. scores.to_csv(evaluation_path, index=False)
  393. stderr('%sSampling time: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  396. def align(
  397. output_dir,
  398. alignment_id=None,
  399. sample_id=None,
  400. mask_path=None,
  401. n_alignments=None,
  402. top_k=None,
  403. sort_by_mi=False,
  404. weight_samples=False,
  405. prealign=False,
  406. scoring_method='corr',
  407. minmax_normalize=True,
  408. compress_outputs=True,
  409. dump_kwargs=True,
  410. indent=0,
  411. **kwargs
  412. ):
  413. """
  414. Align sampled clusterings into a probabilistic parcellation.
  416. :param output_dir: ``str``; Output directory
  417. :param alignment_id: ``str``; Alignment ID
  418. :param sample_id: ``str``; Sample ID
  419. :param mask_path: ``str`` or ``None``; Path to mask image. If ``None``, a mask will be computed from the first
  420. functional image.
  421. :param n_alignments: ``int`` or ``None``; Number of alignments to perform. If ``None``, use the number of samples.
  422. :param top_k: ``int`` or ``None``; Number of samples to retain. If ``None``, use all samples.
  423. :param sort_by_mi: ``bool``; Whether to sort samples by mutual information
  424. :param weight_samples: ``bool``; Whether to weight samples by their scores
  425. :param prealign: ``bool``; Whether to prealign the samples using an initial alignment pass
  426. :param scoring_method: ``str``; Scoring method for alignment
  427. :param minmax_normalize: ``bool``; Whether to minmax normalize the aligned networks
  428. :param compress_outputs: ``bool``; Whether to compress the output files
  429. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  430. :param indent: ``int``; Indentation level for progress reporting
  431. :param kwargs: ``dict``; Unused keyword arguments
  432. :return: ``None``
  433. """
  435. if len(kwargs):
  436. stderr('WARNING: Unused keyword arguments to `align()`: %s\n' % ', '.join(kwargs.keys()))
  438. assert isinstance(alignment_id, str), 'alignment_id must be given as a str'
  439. assert isinstance(sample_id, str), 'sample_id must be given as a str'
  441. t0 = time.time()
  443. stderr('%sAligning (alignment_id=%s)\n' % (' ' * (indent * 2), alignment_id))
  444. indent += 1
  445. scoring_method = scoring_method.lower()
  447. assert isinstance(output_dir, str), 'output_dir must be provided'
  449. alignment_dir = get_path(output_dir, 'subdir', 'align', alignment_id)
  450. if not os.path.exists(alignment_dir):
  451. os.makedirs(alignment_dir)
  452. if dump_kwargs:
  453. kwargs = dict(
  454. alignment_id=alignment_id,
  455. sample_id=sample_id,
  456. mask_path=mask_path,
  457. n_alignments=n_alignments,
  458. top_k=top_k,
  459. sort_by_mi=sort_by_mi,
  460. weight_samples=weight_samples,
  461. prealign=prealign,
  462. scoring_method=scoring_method,
  463. minmax_normalize=minmax_normalize,
  464. compress_outputs=compress_outputs,
  465. output_dir=output_dir,
  466. )
  467. kwargs_path = get_path(output_dir, 'kwargs', 'align', alignment_id)
  468. with open(kwargs_path, 'w') as f:
  469. yaml.safe_dump(kwargs, f, sort_keys=False)
  470. output_path = get_path(output_dir, 'output', 'align', alignment_id, compressed=compress_outputs)
  471. sample_path = get_path(output_dir, 'output', 'sample', sample_id, compressed=compress_outputs)
  472. assert os.path.exists(sample_path), 'Sample file %s not found' % sample_path
  474. data = Data(
  475. nii_ref_path=sample_path,
  476. mask_path=mask_path
  477. )
  478. sample_nii = data.nii_ref
  479. if image.get_data(sample_nii).dtype in (np.uint8, np.uint16):
  480. samples = data.flatten(sample_nii)
  481. samples = samples.T # Shape: <n_samples, n_voxels>, values are integer network indices
  483. # Get sample scores
  484. sample_scores = pd.read_csv(get_path(output_dir, 'evaluation', 'sample', sample_id))['sample_score'].values
  485. sample_scores = minmax_normalize_array(sample_scores) # Lower inertia is better
  487. # Compute sample orders/weights
  488. if sort_by_mi: # By pairwise similarity
  489. # Compute similarity between samples
  490. # MI = np.zeros((samples.shape[0], samples.shape[0]))
  491. # for i in range(samples.shape[0]):
  492. # for j in range(i + 1, samples.shape[0]):
  493. # print(i, j)
  494. # mi = mutual_info_score(samples[i], samples[j])
  495. # # mi = normalized_mutual_info_score(samples[i], samples[j])
  496. # # mi = adjusted_rand_score(samples[i], samples[j])
  497. # MI[i, j] = mi
  498. # MI[j, i] = mi
  499. # MI[np.diag_indices(MI.shape[0])] = 1
  500. # MI_mean = MI.mean(axis=0)
  501. # ix = np.argmax(MI_mean)
  502. # mi = MI[ix]
  503. mi = np.zeros(samples.shape[0])
  504. ix = np.argmin(sample_scores)
  505. for i in range(samples.shape[0]):
  506. mi[i] = mutual_info_score(samples[ix], samples[i])
  507. s_ix = np.argsort(mi)[::-1]
  508. samples = samples[s_ix]
  509. w = mi[s_ix]
  510. else: # By inertia
  511. s_ix = np.argsort(sample_scores)
  512. samples = samples[s_ix]
  513. sample_scores = sample_scores[s_ix]
  514. w = 1 - sample_scores # Flip to upweight lower inertia
  515. w = minmax_normalize_array(w)
  517. # Align to reference
  518. if top_k:
  519. samples = samples[:top_k]
  520. w = w[:top_k]
  522. f_kwargs = dict(
  523. samples=samples,
  524. scoring_method='corr',
  525. n_alignments=n_alignments,
  526. prealign=prealign,
  527. indent=indent + 1
  528. )
  529. if weight_samples:
  530. f_kwargs['w'] = w
  532. parcellation = align_samples(**f_kwargs)
  533. else:
  534. stderr('%sSamples are already continuous, likely due to being a matrix decomposition. '
  535. 'Skipping alignment...\n' % (' ' * (indent * 2)))
  536. parcellation = data.flatten(sample_nii).T
  538. if minmax_normalize:
  539. parcellation = minmax_normalize_array(parcellation)
  541. parcellation = data.unflatten(parcellation.T)
  543. parcellation.to_filename(output_path)
  545. stderr('%sAlignment time: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  548. def label(
  549. output_dir,
  550. reference_atlases='default',
  551. labeling_id=None,
  552. alignment_id=None,
  553. sample_id=None,
  554. xfm_path=None,
  555. mask_path=None,
  556. average_first=True,
  557. scoring_method='corr',
  558. atlas_threshold=None,
  559. max_subnetworks=None,
  560. minmax_normalize=True,
  561. use_poibin=True,
  562. eps=1e-3,
  563. compress_outputs=True,
  564. dump_kwargs=True,
  565. indent=0,
  566. **kwargs
  567. ):
  568. """
  569. Label the parcellation against reference atlases
  571. :param output_dir: ``str``; Output directory
  572. :param reference_atlases: ``str`` or ``list``; Reference atlases
  573. :param labeling_id: ``str``; Labeling ID
  574. :param alignment_id: ``str``; Alignment ID
  575. :param sample_id: ``str``; Sample ID
  576. :param xfm_path: ``str`` or ``None``; if the parcellation is not in MNI space, path to
  577. transformation from MNI to the parcellation space (e.g., native), which will be applied to atlases.
  578. If ``None``, parcellation is assumed to be in MNI space.
  579. :param mask_path: ``str`` or ``None``; path to GM mask image. If ``None``, compute the mask in MNI space.
  580. :param average_first: ``bool``; Whether to average the samples before labeling
  581. :param scoring_method: ``str``; Scoring method for labeling
  582. :param atlas_threshold: ``float`` or ``None``; Threshold for binarizing the atlas. If ``None``, no binarization is
  583. applied.
  584. :param max_subnetworks: ``int`` or ``None``; Maximum number of subnetworks to retain. If ``None``, no bound on
  585. the number of subnetworks.
  586. :param minmax_normalize: ``bool``; Whether to minmax normalize the labeled networks
  587. :param use_poibin: ``bool``; Whether to use the Poisson binomial distribution for computing probabilities in
  588. composite networks
  589. :param eps: ``float``; Epsilon for numerical stability
  590. :param compress_outputs: ``bool``; Whether to compress the output files
  591. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  592. :param indent: ``int``; Indentation level for progress reporting
  593. :param kwargs: ``dict``; Unused keyword arguments
  594. :return: ``None``
  595. """
  597. if len(kwargs):
  598. stderr('WARNING: Unused keyword arguments to `label()`: %s\n' % ', '.join(kwargs.keys()))
  600. assert isinstance(labeling_id, str), 'alignment_id must be given as a str'
  602. t0 = time.time()
  604. stderr('%sLabeling (labeling_id=%s)\n' % (' ' * (indent * 2), labeling_id))
  605. indent += 1
  606. scoring_method = scoring_method.lower()
  608. assert isinstance(output_dir, str), 'output_dir must be provided'
  610. labeling_dir = get_path(output_dir, 'subdir', 'label', labeling_id)
  611. if not os.path.exists(labeling_dir):
  612. os.makedirs(labeling_dir)
  613. if dump_kwargs:
  614. kwargs = dict(
  615. reference_atlases=reference_atlases,
  616. labeling_id=labeling_id,
  617. alignment_id=alignment_id,
  618. sample_id=sample_id,
  619. xfm_path=xfm_path,
  620. mask_path=mask_path,
  621. average_first=average_first,
  622. scoring_method=scoring_method,
  623. atlas_threshold=atlas_threshold,
  624. max_subnetworks=max_subnetworks,
  625. minmax_normalize=minmax_normalize,
  626. use_poibin=use_poibin,
  627. eps=eps,
  628. compress_outputs=compress_outputs,
  629. output_dir=output_dir,
  630. )
  631. kwargs_path = get_path(output_dir, 'kwargs', 'label', labeling_id)
  632. with open(kwargs_path, 'w') as f:
  633. yaml.safe_dump(kwargs, f, sort_keys=False)
  634. output_path = get_path(output_dir, 'output', 'label', labeling_id, compressed=compress_outputs)
  635. if average_first:
  636. assert alignment_id is not None, 'alignment_id must be provided if average_first is True'
  637. input_path = get_path(output_dir, 'output', 'align', alignment_id, compressed=compress_outputs)
  638. assert os.path.exists(input_path), 'Alignment file %s not found' % input_path
  639. else:
  640. assert sample_id is not None, 'sample_id must be provided if average_first is False'
  641. input_path = get_path(output_dir, 'output', 'sample', sample_id, compressed=compress_outputs)
  642. assert os.path.exists(input_path), 'Sample file %s not found' % input_path
  644. input_nii = image.smooth_img(input_path, None)
  646. reference_data = AtlasData(
  647. atlases=reference_atlases,
  648. resampling_target_nii=input_nii,
  649. compress_outputs=compress_outputs,
  650. xfm_path=xfm_path,
  651. mask_path=mask_path
  652. )
  653. reference_atlas_names = reference_data.atlas_names
  654. reference_atlases = reference_data.atlases
  655. v = reference_data.v
  656. reference_data.save_atlases(labeling_dir, prefix=REFERENCE_ATLAS_PREFIX)
  658. input_data = reference_data.flatten(input_nii)
  659. if average_first:
  660. n_networks = input_data.shape[-1]
  661. n_samples = None
  662. else:
  663. n_networks = int(input_data.max() + 1)
  664. n_samples = input_data.shape[-1]
  665. input_data = input_data.T # Shape: <(n_networks | n_samples), v>, values are integer network indices
  666. if not max_subnetworks:
  667. max_subnetworks = n_networks
  669. # Find candidate network(s) for each reference
  670. n_reference_atlases = len(reference_atlas_names)
  671. reference_atlas_scores = np.full((n_reference_atlases,), -np.inf)
  672. candidates = {}
  673. results = []
  674. indent += 1
  675. for j, reference_atlas_name in enumerate(reference_atlas_names):
  676. stderr('%sAtlas: %s\n' % (' ' * (indent * 2), reference_atlas_name))
  677. reference_atlas = reference_atlases[reference_atlas_name]
  678. if atlas_threshold is not None:
  679. _reference_atlas = binarize_array(reference_atlas, threshold=atlas_threshold)
  680. else:
  681. _reference_atlas = reference_atlas
  682. if average_first:
  683. scores = np.zeros(n_networks)
  684. for ni in range(n_networks):
  685. if scoring_method == 'corr':
  686. _score = np.corrcoef(input_data[ni], _reference_atlas)[0, 1]
  687. elif scoring_method == 'avg':
  688. _score = np.dot(input_data[ni], _reference_atlas) / input_data[ni].sum()
  689. else:
  690. raise ValueError('Unrecognized scoring method %s.' % scoring_method)
  691. scores[ni] = _score
  692. reference_ix = np.argsort(scores, axis=-1)[::-1]
  693. else:
  694. indent += 1
  695. stderr('%sDirectly aligning samples to reference\n' % (' ' * (indent * 2)))
  696. samples_relabeled = np.zeros_like(input_data)
  697. scores = np.zeros((n_samples, n_networks))
  698. if scoring_method == 'corr':
  699. _reference_atlas = standardize_array(_reference_atlas)
  700. indent += 1
  701. for si in range(n_samples):
  702. stderr('\r%sSample %d/%d' % (' ' * (indent * 2), si + 1, n_samples))
  703. networks = (input_data[si][None, ...] == np.arange(n_networks)[..., None]).astype(float)
  704. if scoring_method == 'corr':
  705. _networks = standardize_array(networks)
  706. _scores = np.dot(
  707. _networks,
  708. _reference_atlas.T
  709. ) / v
  710. else:
  711. num = np.dot(
  712. networks,
  713. _reference_atlas.T
  714. )
  715. denom = networks.sum(axis=-1)
  716. denom[np.where(denom == 0)] = 1
  717. _scores = num / denom
  718. sort_ix = np.argsort(_scores)[::-1]
  719. ranks = np.argsort(sort_ix)
  720. scores[si] = _scores[sort_ix]
  721. samples_relabeled[si] = ranks[input_data[si]]
  722. input_data = samples_relabeled
  723. reference_ix = np.arange(n_networks)
  724. stderr('\n')
  725. indent -= 2
  726. candidate = None
  727. r = -np.inf
  728. r_prev = -np.inf
  729. candidate_list = []
  730. candidate_scores = []
  731. indent += 1
  732. for ni in range(max_subnetworks):
  733. stderr('\r%sSubnetwork %d' % (' ' * (indent * 2), ni + 1))
  734. ix = reference_ix[ni]
  735. if average_first:
  736. _score = scores[ix]
  737. else:
  738. _score = np.tanh(np.arctanh(scores[:, ix] * (1 - 2 * eps) + eps).mean(axis=-1))
  739. candidate_scores.append(_score)
  740. _candidate = candidate
  741. if average_first:
  742. candidate = input_data[ix]
  743. else:
  744. candidate = input_data == ni
  745. weights = scores[:, ni:ni+1]
  746. weights = minmax_normalize_array(weights)
  747. candidate = (candidate * weights).sum(axis=0) / weights.sum()
  748. candidate = np.clip(candidate, 0, 1)
  749. candidate_list.append(candidate)
  750. if use_poibin and ni > 0:
  751. __candidate = np.zeros(v)
  752. for _v in range(v):
  753. p = 1 - PoiBin([c[_v] for c in candidate_list]).cdf[0]
  754. __candidate[_v] = p
  755. candidate = __candidate
  756. else:
  757. candidate = np.stack(candidate_list, axis=-1).sum(axis=-1)
  758. if scoring_method == 'corr':
  759. r = np.corrcoef(candidate, reference_atlas)[0, 1]
  760. elif scoring_method == 'avg':
  761. r = np.dot(candidate, reference_atlas) / candidate.sum()
  762. else:
  763. raise ValueError('Unrecognized scoring method %s.' % scoring_method)
  765. if r <= r_prev:
  766. candidate_list = candidate_list[:-1]
  767. candidate_scores = candidate_scores[:-1]
  768. candidate = _candidate
  769. r = r_prev
  770. break
  771. r_prev = r
  772. stderr('\n')
  773. indent -= 1
  775. if not use_poibin:
  776. candidate = minmax_normalize_array(candidate)
  778. reference_atlas_scores[j] = r
  779. candidate_list.insert(0, candidate)
  780. candidate_scores.insert(0, r)
  782. for s, candidate in enumerate(candidate_list):
  783. row = {
  784. 'parcel': reference_atlas_name if s == 0 else '%s_sub%d' % (reference_atlas_name, s),
  785. '%sname' % REFERENCE_ATLAS_PREFIX: reference_atlas_name,
  786. '%sscore' % REFERENCE_ATLAS_PREFIX: candidate_scores[s]
  787. }
  788. if s == 0:
  789. row['parcel_type'] = 'network'
  790. else:
  791. row['parcel_type'] = 'subnetwork%d' % s
  792. if minmax_normalize:
  793. candidate = minmax_normalize_array(candidate)
  794. candidate_list[s] = candidate
  795. row['n_voxels'] = candidate.sum()
  796. results.append(row)
  797. candidate = reference_data.unflatten(candidate)
  798. if s == 0:
  799. suffix = ''
  800. else:
  801. suffix = '_sub%d' % s
  802. suffix += get_suffix(compress_outputs)
  803. candidate.to_filename(join(labeling_dir, '%s%s' % (reference_atlas_name, suffix)))
  805. candidates[reference_atlas_name] = candidate_list
  807. indent -= 1
  809. results = pd.DataFrame(results)
  810. results.to_csv(output_path, index=False)
  812. stderr('%sLabeling time: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  815. def evaluate(
  816. output_dir,
  817. evaluation_atlases=None,
  818. evaluation_map=None,
  819. evaluation_id=None,
  820. labeling_id=None,
  821. xfm_path=None,
  822. transform_evaluation_atlases=False,
  823. mask_path=None,
  824. network_threshold=None,
  825. compress_outputs=True,
  826. dump_kwargs=True,
  827. indent=0,
  828. **kwargs
  829. ):
  830. """
  831. Evaluate the labeled atlases against task maps (evaluation atlases)
  833. :param output_dir: ``str``; Output directory
  834. :param evaluation_atlases: ``dict`` or ``None``; Map from evaluation atlas names to paths. If ``None``, no
  835. evaluation will be run.
  836. :param evaluation_map: ``dict`` or ``None``; Map from network names to evaluations to perform. If ``None``,
  837. all evaluations will be performed on each network.
  838. :param evaluation_id: ``str``; Evaluation ID
  839. :param labeling_id: ``str``; Labeling ID
  840. :param xfm_path: ``str`` or ``None``; if the parcellation is not in MNI space, path to
  841. transformation from MNI to the parcellation space (e.g., native), which will be applied to atlases.
  842. If ``None``, parcellation is assumed to be in MNI space.
  843. :param transform_evaluation_atlases: ``bool``; Whether to transform the evaluation atlases from MNI space to
  844. the parcellation space. Ignored unless `xfm` is provided. Defaults to ``False``, i.e., assumes that the
  845. evaluation atlases (subject-specific statmaps) are estimated in the same space as the parcellation. Only use
  846. ``True`` if the evaluation atlases are estimated in MNI space.
  847. :param mask_path: ``str`` or ``None``; path to GM mask image. If ``None``, compute the mask in MNI space.
  848. :param network_threshold: ``float`` or ``None``; Threshold for binarizing the networks. If ``None``, no
  849. binarization is applied.
  850. :param compress_outputs: ``bool``; Whether to compress the output files
  851. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  852. :param indent: ``int``; Indentation level for progress reporting
  853. :param kwargs: ``dict``; Unused keyword arguments
  854. :return: ``None``
  855. """
  857. if len(kwargs):
  858. stderr('WARNING: Unused keyword arguments to `evaluate()`: %s\n' % ', '.join(kwargs.keys()))
  860. assert isinstance(evaluation_id, str), 'evaluation_id must be given as a str'
  861. assert isinstance(labeling_id, str), 'alignment_id must be given as a str'
  863. t0 = time.time()
  865. stderr('%sEvaluating (evaluation_id=%s)\n' % (' ' * (indent * 2), evaluation_id))
  866. indent += 1
  868. assert isinstance(output_dir, str), 'output_dir must be provided'
  870. evaluation_dir = get_path(output_dir, 'subdir', 'evaluate', evaluation_id)
  871. if not os.path.exists(evaluation_dir):
  872. os.makedirs(evaluation_dir)
  873. if dump_kwargs:
  874. kwargs = dict(
  875. output_dir=output_dir,
  876. evaluation_atlases=evaluation_atlases,
  877. evaluation_map=evaluation_map,
  878. evaluation_id=evaluation_id,
  879. labeling_id=labeling_id,
  880. xfm_path=xfm_path,
  881. transform_evaluation_atlases=transform_evaluation_atlases,
  882. mask_path=mask_path,
  883. network_threshold=network_threshold,
  884. compress_outputs=compress_outputs,
  885. )
  886. kwargs_path = get_path(output_dir, 'kwargs', 'evaluate', evaluation_id)
  887. with open(kwargs_path, 'w') as f:
  888. yaml.safe_dump(kwargs, f, sort_keys=False)
  889. output_path = get_path(output_dir, 'output', 'evaluate', evaluation_id, compressed=compress_outputs)
  891. suffix = get_suffix(compress_outputs)
  893. if evaluation_atlases is None:
  894. evaluation_atlases = {}
  896. # Collect references atlases and alignments
  897. labeling_dir = get_path(output_dir, 'subdir', 'label', labeling_id)
  898. label_kwargs = get_cfg(get_path(output_dir, 'kwargs', 'label', labeling_id))
  899. if evaluation_map is None:
  900. evaluation_map = {}
  901. for path in os.listdir(labeling_dir):
  902. if path.startswith(REFERENCE_ATLAS_PREFIX):
  903. reference_atlas_name = path[len(REFERENCE_ATLAS_PREFIX):-len(suffix)]
  904. evaluation_map[reference_atlas_name] = list(evaluation_atlases.keys())
  905. else:
  906. for path in os.listdir(labeling_dir):
  907. if path.startswith(REFERENCE_ATLAS_PREFIX):
  908. reference_atlas_name = path[len(REFERENCE_ATLAS_PREFIX):-len(suffix)]
  909. if reference_atlas_name not in evaluation_map:
  910. evaluation_map[reference_atlas_name] = []
  912. reference_atlas_names = label_kwargs['reference_atlases']
  913. if reference_atlas_names is None:
  914. reference_atlas_names = []
  915. elif isinstance(reference_atlas_names, str):
  916. reference_atlas_names = [reference_atlas_names]
  917. _reference_atlas_names = []
  918. for reference_atlas in reference_atlas_names:
  919. if isinstance(reference_atlas, str) and reference_atlas.lower() in ('default', 'all', 'all_reference'):
  920. reference_atlas = ALL_REFERENCE
  921. elif isinstance(reference_atlas, dict):
  922. reference_atlas = list(reference_atlas.keys())
  923. else:
  924. reference_atlas = [reference_atlas]
  925. _reference_atlas_names.extend(reference_atlas)
  926. reference_atlas_names = _reference_atlas_names
  927. reference_atlases = []
  928. candidates = {}
  929. resampling_target_nii = None
  930. for reference_atlas in reference_atlas_names:
  931. reference_atlas_path = join(labeling_dir, '%s%s%s' % (REFERENCE_ATLAS_PREFIX, reference_atlas, suffix))
  932. assert os.path.exists(reference_atlas_path), 'Reference atlas %s not found' % reference_atlas_path
  933. reference_atlases.append({reference_atlas: reference_atlas_path})
  935. for path in os.listdir(labeling_dir):
  936. if path.startswith(reference_atlas) and path.endswith(suffix):
  937. atlas_name = path[:-len(suffix)]
  938. atlas_name = re.sub(r'_sub\d+$', '', atlas_name)
  939. if atlas_name == reference_atlas:
  940. trim = len(suffix)
  941. name = path[:-trim]
  942. path = join(labeling_dir, path)
  943. if reference_atlas not in candidates:
  944. candidates[reference_atlas] = {}
  945. candidates[reference_atlas][name] = get_nii(path, add_to_cache=False)
  946. if network_threshold:
  947. data = image.get_data(candidates[reference_atlas][name])
  948. data = binarize_array(data, threshold=network_threshold)
  949. candidates[reference_atlas][name] = image.new_img_like(candidates[reference_atlas][name], data)
  950. if resampling_target_nii is None:
  951. resampling_target_nii = candidates[reference_atlas][name]
  953. # Format data
  954. reference_data = AtlasData(
  955. atlases=reference_atlases,
  956. resampling_target_nii=resampling_target_nii,
  957. compress_outputs=compress_outputs,
  958. xfm_path=xfm_path,
  959. mask_path=mask_path
  960. )
  961. reference_atlases = reference_data.atlases
  962. evaluation_kwargs = dict(
  963. atlases=evaluation_atlases,
  964. resampling_target_nii=resampling_target_nii,
  965. compress_outputs=compress_outputs,
  966. mask_path=mask_path
  967. )
  968. if xfm_path is not None and transform_evaluation_atlases:
  969. evaluation_kwargs['xfm_path'] = xfm_path
  970. evaluation_data = AtlasData(**evaluation_kwargs)
  971. evaluation_atlases = evaluation_data.atlases
  972. evaluation_data.save_atlases(evaluation_dir, prefix=EVALUATION_ATLAS_PREFIX)
  973. for x in candidates:
  974. for y in candidates[x]:
  975. candidates[x][y] = reference_data.flatten(candidates[x][y])
  977. stderr(' ' * (indent * 2) + 'Results:\n')
  978. results = []
  979. for reference_atlas_name in reference_atlas_names:
  980. # Score reference atlas as if it were a candidate parcellation (baseline)
  981. reference_atlas = reference_atlases[reference_atlas_name]
  982. _evaluation_atlases = {x: evaluation_atlases[x] for x in evaluation_atlases
  983. if x in evaluation_map[reference_atlas_name]}
  984. atlas = reference_atlas
  985. atlas_name = '%s%s' % (REFERENCE_ATLAS_PREFIX, reference_atlas_name)
  986. row = _get_evaluation_row(
  987. atlas,
  988. atlas_name,
  989. reference_atlas_name=reference_atlas_name,
  990. reference_atlases=reference_atlases,
  991. evaluation_atlases=_evaluation_atlases
  992. )
  993. row['parcel_type'] = 'baseline'
  994. results.append(row)
  995. stderr(_pretty_print_evaluation_row(row, indent=indent + 1) + '\n')
  997. # Score evaluation atlases as if they were candidate parcellations (baseline)
  998. for evaluation_atlas_name in _evaluation_atlases:
  999. atlas = _evaluation_atlases[evaluation_atlas_name]
  1000. atlas_name = '%s%s' % (EVALUATION_ATLAS_PREFIX, evaluation_atlas_name)
  1001. row = _get_evaluation_row(
  1002. atlas,
  1003. atlas_name,
  1004. reference_atlas_name=reference_atlas_name,
  1005. reference_atlases=reference_atlases,
  1006. evaluation_atlases=_evaluation_atlases
  1007. )
  1008. row['parcel_type'] = 'baseline'
  1009. results.append(row)
  1011. # Score candidate parcellations
  1012. candidate_names = sorted(
  1013. list(candidates[reference_atlas_name].keys()),
  1014. key=candidate_name_sort_key
  1015. )
  1016. for c, atlas_name in enumerate(candidate_names):
  1017. atlas = candidates[reference_atlas_name][atlas_name]
  1018. row = _get_evaluation_row(
  1019. atlas,
  1020. atlas_name,
  1021. reference_atlas_name=reference_atlas_name,
  1022. reference_atlases=reference_atlases,
  1023. evaluation_atlases=_evaluation_atlases
  1024. )
  1025. if c == 0:
  1026. row['parcel_type'] = 'network'
  1027. else:
  1028. row['parcel_type'] = 'subnetwork%d' % c
  1029. results.append(row)
  1030. if c == 0 or (len(candidate_names) > 2):
  1031. stderr(_pretty_print_evaluation_row(row, indent=indent + 1) + '\n')
  1033. results = pd.DataFrame(results)
  1034. results.to_csv(output_path, index=False)
  1036. stderr('%sEvaluation time: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  1039. def aggregate(
  1040. output_dir,
  1041. action_sequence,
  1042. grid_params,
  1043. aggregation_id=None,
  1044. evaluation_id=None,
  1045. labeling_id=None,
  1046. subnetwork_id=1,
  1047. exclude='LANA',
  1048. kernel_radius=5,
  1049. eps=1e-3,
  1050. compress_outputs=None,
  1051. dump_kwargs=True,
  1052. indent=0,
  1053. **kwargs
  1054. ):
  1055. """
  1056. Aggregate the results of a grid search and select to top setting.
  1058. :param output_dir: ``str``; Output directory
  1059. :param action_sequence: ``list`` of ``dict``; Action sequence
  1060. :param grid_params: ``dict``; Grid parameters
  1061. :param aggregation_id: ``str``; Aggregation ID
  1062. :param evaluation_id: ``str`` or ``None``; Evaluation ID
  1063. :param labeling_id: ``str`` or ``None``; Labeling ID
  1064. :param subnetwork_id: ``int``; Subnetwork ID to use for scoring
  1065. :param exclude: ``str`` or ``list`` of ``str``; Network(s) to exclude from scoring
  1066. :param kernel_radius: ``int``; Kernel radius for smoothing the grid
  1067. :param eps: ``float``; Epsilon for numerical stability
  1068. :param compress_outputs: ``bool`` or ``None``; Whether to compress the output files
  1069. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  1070. :param indent: ``int``; Indentation level for progress reporting
  1071. :param kwargs: ``dict``; Unused keyword arguments
  1072. :return: ``None``
  1073. """
  1074. if len(kwargs):
  1075. stderr('WARNING: Unused keyword arguments to `aggregate()`: %s\n' % ', '.join(kwargs.keys()))
  1077. assert isinstance(output_dir, str), 'output_dir must be given as a str'
  1078. assert isinstance(grid_params, dict), 'grid_params must be given as a dict'
  1080. t0 = time.time()
  1081. stderr('%sAggregating grid\n' % (' ' * (indent * 2)))
  1082. indent += 1
  1084. _labeling_id = get_action_attr('label', action_sequence, 'id')
  1085. if labeling_id is None:
  1086. labeling_id = _labeling_id
  1087. else:
  1088. assert labeling_id == _labeling_id, ('Mismatch between provided labeling_id (%s) '
  1089. 'and the one contained in action_sequence (%s).' % labeling_id, _labeling_id)
  1091. _evaluation_id = get_action_attr('evaluate', action_sequence, 'id')
  1092. if evaluation_id is None:
  1093. evaluation_id = _evaluation_id
  1094. else:
  1095. assert evaluation_id == _evaluation_id, ('Mismatch between provided evaluation_id (%s) '
  1096. 'and the one contained in action_sequence (%s).' % evaluation_id, _evaluation_id)
  1098. _aggregation_id = get_action_attr('aggregate', action_sequence, 'id')
  1099. if aggregation_id is None:
  1100. aggregation_id = _aggregation_id
  1101. else:
  1102. assert aggregation_id == _aggregation_id, ('Mismatch between provided aggregation_id (%s) '
  1103. 'and the one contained in action_sequence (%s).' % aggregation_id, _aggregation_id)
  1105. aggregation_dir = get_path(output_dir, 'subdir', 'aggregate', aggregation_id)
  1106. if not os.path.exists(aggregation_dir):
  1107. os.makedirs(aggregation_dir)
  1108. if dump_kwargs:
  1109. kwargs = dict(
  1110. output_dir=output_dir,
  1111. action_sequence=action_sequence,
  1112. grid_params=grid_params,
  1113. evaluation_id=evaluation_id,
  1114. aggregation_id=aggregation_id,
  1115. labeling_id=labeling_id,
  1116. subnetwork_id=subnetwork_id,
  1117. exclude=exclude,
  1118. kernel_radius=kernel_radius,
  1119. eps=eps,
  1120. compress_outputs=compress_outputs
  1121. )
  1122. kwargs_path = get_path(output_dir, 'kwargs', 'aggregate', aggregation_id)
  1123. with open(kwargs_path, 'w') as f:
  1124. yaml.safe_dump(kwargs, f, sort_keys=False)
  1125. output_path = get_path(output_dir, 'output', 'aggregate', aggregation_id, compressed=compress_outputs)
  1126. evaluation_path = get_path(output_dir, 'evaluation', 'aggregate', aggregation_id)
  1128. grid_settings = get_iterator_from_grid_params(grid_params)
  1129. grid_array, ix2val = get_grid_array_from_grid_params(grid_params)
  1130. val2ix = {x: {y: i for i, y in enumerate(ix2val[x])} for x in ix2val}
  1131. grid_keys = sorted(list(ix2val.keys()))
  1132. results = []
  1133. grid_ids = []
  1134. scores = []
  1135. for grid_setting in grid_settings:
  1136. grid_id = get_grid_id(grid_setting)
  1137. grid_ids.append(grid_id)
  1138. _output_dir = get_path(output_dir, 'subdir', 'grid', grid_id)
  1139. labeling_dir = get_path(_output_dir, 'subdir', 'label', labeling_id)
  1140. if os.path.exists(labeling_dir):
  1141. results_file_path = get_path(_output_dir, 'output', 'label', labeling_id)
  1142. _results = pd.read_csv(results_file_path)
  1143. _results['grid_id'] = grid_id
  1144. results.append(_results)
  1145. score = _get_atlas_score_from_df(_results, subnetwork_id=subnetwork_id, exclude=exclude, eps=eps)
  1146. else:
  1147. raise ValueError(('No available selection criteria for grid_id %s (no evaluation or labeling '
  1148. 'data found). Aggregation failed.' % grid_id))
  1149. scores.append(score)
  1150. ix = tuple([val2ix[x][grid_setting[x]] for x in grid_keys])
  1151. grid_array[ix] = score
  1152. results = pd.concat(results, axis=0)
  1154. # Weighted average
  1155. if kernel_radius > 1:
  1156. grid_array = smooth(grid_array, kernel_radius=kernel_radius)
  1158. # Select configuration
  1159. best_ix = np.unravel_index(np.argmax(grid_array), grid_array.shape)
  1160. best_setting = {x: ix2val[x][best_ix[i]] for i, x in enumerate(grid_keys)}
  1161. best_id = get_grid_id(best_setting)
  1162. best_score = grid_array[best_ix]
  1164. results['selected'] = (results.grid_id == best_id) & (results.parcel_type != 'baseline')
  1165. results.to_csv(evaluation_path, index=False)
  1167. # Save configuration
  1168. best_grid_dir = get_path(output_dir, 'subdir', 'grid', best_id)
  1169. _action_sequence = []
  1170. for action in action_sequence:
  1171. action_type, action_id = action['type'], action['id']
  1172. if action_type != 'aggregate':
  1173. action = copy.deepcopy(action)
  1174. kwargs_path = get_path(best_grid_dir, 'kwargs', action_type, action_id)
  1175. exists = os.path.exists(kwargs_path)
  1176. if action_type in ('sample', 'label', 'aggregate'):
  1177. assert exists, '%s does not exist' % kwargs_path
  1178. kwargs = get_cfg(kwargs_path)
  1179. kwargs.update(dict(
  1180. output_dir=output_dir
  1181. ))
  1182. action['kwargs'] = kwargs
  1183. _action_sequence.append(action)
  1185. parcellate_kwargs = dict(
  1186. output_dir=output_dir,
  1187. action_sequence=_action_sequence,
  1188. grid_params=None,
  1189. eps=eps,
  1190. )
  1192. with open(output_path, 'w') as f:
  1193. yaml.safe_dump(parcellate_kwargs, f, sort_keys=False)
  1195. stderr('%sBest grid_id: %s | atlas score: %0.3f\n' % (' ' * (indent * 2), best_id, best_score))
  1197. stderr('%sAggregation time: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  1200. def parcellate(
  1201. output_dir,
  1202. action_sequence,
  1203. grid_params=None,
  1204. grid_only=False,
  1205. eps=1e-3,
  1206. compress_outputs=True,
  1207. overwrite=False,
  1208. dump_kwargs=True,
  1209. indent=0,
  1210. **kwargs
  1211. ):
  1212. """
  1213. Parcellate the data (run an entire action sequence and/or grid search)
  1215. :param output_dir: ``str``; Output directory
  1216. :param action_sequence: ``list`` of ``dict``; Action sequence
  1217. :param grid_params: ``dict`` or ``None``; Grid parameters
  1218. :param grid_only: ``bool``; Whether to only run the grid search
  1219. :param eps: ``float``; Epsilon for numerical stability
  1220. :param compress_outputs: ``bool``; Whether to compress the output files
  1221. :param overwrite: ``bool`` or ``dict``; Overwrite parameters
  1222. :param dump_kwargs: ``bool``; Whether to dump the keyword arguments to a YAML file
  1223. :param indent: ``int``; Indentation level for progress reporting
  1224. :param kwargs: ``dict``; Unused keyword arguments
  1225. :return: ``None``
  1226. """
  1227. if len(kwargs):
  1228. stderr('WARNING: Unused keyword arguments to `parcellate()`: %s\n' % ', '.join(kwargs.keys()))
  1230. assert isinstance(output_dir, str), 'output_dir is required, must be given as a str'
  1231. assert isinstance(action_sequence, list), ('action_sequence is required, must be given as a list of dict'
  1232. 'and grid_params must be provided as dicts, or neither can be.')
  1233. assert grid_params is not None or not grid_only, ('grid_params must be provided if grid_only is True.')
  1235. validate_action_sequence(action_sequence)
  1237. t0 = time.time()
  1238. stderr('%sParcellating\n' % (' ' * (indent * 2)))
  1239. indent += 1
  1241. sample_id = get_action_attr('sample', action_sequence, 'id')
  1242. labeling_id = get_action_attr('label', action_sequence, 'id')
  1243. aggregation_id = get_action_attr('aggregate', action_sequence, 'id')
  1244. parcellation_id = get_action_attr('parcellate', action_sequence, 'id')
  1246. assert isinstance(sample_id, str), 'sample_id is required, must be given as a str.'
  1247. assert isinstance(labeling_id, str), 'labeling_id is required, must be given as a str.'
  1248. assert isinstance(parcellation_id, str), 'parcellation_id is required, must be given as a str.'
  1250. overwrite = get_overwrite(overwrite)
  1252. use_grid = aggregation_id is not None
  1253. grid_optimized = get_action('parcellate', action_sequence)['kwargs'].get('grid_optimized', True)
  1255. parcellation_dir = get_path(output_dir, 'subdir', 'parcellate', parcellation_id)
  1256. if not os.path.exists(parcellation_dir):
  1257. os.makedirs(parcellation_dir)
  1258. if dump_kwargs:
  1259. kwargs = dict(
  1260. output_dir=output_dir,
  1261. action_sequence=action_sequence,
  1262. grid_params=grid_params,
  1263. grid_optimized=grid_optimized,
  1264. eps=eps,
  1265. compress_outputs=compress_outputs,
  1266. )
  1267. kwargs_path = get_path(output_dir, 'kwargs', 'parcellate', parcellation_id)
  1268. with open(kwargs_path, 'w') as f:
  1269. yaml.safe_dump(kwargs, f, sort_keys=False)
  1270. output_path = get_path(output_dir, 'output', 'parcellate', parcellation_id, compressed=compress_outputs)
  1272. for action in action_sequence:
  1273. action['kwargs']['output_dir'] = output_dir
  1274. action['kwargs']['compress_outputs'] = compress_outputs
  1276. suffix = get_suffix(compress_outputs)
  1278. # Grid search
  1279. if use_grid and grid_params:
  1280. stderr('%sGrid searching\n' % (' ' * (indent * 2)))
  1282. # Core loop
  1283. indent += 1
  1284. grid_settings = get_iterator_from_grid_params(grid_params)
  1285. for grid_setting in grid_settings:
  1286. grid_id = get_grid_id(grid_setting)
  1288. stderr('%sGrid id: %s\n' % (' ' * (indent * 2), grid_id))
  1290. # Update kwargs
  1291. _output_dir = get_path(output_dir, 'subdir', 'grid', grid_id)
  1292. _action_sequence = []
  1293. for action in action_sequence:
  1294. if action['type'] != 'aggregate':
  1295. action = copy.deepcopy(action)
  1296. _kwargs = action['kwargs']
  1297. if action['type'] == 'parcellate':
  1298. # Don't pass down the top-level parcellation kwargs
  1299. _kwargs = {x: _kwargs[x] for x in _kwargs if x in
  1300. ('output_dir', 'compress_outputs', 'parcellation_id')}
  1301. action['kwargs'] = _kwargs
  1302. _kwarg_keys = set(inspect.signature(ACTIONS[action['type']]).parameters.keys())
  1303. _grid_setting = {x: grid_setting[x] for x in grid_setting if x in _kwarg_keys}
  1304. _kwargs.update(_grid_setting)
  1305. _action_sequence.append(action)
  1307. # Recursion bottoms out since grid_params is None
  1308. parcellate(
  1309. _output_dir,
  1310. _action_sequence,
  1311. eps=eps,
  1312. compress_outputs=compress_outputs,
  1313. overwrite=overwrite,
  1314. indent=indent + 1
  1315. )
  1317. indent -= 1
  1319. if grid_only:
  1320. stderr('%sTotal time elapsed: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  1321. return
  1323. # Aggregate
  1324. action = None
  1325. _action_sequence = []
  1326. for action in action_sequence:
  1327. if action['type'] == 'evaluate':
  1328. # Evaluation is not a dependency of aggregate
  1329. continue
  1330. _action_sequence.append(action)
  1331. if action['type'] == 'aggregate':
  1332. break
  1333. assert action is not None and action['type'] == 'aggregate', ('action type "aggregate" not found in '
  1334. 'action_sequence')
  1335. mtime, exists = check_deps(
  1336. output_dir,
  1337. _action_sequence,
  1338. compressed=True
  1339. )
  1340. stale = mtime == 1
  1341. if overwrite['aggregate'] or stale or not exists:
  1342. # Recursion bottoms out since grid_params is None
  1343. _aggregate_kwargs = copy.deepcopy(action['kwargs'])
  1344. # _action_sequence.append(get_action('aggregate', action_sequence))
  1345. _aggregate_kwargs.update(dict(
  1346. output_dir=output_dir,
  1347. action_sequence=action_sequence,
  1348. grid_params=grid_params,
  1349. eps=eps,
  1350. compress_outputs=compress_outputs,
  1351. indent=indent,
  1352. ))
  1353. aggregate(**_aggregate_kwargs)
  1354. else:
  1355. stderr('%sAggregation exists. Skipping. To re-aggregate, run with overwrite=True.\n' %
  1356. (' ' * (indent * 2)))
  1357. aggregation_output_path = get_path(output_dir, 'output', 'aggregate', aggregation_id)
  1359. # Parcellate
  1360. action = None
  1361. _action_sequence = []
  1362. with open(aggregation_output_path, 'r') as f:
  1363. parcellate_kwargs = yaml.safe_load(f)
  1364. parcellate_override_kwargs = get_action('parcellate', action_sequence)
  1365. if parcellate_override_kwargs is not None and 'kwargs' in parcellate_override_kwargs:
  1366. parcellate_override_kwargs = parcellate_override_kwargs['kwargs']
  1367. for action in action_sequence:
  1368. if action['type'] == 'aggregate':
  1369. continue
  1370. action_ = get_action(action['type'], parcellate_kwargs['action_sequence'])
  1371. if grid_optimized:
  1372. kwargs_update = {x: action['kwargs'][x] for x in action['kwargs'] if x not in grid_params}
  1373. else:
  1374. kwargs_update = action['kwargs']
  1375. # Change any kwargs that differ between the config and the saved kwargs file
  1376. if action['type'] != 'parcellate':
  1377. _kwarg_keys = set(inspect.signature(ACTIONS[action['type']]).parameters.keys())
  1378. kwargs_update = {x: kwargs_update[x] for x in kwargs_update if x in _kwarg_keys}
  1379. action_['kwargs'].update(kwargs_update)
  1380. # Change any kwargs that are over-ridden by the settings in the parcellate action
  1381. if action['type'] in parcellate_override_kwargs:
  1382. action_['kwargs'].update(parcellate_override_kwargs[action['type']])
  1383. assert action is not None and action['type'] == 'parcellate', ('Final action type in '
  1384. 'grid searched "action_sequence" must be "parcellate"')
  1385. action_prefix = []
  1386. if grid_optimized:
  1387. for action in action_sequence[:-1]: # Add dependencies to grid, ignoring last ('parcellate') action
  1388. if action['type'] != 'evaluate': # Changing the evaluation doesn't make the aggregation stale
  1389. action_prefix.append(dict(
  1390. type=action['type'],
  1391. id=action['id'],
  1392. kwargs={}
  1393. ))
  1394. parcellate_kwargs['action_sequence'] = action_prefix + parcellate_kwargs['action_sequence']
  1395. parcellate_kwargs['dump_kwargs'] = False # Don't let recursive call overwrite top-level kwargs file
  1396. parcellate_kwargs['overwrite'] = overwrite
  1397. parcellate(**parcellate_kwargs)
  1398. else:
  1399. action_sequence_full = action_sequence
  1400. _action_sequence = []
  1401. for a in range(len(action_sequence) - 1, -1, -1):
  1402. action = action_sequence[a]
  1403. _action_sequence.insert(0, action)
  1404. if action['type'] == 'sample':
  1405. break
  1406. action_sequence = _action_sequence
  1408. sample_id = get_action_attr('sample', action_sequence, 'id')
  1409. alignment_id = get_action_attr('align', action_sequence, 'id')
  1410. labeling_id = get_action_attr('label', action_sequence, 'id')
  1411. evaluation_id = get_action_attr('evaluate', action_sequence, 'id')
  1412. aggregation_id = get_action_attr('aggregate', action_sequence_full, 'id')
  1414. parcellate_kwargs = get_action_attr('parcellate', action_sequence, 'kwargs')
  1416. n = len(action_sequence)
  1417. N = len(action_sequence_full)
  1418. for a, action in enumerate(action_sequence):
  1419. action_type = action['type']
  1420. action_kwargs = action['kwargs']
  1422. e = N - n + a + 1
  1424. mtime, exists = check_deps(
  1425. output_dir,
  1426. action_sequence_full[:e],
  1427. compressed=True
  1428. )
  1429. stale = mtime == 1
  1430. if overwrite[action_type] or stale or not exists:
  1431. do_action = True
  1432. else:
  1433. do_action = False
  1435. if action_type == 'sample':
  1436. if do_action:
  1437. action_kwargs.update(dict(
  1438. output_dir=output_dir,
  1439. sample_id=sample_id
  1440. ))
  1441. if action_type in parcellate_kwargs:
  1442. for key in parcellate_kwargs[action_type]:
  1443. action_kwargs[key] = parcellate_kwargs[action_type][key]
  1444. sample(**action_kwargs, indent=indent)
  1445. else:
  1446. stderr('%sSample exists. Skipping. To resample, run with overwrite=True.\n' %
  1447. (' ' * (indent * 2)))
  1448. elif action_type == 'align':
  1449. if do_action:
  1450. action_kwargs.update(dict(
  1451. output_dir=output_dir,
  1452. alignment_id=alignment_id,
  1453. sample_id=sample_id
  1454. ))
  1455. if action_type in parcellate_kwargs:
  1456. for key in parcellate_kwargs[action_type]:
  1457. action_kwargs[key] = parcellate_kwargs[action_type][key]
  1458. align(**action_kwargs, indent=indent)
  1459. else:
  1460. stderr('%sAlignment exists. Skipping. To re-align, run with overwrite=True.\n' %
  1461. (' ' * (indent * 2)))
  1462. elif action_type == 'label':
  1463. if do_action:
  1464. action_kwargs.update(dict(
  1465. output_dir=output_dir,
  1466. labeling_id=labeling_id,
  1467. alignment_id=alignment_id,
  1468. sample_id=sample_id
  1469. ))
  1470. if action_type in parcellate_kwargs:
  1471. for key in parcellate_kwargs[action_type]:
  1472. action_kwargs[key] = parcellate_kwargs[action_type][key]
  1473. label(**action_kwargs, indent=indent)
  1474. else:
  1475. stderr('%sLabeling exists. Skipping. To re-label, run with overwrite=True.\n' %
  1476. (' ' * (indent * 2)))
  1477. elif action_type == 'evaluate':
  1478. if do_action:
  1479. action_kwargs.update(dict(
  1480. output_dir=output_dir,
  1481. evaluation_id=evaluation_id,
  1482. labeling_id=labeling_id
  1483. ))
  1484. if action_type in parcellate_kwargs:
  1485. for key in parcellate_kwargs[action_type]:
  1486. action_kwargs[key] = parcellate_kwargs[action_type][key]
  1487. evaluate(**action_kwargs, indent=indent)
  1488. else:
  1489. stderr('%sEvaluation exists. Skipping. To re-evaluate, run with overwrite=True.\n' %
  1490. (' ' * (indent * 2)))
  1491. elif action_type == 'parcellate':
  1492. # Copy final files to destination
  1493. if do_action:
  1494. results_copied = False
  1495. if aggregation_id is not None:
  1496. parcellation_kwargs_path = get_path(output_dir, 'output', 'aggregate', aggregation_id)
  1497. assert os.path.exists(parcellation_kwargs_path), ('Aggregation output %s not found' %
  1498. parcellation_kwargs_path)
  1499. if grid_optimized:
  1500. shutil.copy(
  1501. parcellation_kwargs_path, join(parcellation_dir, 'parcellate_kwargs_optimized.yml')
  1502. )
  1504. if evaluation_id is not None:
  1505. evaluation_dir = get_path(output_dir, 'subdir', 'evaluate', evaluation_id)
  1506. for filename in os.listdir(evaluation_dir):
  1507. if filename.endswith(suffix) or filename == PATHS['evaluate']['output']:
  1508. shutil.copy(join(evaluation_dir, filename), join(parcellation_dir, filename))
  1509. if filename == PATHS['evaluate']['output']:
  1510. results_copied = True
  1512. if alignment_id is not None:
  1513. alignment_dir = get_path(output_dir, 'subdir', 'align', alignment_id)
  1514. for filename in os.listdir(alignment_dir):
  1515. if filename.endswith(suffix):
  1516. shutil.copy(join(alignment_dir, filename), join(parcellation_dir, filename))
  1518. labeling_dir = get_path(output_dir, 'subdir', 'label', labeling_id)
  1519. for filename in os.listdir(labeling_dir):
  1520. if filename.endswith(suffix) or (not results_copied and filename == PATHS['label']['output']):
  1521. shutil.copy(join(labeling_dir, filename), join(parcellation_dir, filename))
  1522. else:
  1523. stderr('%sParcellation exists. Skipping. To re-parcellate, run with overwrite=True.\n' %
  1524. (' ' * (indent * 2)))
  1525. else:
  1526. raise ValueError('Unrecognized action_type %s' % action_type)
  1528. with open(output_path, 'w') as f:
  1529. f.write(datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
  1530. assert os.path.exists(output_path)
  1532. stderr('%sTotal time elapsed: %ds\n' % (' ' * (indent * 2), time.time() - t0))
  1543. ######################################
  1544. #
  1545. # PRIVATE HELPER METHODS
  1546. #
  1547. ######################################
  1550. def _get_n_voxels(
  1551. atlas,
  1552. ):
  1553. m, M = atlas.min(), atlas.max()
  1554. if m < 0 or M > 1:
  1555. atlas = minmax_normalize_array(np.clip(atlas, 0, np.inf))
  1556. n_voxels = atlas.sum()
  1558. row = dict(
  1559. n_voxels=n_voxels
  1560. )
  1562. return row
  1565. def _get_atlas_score(
  1566. atlas,
  1567. reference_atlas,
  1568. ):
  1569. with np.errstate(divide='ignore', invalid='ignore'):
  1570. r = np.corrcoef(reference_atlas, atlas)[0, 1]
  1572. row = {
  1573. '%sscore' % REFERENCE_ATLAS_PREFIX: r
  1574. }
  1576. return row
  1579. def _get_evaluation_spcorr(
  1580. atlas,
  1581. evaluation_atlases,
  1582. ):
  1583. row = {}
  1584. evaluation_atlas_names = list(evaluation_atlases.keys())
  1585. for evaluation_atlas_name in evaluation_atlas_names:
  1586. evaluation_atlas = evaluation_atlases[evaluation_atlas_name]
  1587. with np.errstate(divide='ignore', invalid='ignore'):
  1588. r = np.corrcoef(atlas, evaluation_atlas)[0, 1]
  1589. row['%s_score' % evaluation_atlas_name] = r
  1591. return row
  1594. def _get_evaluation_contrasts(
  1595. atlas,
  1596. evaluation_atlases
  1597. ):
  1598. m, M = atlas.min(), atlas.max()
  1599. if m < 0 or M > 1:
  1600. atlas = minmax_normalize_array(np.clip(atlas, 0, np.inf))
  1601. row = {}
  1602. evaluation_atlas_names = list(evaluation_atlases.keys())
  1603. for evaluation_atlas_name in evaluation_atlas_names:
  1604. evaluation_atlas = evaluation_atlases[evaluation_atlas_name]
  1605. denom = atlas.sum()
  1606. if denom:
  1607. contrast = (atlas * evaluation_atlas).sum() / denom
  1608. else:
  1609. contrast = 0
  1610. row['%s_contrast' % evaluation_atlas_name] = contrast
  1612. return row
  1615. def _get_evaluation_row(
  1616. atlas,
  1617. atlas_name,
  1618. reference_atlas_name,
  1619. reference_atlases=None,
  1620. evaluation_atlases=None
  1621. ):
  1622. row = {
  1623. 'parcel': atlas_name,
  1624. '%sname' % REFERENCE_ATLAS_PREFIX: reference_atlas_name
  1625. }
  1626. row.update(_get_n_voxels(atlas))
  1627. if reference_atlases is not None:
  1628. row.update(_get_atlas_score(
  1629. atlas,
  1630. reference_atlases[reference_atlas_name]
  1631. ))
  1632. _reference_atlases = {'%s%s' % (REFERENCE_ATLAS_PREFIX, x): reference_atlases[x] for x in reference_atlases}
  1633. row.update(_get_evaluation_spcorr(
  1634. atlas,
  1635. _reference_atlases
  1636. ))
  1637. if evaluation_atlases is not None:
  1638. _evaluation_atlases = {'%s%s' % (EVALUATION_ATLAS_PREFIX, x): evaluation_atlases[x] for x in evaluation_atlases}
  1639. row.update(_get_evaluation_spcorr(
  1640. atlas,
  1641. _evaluation_atlases
  1642. ))
  1643. row.update(_get_evaluation_contrasts(
  1644. atlas,
  1645. _evaluation_atlases
  1646. ))
  1648. return row
  1651. def _pretty_print_evaluation_row(
  1652. row,
  1653. max_evals=None,
  1654. indent=0
  1655. ):
  1656. to_print = []
  1657. scores = set()
  1658. contrasts = set()
  1659. for col in row:
  1660. if col == 'parcel':
  1661. to_print.append(row[col])
  1662. elif col == 'n_voxels':
  1663. to_print.append('n voxels: %d' % row[col])
  1664. elif (col == ('%sscore' % REFERENCE_ATLAS_PREFIX) or
  1665. (col.startswith(EVALUATION_ATLAS_PREFIX) and col.endswith('_score'))):
  1666. _col = '_'.join(col.split('_')[:-1])
  1667. if max_evals is None or len(scores) < max_evals:
  1668. if col != 'ref_score':
  1669. scores.add(col)
  1670. to_print.append('%s score: %0.3f' % (_col, row[col]))
  1671. elif col.endswith('contrast'):
  1672. _col = '_'.join(col.split('_')[:-1])
  1673. if max_evals is None or len(contrasts) < max_evals:
  1674. contrasts.add(col)
  1675. to_print.append('%s contrast: %0.3f' % (_col, row[col]))
  1677. to_print = ('\n%s' % (' ' * ((indent + 1) * 2))).join(to_print)
  1678. to_print = ' ' * (indent * 2) + to_print
  1680. return to_print
  1683. def _get_atlas_score_from_df(df_scores, subnetwork_id=None, exclude=None, eps=1e-3):
  1684. if exclude is None:
  1685. exclude = []
  1686. if isinstance(exclude, str):
  1687. exclude = [exclude]
  1688. exclude = set(exclude)
  1689. reference_atlas_names = df_scores['%sname' % REFERENCE_ATLAS_PREFIX].unique().tolist()
  1690. if exclude:
  1691. reference_atlas_names = [x for x in reference_atlas_names if not x in exclude]
  1692. parcel_names = df_scores.parcel
  1693. scores = df_scores['%sscore' % REFERENCE_ATLAS_PREFIX]
  1694. target_parcels = reference_atlas_names
  1695. if subnetwork_id:
  1696. target_parcels = [x + '_sub%d' % subnetwork_id for x in target_parcels]
  1697. sel = parcel_names.isin(target_parcels)
  1698. parcel_names = parcel_names[sel].tolist()
  1700. assert set(target_parcels) == set(parcel_names), (f'Parcel names and reference atlas names mismatch. '
  1701. f'Parcel names: {parcel_names}. Reference atlas names: {reference_atlas_names}')
  1703. scores = scores[sel]
  1704. score = np.tanh(np.arctanh(scores * (1 - 2 * eps) + eps).mean())
  1706. return score
  1715. ######################################
  1716. #
  1717. # CONSTANTS
  1718. #
  1719. ######################################
  1721. ACTIONS = dict(
  1722. sample=sample,
  1723. align=align,
  1724. label=label,
  1725. evaluate=evaluate,
  1726. aggregate=aggregate,
  1727. parcellate=parcellate,
  1728. )

model.py at commit fbb979b, no license · at the source

Overview

Authors: Cory Shain1, Evelina Fedorenko2
ORCID iDs: Cory Shain, Evelina Fedorenko
  1. Department of Linguistics, Stanford University, Stanford, CA USA
  2. Department of Brain & Cognitive Sciences and McGovern Institute for Brain Research, Massachusetts Institute of Technology, Cambridge, MA USA
Institutions: Stanford University (United States); McGovern Institute for Brain Research (United States); Massachusetts Institute of Technology (United States)
Journal: Nature communications, volume 17, issue 1, article 8180
Dates: received 9 February 2026; accepted 7 July 2026; published online 13 August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-75745-8 · PMID 42595752 · PMCID PMC13473096 · OpenAlex W7163357038
Open access: gold, a free copy (OpenAlex)
Preprint: osf.io/nd5g3
Status: code verified
Categories: human (organism), cognitive (subfield)
Methods: Spectral & time-frequency, Connectivity, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Statistics, fMRI & imaging
Keywords: Language, Network models
MeSH: Brain*, Connectome*, Language*, Nerve Net*, Adult, Brain Mapping, Female, Humans, Magnetic Resonance Imaging, Male, Young Adult (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: U.S. Department of Health &amp; Human Services | NIH | National Institute of Neurological Disorders and Stroke (NS121471); U.S. Department of Health & Human Services | NIH | National Institute of Neurological Disorders and Stroke (NINDS) (NS121471); NINDS NIH HHS (U01 NS121471)
Citations: cited by 6 papers (Europe PMC); 172 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 10 matches between paragraphs and lines of code.

coryshain/parcellate

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: fbb979b21592f7c2d62d243a759af5d006cf5660, 9 February 2026
Languages: Python (13)
Size: 41 files, 13 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (conda.yml, source/requirements.txt), documentation
Not found: license file, CITATION.cff, tests, continuous integration
Tools: NumPy (5 files), Nilearn (3 files), SciPy (3 files), pandas (2 files), Matplotlib (1 file), Pillow (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
14 files
At the source: github.com/coryshain/parcellate

coryshain/langlocFC

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 4f1100e91ebd11cf40fb2f8c880d410d23830bbe, 8 June 2026
Languages: Python (36)
Size: 51 files, 36 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (conda.yml)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: pandas (27 files), NumPy (21 files), Matplotlib (9 files), Nilearn (9 files), seaborn (6 files), Pillow (4 files), SciPy (3 files), h5py (2 files), NiBabel (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
37 files
At the source: github.com/coryshain/langlocFC

Zenodo 20520281

License: CC-BY-4.0
State: the link answers, verified on 27 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 (5 files), Nilearn (3 files), SciPy (3 files), pandas (2 files), Matplotlib (1 file), Pillow (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
13 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-75745-8.

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:

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

Code and data availability statement

The paper has a code and data 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-75745-8.

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

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 2 keywords, 11 MeSH terms, 3 funders, 151 references.

Cite

This paper

Shain, C., & Fedorenko, E. (2026). A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks. Nature communications, 17(1), 8180. https://doi.org/10.1038/s41467-026-75745-8

BibTeX

@article{shain2026language,
author = {Shain, Cory and Fedorenko, Evelina},
title = {{A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks}},
journal = {Nature communications},
year = {2026},
month = aug,
volume = {17},
number = {1},
pages = {8180},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75745-8},
url = {https://doi.org/10.1038/s41467-026-75745-8},
pmid = {42595752},
pmcid = {PMC13473096}
}

RIS

TY - JOUR
AU - Shain, Cory
AU - Fedorenko, Evelina
TI - A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/08/13
VL - 17
IS - 1
SP - 8180
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75745-8
UR - https://doi.org/10.1038/s41467-026-75745-8
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75745-8",
"type": "article-journal",
"title": "A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks",
"container-title": "Nature communications",
"author": [
{
"family": "Shain",
"given": "Cory"
},
{
"family": "Fedorenko",
"given": "Evelina"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8180",
"DOI": "10.1038/s41467-026-75745-8",
"PMID": "42595752",
"PMCID": "PMC13473096",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75745-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
13
]
]
}
}

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.1162/imag.a.1246 [code]
A 3.5-minute-long reading-based fMRI localizer for the language network.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: h5py, NiBabel, seaborn, 4 other tools, 34 references, author Evelina Fedorenko
[2] doi:10.1038/s41467-026-76598-x
Preserved topography, lateralization, selectivity, and functional connectivity of the language network in older brains.
Journal: Nature communications
In common: cognitive, 37 references, author Evelina Fedorenko
[3] doi:10.1038/s41467-026-72916-5 [code]
Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age.
Journal: Nature communications
In common: 24 references, author Evelina Fedorenko
[4] doi:10.1162/imag.a.1283 [code]
The language network responds robustly to sentences across tasks.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, h5py, seaborn, 5 other tools, 22 references
[5] doi:10.1523/jneurosci.0638-25.2026 [code]
The Extended Language Network: Language-Responsive Brain Areas Whose Contributions to Language Remain To Be Discovered.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: cognitive, 22 references
[6] doi:10.1016/j.neuron.2026.04.011 [code]
Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.
Journal: Neuron
In common: Nilearn, NiBabel, scikit-learn, 3 other tools, 18 references
[7] 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: Nilearn, h5py, Pillow, 7 other tools, cognitive, 12 references
[8] 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: Nilearn, h5py, Pillow, 7 other tools, cognitive, 12 references
[9] doi:10.1038/s41586-026-10691-5 [code]
Mapping the neuronal building blocks of human language with language models.
Journal: Nature
In common: scikit-learn, pandas, SciPy, 2 other tools, cognitive, 12 references
[10] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, NiBabel, scikit-learn, 3 other tools, 10 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.