OSCR

Contrastive learning to fine-tune feature extraction models for the visual cortex.

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] § Results ↔ code/results_utils.py, lines 739–824 · score 0.66 · mfs words, VWFA, FBA, FFA, OFA, OWFA
  2. [2] § Results ↔ code/utils.py, lines 314–355 · score 0.65 · mfs words, VWFA, FBA, FFA, OFA, OWFA
  3. [3] § Methodology › Implementation details ↔ code/models.py, lines 28–77 · score 0.59 · MLP projection head, ReLU, batch, layers, voxels
  4. [4] § Methodology › Baseline feature extraction models ↔ code/results_utils.py, lines 22–84 · score 0.55 · ImageNet, channels, resized, AlexNet
  5. [5] § Methodology › External dataset validation ↔ code/fit_encoding_models_nod.py, lines 212–271 · score 0.51 · AlexNet model, pooled model, Encoding Model, CL tuned, V1, NOD

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 · 890 lines · 43 KB · no license · 2 matches

  2. from models import CLR_model, fmri_reg, get_pooled_CL_model
  3. from utils import get_dataloaders
  4. from sklearn.metrics.pairwise import cosine_similarity
  5. import torchextractor as tx
  6. from tqdm import tqdm
  7. from sklearn.metrics import accuracy_score
  8. from sklearn.preprocessing import StandardScaler
  9. from sklearn.linear_model import LogisticRegression
  10. from torchvision.models.feature_extraction import create_feature_extractor
  11. from torchvision.models import AlexNet_Weights
  12. import torchvision
  13. import joblib
  14. from torchvision.models import AlexNet_Weights
  15. from torchvision import transforms
  16. import torch
  17. import numpy as np
  18. import os
  19. os.environ["OMP_NUM_THREADS"] = '1'
  22. # Get results for image classification task
  23. def image_classification_results(project_dir, subj_num, hemisphere, rois, device, tuning_method='cl', dataset_name='caltech256', save=False,
  24. pooled=False, pooled_h_method='const', save_probs=False):
  26. hemisphere_abbr = 'l' if hemisphere == 'left' else 'r'
  28. if pooled and pooled_h_method == 'avg':
  29. save_path = os.path.join(project_dir, "results", hemisphere_abbr + "h_" +
  30. dataset_name + "_" + tuning_method + "_results_pooled_havg.joblib")
  31. if save_probs:
  32. probs_save_folder = os.path.join(
  33. project_dir, "classification_preds", "Pooled")
  34. elif pooled and pooled_h_method == 'const':
  35. save_path = os.path.join(project_dir, "results", hemisphere_abbr + "h_" +
  36. dataset_name + "_" + tuning_method + "_results_pooled_hconst.joblib")
  37. if save_probs:
  38. probs_save_folder = os.path.join(
  39. project_dir, "classification_preds", "Pooled")
  40. else:
  41. save_path = project_dir + "/results/Subj" + str(subj_num) + "/subj" + str(
  42. subj_num) + "_" + hemisphere_abbr + "h_" + dataset_name + "_" + tuning_method + "_results.joblib"
  43. probs_save_folder = os.path.join(
  44. project_dir, "classification_preds", "Subj" + str(subj_num))
  46. # Seed RNG, define image transforms for alexnet
  47. torch.manual_seed(0)
  48. alex_transform = transforms.Compose([
  49. transforms.Resize(256),
  50. transforms.CenterCrop(224),
  51. transforms.ToTensor(), # convert the images to a PyTorch tensor
  52. # normalize the images color channels
  53. transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])
  54. ])
  56. # Load the Data
  57. if (dataset_name == 'caltech256'):
  58. image_dir = os.path.join(project_dir, "caltech256")
  59. dataset = torchvision.datasets.ImageFolder(
  60. root=image_dir, transform=alex_transform)
  61. elif (dataset_name == 'places365'):
  62. image_dir = os.path.join(project_dir, "places365")
  63. dataset = torchvision.datasets.Places365(
  64. root=image_dir, split='val', small=True, transform=alex_transform)
  65. elif (dataset_name == 'sun397'):
  66. image_dir = os.path.join(project_dir, "sun397")
  67. dataset = torchvision.datasets.SUN397(
  68. root=image_dir, transform=alex_transform, download=True)
  69. elif (dataset_name == 'imagenet'):
  70. image_dir = os.path.join(project_dir, "imagenet")
  71. dataset = torchvision.datasets.ImageNet(
  72. root=image_dir, split='val', transform=alex_transform)
  74. total_num_images = len(dataset)
  75. generator = torch.Generator()
  76. generator.manual_seed(0)
  77. shuffled_idxs = torch.randperm(total_num_images, generator=generator)
  79. train_size = int(0.85 * total_num_images)
  80. test_size = total_num_images - train_size
  81. print(train_size, test_size)
  83. train_idxs = shuffled_idxs[:train_size]
  84. test_idxs = shuffled_idxs[train_size:train_size+test_size]
  86. # Create train and test dataloaders
  87. train_dataset = torch.utils.data.Subset(dataset, train_idxs)
  88. test_dataset = torch.utils.data.Subset(dataset, test_idxs)
  89. train_dataloader = torch.utils.data.DataLoader(
  90. train_dataset, batch_size=256, shuffle=False)
  91. test_dataloader = torch.utils.data.DataLoader(
  92. test_dataset, batch_size=256, shuffle=False)
  94. # Get image features for untuned AlexNet
  95. alex = torch.hub.load('pytorch/vision:v0.10.0', 'alexnet',
  96. weights=AlexNet_Weights.IMAGENET1K_V1)
  97. alex.to(device)
  98. alex.eval()
  100. feature_extractor = create_feature_extractor(
  101. alex, return_nodes=['classifier.5']).to(device)
  102. del alex
  104. # Get untuned alexnet features
  105. train_features_untuned = np.zeros((train_size, 4096))
  106. train_labels = np.zeros(train_size)
  107. for batch_index, data in tqdm(enumerate(train_dataloader), total=len(train_dataloader)):
  108. batch_size = data[0].shape[0]
  109. if batch_index == 0:
  110. low_idx = 0
  111. high_idx = batch_size
  112. else:
  113. low_idx = high_idx
  114. high_idx += batch_size
  115. # Extract features
  116. with torch.no_grad():
  117. ft = feature_extractor(data[0].to(device))
  118. # Flatten the features
  119. ft = torch.hstack([torch.flatten(l, start_dim=1)
  120. for l in ft.values()]).cpu().detach().numpy()
  121. train_features_untuned[low_idx:high_idx] = ft
  122. train_labels[low_idx:high_idx] = data[1]
  123. del ft
  125. test_features_untuned = np.zeros((test_size, 4096))
  126. test_labels = np.zeros(test_size)
  127. for batch_index, data in tqdm(enumerate(test_dataloader), total=len(test_dataloader)):
  128. batch_size = data[0].shape[0]
  129. if batch_index == 0:
  130. low_idx = 0
  131. high_idx = batch_size
  132. else:
  133. low_idx = high_idx
  134. high_idx += batch_size
  135. # Extract features
  136. with torch.no_grad():
  137. ft = feature_extractor(data[0].to(device))
  138. # Flatten the features
  139. ft = torch.hstack([torch.flatten(l, start_dim=1)
  140. for l in ft.values()]).cpu().detach().numpy()
  141. test_features_untuned[low_idx:high_idx] = ft
  142. del ft
  143. test_labels[low_idx:high_idx] = data[1]
  145. del feature_extractor
  147. # Save labels
  148. if save_probs:
  149. np.save(os.path.join(project_dir, "classification_preds",
  150. dataset_name + "_test_labels.npy"), test_labels)
  152. scaler = StandardScaler()
  153. fit_scaler = scaler.fit(train_features_untuned)
  154. train_features_untuned = fit_scaler.transform(train_features_untuned)
  155. test_features_untuned = fit_scaler.transform(test_features_untuned)
  157. print("Fitting linear classifier...")
  158. classifier = LogisticRegression(max_iter=5000).fit(
  159. train_features_untuned, train_labels)
  160. preds = classifier.predict(test_features_untuned)
  161. if save_probs:
  162. untuned_pred_probs = classifier.predict_proba(test_features_untuned)
  163. untuned_pred_probs_save_path = os.path.join(
  164. project_dir, "classification_preds", "Untuned", "untuned_" + dataset_name + "_test_probs.npy")
  165. np.save(untuned_pred_probs_save_path, untuned_pred_probs)
  166. acc = accuracy_score(test_labels, preds) * 100
  167. print("Untuned", acc)
  169. if rois[0] == 'all' and hemisphere == 'left':
  170. rois = ["V1v", "V1d", "V2v", "V2d", "V3v", "V3d", "hV4", "EBA", "FBA-1", "FBA-2",
  171. "OFA", "FFA-1", "FFA-2", "OPA", "PPA", "RSC", "OWFA", "VWFA-1", "VWFA-2", "mfs-words"]
  172. elif rois[0] == 'all' and hemisphere == 'right':
  173. rois = ["V1v", "V1d", "V2v", "V2d", "V3v", "V3d", "hV4", "EBA", "FBA-1", "FBA-2",
  174. "mTL-bodies", "OFA", "FFA-1", "FFA-2", "OPA",
  175. "PPA", "RSC", "OWFA", "VWFA-1", "VWFA-2", "mfs-words", "mTL-words"]
  177. # Go through list of rois
  178. for roi in rois:
  180. print(roi)
  181. found_subj = False
  183. # Get number of voxels for this model
  184. _, _, _, _, num_voxels = get_dataloaders(
  185. project_dir, device, subj_num, hemisphere, roi, 1024, shuffle=False)
  187. if num_voxels < 20 and not pooled:
  188. print(roi, "is too small or empty")
  189. elif num_voxels >= 20 and not pooled:
  190. found_subj = True
  192. # If using pooled models, may need dummy fmri input data matching shape of some other subject's acitvations
  193. elif num_voxels >= 20 and pooled:
  194. found_subj = True
  195. elif num_voxels < 20 and pooled:
  196. for subj_idx in range(2, 9):
  197. _, _, _, _, num_voxels = get_dataloaders(
  198. project_dir, device, subj_idx, hemisphere, roi, 1024, shuffle=False)
  199. if num_voxels > 0:
  200. found_subj = True
  201. break
  203. if found_subj:
  204. if tuning_method == 'cl':
  205. if pooled == True:
  206. # Get voxel counts for ROI across all subjects with the ROI from previously created dictionary
  207. voxel_counts_file = os.path.join(
  208. project_dir, hemisphere_abbr + "h_voxel_counts_rois.joblib")
  209. voxel_counts = joblib.load(voxel_counts_file)
  210. num_voxels_subjs = np.array(voxel_counts[roi])
  211. num_voxels_subjs = np.where(num_voxels_subjs == 0, np.nan, num_voxels_subjs)
  212. avg_voxel_dim = int(np.nanmean(num_voxels_subjs))
  213. if pooled_h_method == 'avg':
  214. # Load pooled CL model
  215. model, _ = get_pooled_CL_model(
  216. num_voxels_subjs, device, avg_voxel_dim)
  217. model_path = os.path.join(
  218. project_dir, "cl_models", hemisphere_abbr + "h_" + roi + "_pooled_model_e30_havg.pt")
  219. elif pooled_h_method == 'const':
  220. h_max_dim = 5741
  221. # Load pooled model
  222. model, _ = get_pooled_CL_model(
  223. num_voxels_subjs, device, h_max_dim)
  224. model_path = os.path.join(
  225. project_dir, "cl_models", hemisphere_abbr + "h_" + roi + "_pooled_model_e30_hconst.pt")
  226. else:
  227. model_dir = os.path.join(project_dir, "cl_models", "Subj" + str(subj_num))
  228. model_path = os.path.join(model_dir, "subj" + \
  229. str(subj_num) + "_" + hemisphere_abbr + \
  230. "h_" + roi + "_model_e30.pt")
  231. h_dim = int(num_voxels*0.8)
  232. z_dim = int(num_voxels*0.2)
  233. model = CLR_model(num_voxels, h_dim, z_dim)
  235. elif tuning_method == 'reg':
  236. model_dir = os.path.join(project_dir, "baseline_models", "nn_reg", "Subj" + str(subj_num))
  237. model_path = os.path.join(model_dir, "subj" + \
  238. str(subj_num) + "_" + hemisphere_abbr + \
  239. "h_" + roi + "_reg_model_e75.pt")
  240. model = fmri_reg(num_voxels)
  242. # Some models are saved differently
  243. try:
  244. model.load_state_dict(torch.load(
  245. model_path, map_location=torch.device('cpu'))[0].state_dict())
  246. except:
  247. try:
  248. model.load_state_dict(torch.load(
  249. model_path, map_location=torch.device('cpu')).state_dict())
  250. except:
  251. model.load_state_dict(torch.load(
  252. model_path, map_location=torch.device('cpu')))
  254. model.to(device)
  255. model.eval()
  256. feature_extractor = tx.Extractor(
  257. model, ["alex.classifier.5"]).to(device)
  259. train_features_tuned = np.zeros((train_size, 4096))
  260. train_labels = np.zeros(train_size)
  261. for batch_index, data in tqdm(enumerate(train_dataloader), total=len(train_dataloader)):
  262. batch_size = data[0].shape[0]
  263. if batch_index == 0:
  264. low_idx = 0
  265. high_idx = batch_size
  266. else:
  267. low_idx = high_idx
  268. high_idx += batch_size
  269. # Extract features
  270. with torch.no_grad():
  271. if tuning_method == 'cl':
  272. fmri_dummy = torch.zeros(
  273. (batch_size, num_voxels)).to(device)
  274. if pooled:
  275. _, alex_out_dict = feature_extractor(
  276. fmri_dummy, data[0], subj_num)
  277. else:
  278. _, alex_out_dict = feature_extractor(
  279. fmri_dummy, data[0])
  280. # _, alex_out_dict = feature_extractor(fmri_dummy, data[0].to(device))
  281. elif tuning_method == 'reg':
  282. _, alex_out_dict = feature_extractor(
  283. data[0].to(device))
  284. ft = alex_out_dict['alex.classifier.5'].detach().cpu().numpy()
  285. train_features_tuned[low_idx:high_idx] = ft
  286. train_labels[low_idx:high_idx] = data[1]
  287. del ft
  289. test_features_tuned = np.zeros((test_size, 4096))
  290. test_labels = np.zeros(test_size)
  291. for batch_index, data in tqdm(enumerate(test_dataloader), total=len(test_dataloader)):
  292. batch_size = data[0].shape[0]
  293. if batch_index == 0:
  294. low_idx = 0
  295. high_idx = batch_size
  296. else:
  297. low_idx = high_idx
  298. high_idx += batch_size
  299. # Extract features
  300. with torch.no_grad():
  301. if tuning_method == 'cl':
  302. fmri_dummy = torch.zeros(
  303. (batch_size, num_voxels)).to(device)
  304. if pooled:
  305. _, alex_out_dict = feature_extractor(
  306. fmri_dummy, data[0].to(torch.float), subj_num)
  307. else:
  308. _, alex_out_dict = feature_extractor(
  309. fmri_dummy, data[0])
  310. # _, alex_out_dict = feature_extractor(fmri_dummy, data[0].to(device))
  311. elif tuning_method == 'reg':
  312. _, alex_out_dict = feature_extractor(
  313. data[0].to(device))
  314. ft = alex_out_dict['alex.classifier.5'].detach().cpu().numpy()
  315. test_features_tuned[low_idx:high_idx] = ft
  316. test_labels[low_idx:high_idx] = data[1]
  317. del ft
  319. scaler = StandardScaler()
  320. fit_scaler = scaler.fit(train_features_tuned)
  321. train_features_tuned = fit_scaler.transform(train_features_tuned)
  322. test_features_tuned = fit_scaler.transform(test_features_tuned)
  324. print("Fitting linear classifier...")
  325. classifier = LogisticRegression(max_iter=5000).fit(
  326. train_features_tuned, train_labels)
  327. preds = classifier.predict(test_features_tuned)
  328. if save_probs:
  329. tuned_pred_probs = classifier.predict_proba(
  330. test_features_tuned)
  331. if pooled:
  332. if pooled_h_method == 'avg':
  333. tuned_pred_probs_save_path = os.path.join(
  334. probs_save_folder, hemisphere_abbr + "h_" + roi + "_pooled_havg_" + dataset_name + "_test_probs.npy")
  335. elif pooled_h_method == 'const':
  336. tuned_pred_probs_save_path = os.path.join(
  337. probs_save_folder, hemisphere_abbr + "h_" + roi + "_pooled_hconst_" + dataset_name + "_test_probs.npy")
  338. else:
  339. if tuning_method == 'cl':
  340. tuned_pred_probs_save_path = os.path.join(probs_save_folder, "subj" + str(
  341. subj_num) + "_" + hemisphere_abbr + "h_" + roi + "_" + dataset_name + "_test_probs.npy")
  342. else:
  343. tuned_pred_probs_save_path = os.path.join(probs_save_folder, "subj" + str(
  344. subj_num) + "_" + hemisphere_abbr + "h_" + roi + "_" + dataset_name + "_reg_test_probs.npy")
  345. np.save(tuned_pred_probs_save_path, tuned_pred_probs)
  346. acc = accuracy_score(test_labels, preds) * 100
  347. # results[roi] = acc
  348. print(roi, acc)
  350. if save:
  351. try:
  352. # Load existing results, add result for roi if not already in the existing results
  353. existing_results = joblib.load(save_path)
  354. if roi not in existing_results.keys():
  355. existing_results[roi] = acc
  356. joblib.dump(existing_results, save_path)
  357. except:
  358. results = {}
  359. results[roi] = acc
  360. joblib.dump(results, save_path)
  363. # Compute lower bound on mutual information between CNN features and ROI response using the testing data
  364. def compute_mi_lower_bound(project_dir, device, subj_num, roi, hemisphere, pooled=False, h_method='const'):
  366. # Temperature tau used to tune CL models
  367. tau = 0.3
  369. hemisphere_abbr = 'l' if hemisphere == 'left' else 'r'
  370. print(roi, hemisphere)
  371. # Get test dataloader
  372. _, test_dataloader, _, test_size, num_voxels = get_dataloaders(
  373. project_dir, device, subj_num, hemisphere, roi, 1024, shuffle=False)
  375. if (num_voxels == 0):
  376. print("Empty ROI")
  377. return -1, -1, -1, -1, -1, -1, -1
  378. elif (num_voxels < 20):
  379. print("Too few voxels")
  380. return -1, -1, -1, -1, -1, -1, -1
  382. if pooled:
  383. # Get list of number of voxels for each subj
  384. present_subjs = []
  385. num_voxels_subjs = []
  386. for subj_idx in range(1, 9):
  387. _, _, _, _, num_voxels = get_dataloaders(
  388. project_dir, device, subj_idx, hemisphere, roi, batch_size=1024)
  389. if num_voxels != 0:
  390. present_subjs.append(subj_idx)
  391. num_voxels_subjs.append(num_voxels)
  393. # Get voxel counts for ROI across all subjects with the ROI from previously created dictionary
  394. voxel_counts_file = os.path.join(
  395. project_dir, hemisphere_abbr + "h_voxel_counts_rois.joblib")
  396. voxel_counts = joblib.load(voxel_counts_file)
  397. num_voxels_subjs = np.array(voxel_counts[roi])
  398. num_voxels_subjs = np.where(num_voxels_subjs == 0, np.nan, num_voxels_subjs)
  399. avg_voxel_dim = int(np.nanmean(num_voxels_subjs))
  401. # Get pooled model
  402. print("Getting CL predictions...")
  403. cl_model_dir = os.path.join(project_dir, "cl_models")
  404. if h_method == 'avg':
  405. h_dim = avg_voxel_dim
  406. # Load pooled CL model
  407. cl_model, _ = get_pooled_CL_model(num_voxels_subjs, device, h_dim)
  408. cl_model_path = os.path.join(
  409. cl_model_dir, hemisphere_abbr + "h_" + roi + "_pooled_model_e30_havg.pt")
  410. elif h_method == 'const':
  411. h_dim = 5741
  412. # Load pooled model
  413. cl_model, _ = get_pooled_CL_model(num_voxels_subjs, device, h_dim)
  414. cl_model_path = os.path.join(
  415. cl_model_dir, hemisphere_abbr + "h_" + roi + "_pooled_model_e30_hconst.pt")
  416. z_dim = int(h_dim * 0.25)
  418. # Some models are saved differently
  419. try:
  420. cl_model.load_state_dict(torch.load(
  421. cl_model_path, map_location=torch.device('cpu'))[0].state_dict())
  422. except:
  423. try:
  424. cl_model.load_state_dict(torch.load(
  425. cl_model_path, map_location=torch.device('cpu')).state_dict())
  426. except:
  427. cl_model.load_state_dict(torch.load(
  428. cl_model_path, map_location=torch.device('cpu')))
  429. else:
  430. # Load CL-tuned model
  431. cl_model_dir = os.path.join(project_dir, "Subj" + str(subj_num))
  432. cl_model_path = os.path.join(cl_model_dir, "subj" + str(subj_num) + "_" + hemisphere_abbr +
  433. "h_" + roi + "_model_e30.pt")
  434. h_dim = int(num_voxels*0.8)
  435. z_dim = int(num_voxels*0.2)
  436. cl_model = CLR_model(num_voxels, h_dim, z_dim)
  437. # Some models are saved differently
  438. try:
  439. cl_model.load_state_dict(torch.load(
  440. cl_model_path, map_location=torch.device('cpu'))[0].state_dict())
  441. except:
  442. try:
  443. cl_model.load_state_dict(torch.load(
  444. cl_model_path, map_location=torch.device('cpu')).state_dict())
  445. except:
  446. cl_model.load_state_dict(torch.load(
  447. cl_model_path, map_location=torch.device('cpu')))
  448. cl_model.to(device)
  449. cl_model.eval()
  451. # Use just 1 batch
  452. K = test_size
  453. # Get outputs (after projection head) for CL model and fMRI data
  454. nn_z = np.zeros((K, z_dim))
  455. fmri_z = np.zeros((K, z_dim))
  456. for batch_index, data in tqdm(enumerate(test_dataloader), total=len(test_dataloader)):
  457. batch_size = data[0].shape[0]
  458. if batch_index == 0:
  459. low_idx = 0
  460. high_idx = batch_size
  461. else:
  462. low_idx = high_idx
  463. high_idx += batch_size
  464. # Extract features
  465. with torch.no_grad():
  466. if pooled:
  467. # Forward function for pooled CL models expects subj num = (1-based) position in list of subjects with roi present (exclusing subjects without roi),
  468. # so need to correct indexing
  469. subj_num_adjusted = present_subjs.index(subj_num) + 1
  470. # print(subj_num_adjusted)
  471. ft_fmri, ft_nn = cl_model(data[0], data[1], subj_num_adjusted)
  472. else:
  473. ft_fmri, ft_nn = cl_model(data[0], data[1])
  474. # Flatten the features, collect them
  475. ft_fmri = ft_fmri.cpu().detach().numpy()
  476. ft_nn = ft_nn.cpu().detach().numpy()
  477. fmri_z[low_idx:high_idx] = ft_fmri
  478. nn_z[low_idx:high_idx] = ft_nn
  480. critic_out = (1 / tau) * cosine_similarity(nn_z, fmri_z)
  481. exp_critic = np.exp(critic_out)
  482. lower_bound_mi = 0
  483. for i in range(K):
  484. lower_bound_mi += np.log(exp_critic[i, i] /
  485. ((1 / K) * np.sum(exp_critic[i, :])))
  486. lower_bound_mi_unscaled = (1 / K) * lower_bound_mi
  488. # Do post-hoc scaling for temp
  489. betas = np.logspace(-2, 2, 1000)
  490. best_beta = 1
  491. best_lower_bound = lower_bound_mi_unscaled
  492. for beta in betas:
  493. tau_new = tau * beta
  494. critic_out = (1 / tau_new) * cosine_similarity(nn_z, fmri_z)
  495. exp_critic = np.exp(critic_out)
  496. lower_bound_mi = 0
  497. for i in range(K):
  498. lower_bound_mi += np.log(exp_critic[i, i] /
  499. ((1 / K) * np.sum(exp_critic[i, :])))
  500. lower_bound_mi = (1 / K) * lower_bound_mi
  501. if lower_bound_mi > best_lower_bound:
  502. best_lower_bound = lower_bound_mi
  503. best_beta = beta
  504. print(best_beta, best_lower_bound)
  506. # Save results
  507. if pooled:
  508. save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(
  509. subj_num) + "_" + hemisphere_abbr + "h_mi_lower_bound_pooled_h" + h_method + ".joblib")
  510. else:
  511. save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  512. hemisphere_abbr + "h_mi_lower_bound.joblib")
  513. try:
  514. results = joblib.load(save_path)
  515. except:
  516. results = {}
  517. results[roi] = lower_bound_mi_unscaled, np.log(
  518. K), best_beta, best_lower_bound
  519. joblib.dump(results, save_path)
  522. # Load test cv results for single subject untuned for all layers and subjects to create figure 2 (matrix of encoding scores)
  523. def load_test_cv_single_subj_untuned_all_layers(project_dir):
  524. num_rows = 46
  525. results_matrix = np.zeros((num_rows, 8)) # Average across subjects
  526. # Keep track of how many subjs have each ROI
  527. row_counters = np.zeros((num_rows))
  528. for subj_num in range(1, 9):
  529. results_folder_path = os.path.join(
  530. project_dir, "best_alex_out_layers_test_cv", "Subj" + str(subj_num))
  531. results_file = os.path.join(
  532. results_folder_path, "best_alex_layers_mat_untuned.npy")
  533. subj_results = np.load(results_file)
  534. for row_idx in range(num_rows):
  535. if subj_results[row_idx, 0] > 0:
  536. row_counters[row_idx] += 1
  537. results_matrix[row_idx, :] += subj_results[row_idx, :]
  538. for row_idx in range(num_rows):
  539. if results_matrix[row_idx, 0] > 0:
  540. results_matrix[row_idx, :] /= row_counters[row_idx]
  541. return results_matrix
  544. # Load test cv results for single subject (untuned, cl-tuned, reg-tuned, pooled-avg, or pooled-specific)
  545. # split options are all, early, higher
  546. def load_test_cv_single_subj_results(project_dir, subj_num, split='all', return_pooled_results=False, roi=None):
  548. results_folder_path = os.path.join(
  549. project_dir, "best_alex_out_layers_test_cv", "Subj" + str(subj_num))
  550. untuned_results_path = os.path.join(
  551. results_folder_path, "best_alphas_voxel_accs_dict_untuned.joblib")
  552. cl_tuned_results_path = os.path.join(
  553. results_folder_path, "best_alphas_voxel_accs_dict_cl_tuned.joblib")
  554. reg_tuned_results_path = os.path.join(
  555. results_folder_path, "best_alphas_voxel_accs_dict_reg_tuned.joblib")
  556. pooled_avg_results_path = os.path.join(
  557. results_folder_path, "best_alphas_voxel_accs_dict_pooled_avg.joblib")
  558. pooled_const_results_path = os.path.join(
  559. results_folder_path, "best_alphas_voxel_accs_dict_pooled_const.joblib")
  561. untuned_results = joblib.load(untuned_results_path)
  562. cl_tuned_results = joblib.load(cl_tuned_results_path)
  563. reg_tuned_results = joblib.load(reg_tuned_results_path)
  564. pooled_avg_results = joblib.load(pooled_avg_results_path)
  565. pooled_const_results = joblib.load(pooled_const_results_path)
  567. if roi is not None:
  568. try:
  569. return untuned_results[roi], reg_tuned_results[roi], cl_tuned_results[roi], pooled_avg_results[roi], pooled_const_results[roi]
  570. except:
  571. return 0, 0, 0, 0, 0
  572. else:
  573. if split == 'all':
  574. rois = ["V1v", "V1d", "V2v", "V2d", "V3v", "V3d", "hV4", "EBA", "FBA-1", "FBA-2",
  575. "mTL-bodies", "OFA", "FFA-1", "FFA-2", "mTL-faces", "aTL-faces", "OPA",
  576. "PPA", "RSC", "OWFA", "VWFA-1", "VWFA-2", "mfs-words", "mTL-words"]
  577. elif split == 'early':
  578. rois = ["V1v", "V1d", "V2v", "V2d", "V3v", "V3d", "hV4"]
  579. elif split == 'higher':
  580. rois = ["EBA", "FBA-1", "FBA-2", "mTL-bodies", "OFA", "FFA-1", "FFA-2", "mTL-faces", "aTL-faces", "OPA",
  581. "PPA", "RSC", "OWFA", "VWFA-1", "VWFA-2", "mfs-words", "mTL-words"]
  582. else:
  583. print("Unsupported Split!")
  584. return
  586. untuned_mean_acc = 0
  587. reg_tuned_mean_acc = 0
  588. cl_tuned_mean_acc = 0
  589. pooled_avg_mean_acc = 0
  590. pooled_const_mean_acc = 0
  591. cl_over_untuned_voxel_mean_percentage = 0
  592. cl_over_reg_tuned_voxel_mean_percentage = 0
  593. pooled_avg_over_cl_voxel_mean_percentage = 0
  594. pooled_const_over_cl_voxel_mean_percentage = 0
  595. num_rois_present = 0
  597. for hemi in ["lh", "rh"]:
  598. for roi in rois:
  599. try:
  600. untuned_hemi_roi_results = untuned_results[hemi + '_' + roi]
  601. reg_tuned_hemi_roi_results = reg_tuned_results[hemi + '_' + roi]
  602. cl_tuned_hemi_roi_results = cl_tuned_results[hemi + '_' + roi]
  603. pooled_avg_hemi_roi_results = pooled_avg_results[hemi + '_' + roi]
  604. pooled_const_hemi_roi_results = pooled_const_results[hemi + '_' + roi]
  606. untuned_mean_acc += untuned_hemi_roi_results.mean()
  607. reg_tuned_mean_acc += reg_tuned_hemi_roi_results.mean()
  608. cl_tuned_mean_acc += cl_tuned_hemi_roi_results.mean()
  609. pooled_avg_mean_acc += pooled_avg_hemi_roi_results.mean()
  610. pooled_const_mean_acc += pooled_const_hemi_roi_results.mean()
  612. cl_over_untuned_voxel_mean_percentage += np.count_nonzero(
  613. cl_tuned_hemi_roi_results - untuned_hemi_roi_results > 0) / cl_tuned_hemi_roi_results.shape[0]
  614. cl_over_reg_tuned_voxel_mean_percentage += np.count_nonzero(
  615. cl_tuned_hemi_roi_results - reg_tuned_hemi_roi_results > 0) / cl_tuned_hemi_roi_results.shape[0]
  617. pooled_avg_over_cl_voxel_mean_percentage += np.count_nonzero(
  618. pooled_avg_hemi_roi_results - cl_tuned_hemi_roi_results > 0) / cl_tuned_hemi_roi_results.shape[0]
  619. pooled_const_over_cl_voxel_mean_percentage += np.count_nonzero(
  620. pooled_const_hemi_roi_results - cl_tuned_hemi_roi_results > 0) / cl_tuned_hemi_roi_results.shape[0]
  622. num_rois_present += 1
  623. except:
  624. pass # Skip missing ROIs
  626. untuned_mean_acc /= num_rois_present
  627. reg_tuned_mean_acc /= num_rois_present
  628. cl_tuned_mean_acc /= num_rois_present
  629. pooled_avg_mean_acc /= num_rois_present
  630. pooled_const_mean_acc /= num_rois_present
  632. cl_over_untuned_voxel_mean_percentage /= num_rois_present
  633. cl_over_reg_tuned_voxel_mean_percentage /= num_rois_present
  634. pooled_avg_over_cl_voxel_mean_percentage /= num_rois_present
  635. pooled_const_over_cl_voxel_mean_percentage /= num_rois_present
  637. if return_pooled_results:
  638. return pooled_avg_mean_acc, pooled_const_mean_acc, pooled_avg_over_cl_voxel_mean_percentage, pooled_const_over_cl_voxel_mean_percentage
  639. else:
  640. return untuned_mean_acc, reg_tuned_mean_acc, cl_tuned_mean_acc, cl_over_untuned_voxel_mean_percentage, cl_over_reg_tuned_voxel_mean_percentage
  643. # Return matrix of results for each layer/roi
  644. def load_test_cv_single_subj_results_all_layers(project_dir, subj_num):
  646. results_folder_path = os.path.join(
  647. project_dir, "best_alex_out_layers_test_cv", "Subj" + str(subj_num))
  648. untuned_results_path = os.path.join(
  649. results_folder_path, "best_alex_layers_mat_untuned.npy")
  650. cl_tuned_results_path = os.path.join(
  651. results_folder_path, "best_alex_layers_mat_cl_tuned.npy")
  652. reg_tuned_results_path = os.path.join(
  653. results_folder_path, "best_alex_layers_mat_reg_tuned.npy")
  654. pooled_avg_results_path = os.path.join(
  655. results_folder_path, "best_alex_layers_mat_pooled_avg.npy")
  656. pooled_const_results_path = os.path.join(
  657. results_folder_path, "best_alex_layers_mat_pooled_const.npy")
  659. untuned_results = np.load(untuned_results_path)
  660. cl_tuned_results = np.load(cl_tuned_results_path)
  661. reg_tuned_results = np.load(reg_tuned_results_path)
  662. pooled_avg_results = np.load(pooled_avg_results_path)
  663. pooled_const_results = np.load(pooled_const_results_path)
  665. return untuned_results, cl_tuned_results, reg_tuned_results, pooled_avg_results, pooled_const_results
  669. # Generate embeddings for test images from untuned, CL, or regression-tuned models. Save corresponding NSD IDs of images.
  670. # Options for tuning_method are 'untuned', 'CL' or 'reg'
  671. # Optionally get embeddings from best layer for encoding instead of final layer
  672. def save_embeddings(project_dir, subj_num, hemisphere, roi, device, tuning_method='CL', use_best_intermediate_layer=False, use_other_subj_images=False, cross_subj_num=1):
  674. # Strip whitespace from roi to handle cases where it comes from files with trailing spaces
  675. roi = roi.strip()
  677. hemisphere_abbr = 'l' if hemisphere == 'left' else 'r'
  679. if tuning_method == 'CL':
  680. if use_best_intermediate_layer:
  681. if use_other_subj_images:
  682. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  683. hemisphere_abbr + "h_" + roi + "_cl_embeddings_best_encoding_layer_cross_subj" + str(cross_subj_num) + ".npy")
  684. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  685. hemisphere_abbr + "h_" + roi + "_cl_embeddings_best_encoding_layer_cross_subj" + str(cross_subj_num) + "_img_ids.npy")
  686. else:
  687. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  688. hemisphere_abbr + "h_" + roi + "_cl_embeddings_best_encoding_layer.npy")
  689. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  690. hemisphere_abbr + "h_" + roi + "_cl_embeddings_best_encoding_layer_img_ids.npy")
  691. else:
  692. if use_other_subj_images:
  693. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  694. hemisphere_abbr + "h_" + roi + "_cl_embeddings_cross_subj" + str(cross_subj_num) + ".npy")
  695. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  696. hemisphere_abbr + "h_" + roi + "_cl_embeddings_cross_subj" + str(cross_subj_num) + "_img_ids.npy")
  697. else:
  698. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  699. hemisphere_abbr + "h_" + roi + "_cl_embeddings.npy")
  700. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  701. hemisphere_abbr + "h_" + roi + "_cl_embeddings_img_ids.npy")
  702. elif tuning_method == 'reg':
  703. if use_best_intermediate_layer:
  704. if use_other_subj_images:
  705. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  706. hemisphere_abbr + "h_" + roi + "_reg_embeddings_best_encoding_layer_cross_subj" + str(cross_subj_num) + ".npy")
  707. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  708. hemisphere_abbr + "h_" + roi + "_reg_embeddings_best_encoding_layer_cross_subj" + str(cross_subj_num) + "_img_ids.npy")
  709. else:
  710. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  711. hemisphere_abbr + "h_" + roi + "_reg_embeddings_best_encoding_layer.npy")
  712. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  713. hemisphere_abbr + "h_" + roi + "_reg_embeddings_best_encoding_layer_img_ids.npy")
  714. else:
  715. if use_other_subj_images:
  716. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  717. hemisphere_abbr + "h_" + roi + "_reg_embeddings_cross_subj" + str(cross_subj_num) + ".npy")
  718. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  719. hemisphere_abbr + "h_" + roi + "_reg_embeddings_cross_subj" + str(cross_subj_num) + "_img_ids.npy")
  720. else:
  721. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  722. hemisphere_abbr + "h_" + roi + "_reg_embeddings.npy")
  723. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  724. hemisphere_abbr + "h_" + roi + "_reg_embeddings_img_ids.npy")
  725. elif tuning_method == 'untuned':
  726. features_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_untuned_embeddings.npy")
  727. ids_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_untuned_embeddings_img_ids.npy")
  729. if use_other_subj_images:
  730. _, _, _, _, num_voxels = get_dataloaders(project_dir,
  731. device, subj_num, hemisphere, roi, batch_size=1024, use_all_data=False, shuffle=False, return_nsd_id=True)
  732. _, test_dataloader, _, test_size, _ = get_dataloaders(project_dir,
  733. device, cross_subj_num, hemisphere, roi, batch_size=1024, use_all_data=False, shuffle=False, return_nsd_id=True)
  735. else:
  736. _, test_dataloader, _, test_size, num_voxels = get_dataloaders(project_dir,
  737. device, subj_num, hemisphere, roi, batch_size=1024, use_all_data=False, shuffle=False, return_nsd_id=True)
  739. # Load tuned model
  740. if tuning_method == 'CL':
  741. model_dir = os.path.join(project_dir, "cl_models", "Subj" + str(subj_num))
  742. model_path = os.path.join(model_dir, "subj" + str(subj_num) + "_" + hemisphere_abbr + "h_" + roi + "_model_e30.pt")
  743. h_dim = int(num_voxels*0.8)
  744. z_dim = int(num_voxels*0.2)
  745. model = CLR_model(num_voxels, h_dim, z_dim)
  746. elif tuning_method == 'reg':
  747. model_dir = os.path.join(project_dir, "baseline_models", "nn_reg", "Subj" + str(subj_num))
  748. model_path = os.path.join(model_dir, "subj" + \
  749. str(subj_num) + "_" + hemisphere_abbr + \
  750. "h_" + roi + "_reg_model_e75.pt")
  751. model = fmri_reg(num_voxels)
  752. elif tuning_method == 'untuned':
  753. model = torch.hub.load('pytorch/vision:v0.10.0', 'alexnet',
  754. weights=AlexNet_Weights.IMAGENET1K_V1)
  756. # Some tuned models are saved differently
  757. if tuning_method == 'CL' or tuning_method == 'reg':
  758. try:
  759. model.load_state_dict(torch.load(
  760. model_path, map_location=torch.device('cpu'), weights_only=False)[0].state_dict())
  761. except:
  762. try:
  763. model.load_state_dict(torch.load(
  764. model_path, map_location=torch.device('cpu'), weights_only=False).state_dict())
  765. except:
  766. model.load_state_dict(torch.load(
  767. model_path, map_location=torch.device('cpu'), weights_only=False))
  768. model.to(device)
  769. model.eval()
  771. roi_names = ["V1v", "V1d", "V2v", "V2d", "V3v", "V3d", "hV4", "EBA", "FBA-1", "FBA-2",
  772. "mTL-bodies", "OFA", "FFA-1", "FFA-2", "mTL-faces", "OPA",
  773. "PPA", "RSC", "OWFA", "VWFA-1", "VWFA-2", "mfs-words", "mTL-words"]
  774. hemis = ["lh", "rh"]
  775. layers = [
  776. "features.2", "features.5", "features.7", "features.9",
  777. "features.12", "classifier.2", "classifier.5", "classifier.6"
  778. ]
  779. layer_dim_dict = {"features.2": 46656, "features.5": 32448, "features.7": 64896, "features.9": 43264,
  780. "features.12": 9216, "classifier.2": 4096, "classifier.5": 4096, "classifier.6": 1000}
  781. if tuning_method == 'CL' or tuning_method == 'reg':
  782. if use_best_intermediate_layer:
  783. _, results, _, _, _ = load_test_cv_single_subj_results_all_layers(project_dir, subj_num)
  784. results_dict = {}
  785. counter = 0
  786. for hemi in hemis:
  787. for roi_name in roi_names:
  788. current_roi = hemi + "_" + roi_name
  789. if not np.isnan(results[counter]).any():
  790. results_dict[current_roi] = results[counter]
  791. counter += 1
  792. best_layer = layers[np.argmax(results_dict[hemisphere_abbr + "h_" + roi])]
  793. print("Using layer:", best_layer)
  794. feature_extractor = tx.Extractor(
  795. model, ["alex." + best_layer]).to(device)
  796. output_dim = layer_dim_dict[best_layer]
  797. else:
  798. feature_extractor = tx.Extractor(
  799. model, ["alex.classifier.5"]).to(device)
  800. output_dim = 4096
  801. elif tuning_method == 'reg':
  802. if use_best_intermediate_layer:
  803. _, _, results, _, _ = load_test_cv_single_subj_results_all_layers(project_dir, subj_num)
  804. results_dict = {}
  805. counter = 0
  806. for hemi in hemis:
  807. for roi_name in roi_names:
  808. current_roi = hemi + "_" + roi_name
  809. if not np.isnan(results[counter]).any():
  810. results_dict[current_roi] = results[counter]
  811. counter += 1
  812. best_layer = layers[np.argmax(results_dict[hemisphere_abbr + "h_" + roi])]
  813. print("Using layer:", best_layer)
  814. feature_extractor = tx.Extractor(
  815. model, ["alex." + best_layer]).to(device)
  816. output_dim = layer_dim_dict[best_layer]
  817. else:
  818. feature_extractor = tx.Extractor(
  819. model, ["alex.classifier.5"]).to(device)
  820. output_dim = 4096
  821. elif tuning_method == 'untuned':
  822. feature_extractor = tx.Extractor(
  823. model, ["classifier.5"]).to(device)
  824. output_dim = 4096
  826. features = np.zeros((test_size, output_dim))
  827. ids = np.zeros(test_size)
  828. for batch_index, data in tqdm(enumerate(test_dataloader), total=len(test_dataloader)):
  829. batch_size = data[0].shape[0]
  830. if batch_index == 0:
  831. low_idx = 0
  832. high_idx = batch_size
  833. else:
  834. low_idx = high_idx
  835. high_idx += batch_size
  836. # Extract features
  837. with torch.no_grad():
  838. if tuning_method == 'CL':
  839. fmri_dummy = torch.zeros(
  840. (batch_size, num_voxels)).to(device)
  841. _, alex_out_dict = feature_extractor(fmri_dummy, data[1])
  842. # _, alex_out_dict = feature_extractor(fmri_dummy, data[0].to(device))
  843. elif tuning_method == 'reg':
  844. _, alex_out_dict = feature_extractor(data[1].to(device))
  845. elif tuning_method == 'untuned':
  846. _, alex_out_dict = feature_extractor(data[1].to(device))
  847. if tuning_method == 'CL' or tuning_method == 'reg':
  848. if use_best_intermediate_layer:
  849. ft = alex_out_dict["alex." + best_layer].detach().cpu().numpy().reshape(high_idx-low_idx, -1)
  850. else:
  851. ft = alex_out_dict['alex.classifier.5'].detach().cpu().numpy()
  852. elif tuning_method == 'untuned':
  853. ft = alex_out_dict['classifier.5'].detach().cpu().numpy()
  854. features[low_idx:high_idx] = ft
  855. ids[low_idx:high_idx] = data[3]
  856. del ft
  858. # Save features and ids
  859. np.save(features_save_path, features)
  860. np.save(ids_save_path, ids)
  863. # Save fmri responses for test images. Save corresponding NSD IDs of images.
  864. # Use to match format of save_embeddings function.
  865. def save_test_fmri_responses(project_dir, subj_num, hemisphere, roi, device):
  867. # Strip whitespace from roi to handle cases where it comes from files with trailing spaces
  868. roi = roi.strip()
  870. hemisphere_abbr = 'l' if hemisphere == 'left' else 'r'
  872. fmri_responses_save_path = os.path.join(project_dir, "results", "Subj" + str(subj_num), "subj" + str(subj_num) + "_" +
  873. hemisphere_abbr + "h_" + roi + "_test_fmri_responses.npy")
  875. _, test_dataloader, _, test_size, num_voxels = get_dataloaders(project_dir,
  876. device, subj_num, hemisphere, roi, batch_size=1024, use_all_data=False, shuffle=False, return_nsd_id=True)
  878. fmri_responses = np.zeros((test_size, num_voxels))
  879. for batch_index, data in tqdm(enumerate(test_dataloader), total=len(test_dataloader)):
  880. batch_size = data[0].shape[0]
  881. if batch_index == 0:
  882. low_idx = 0
  883. high_idx = batch_size
  884. else:
  885. low_idx = high_idx
  886. high_idx += batch_size
  887. fmri_responses[low_idx:high_idx] = data[0]
  888. del data
  890. np.save(fmri_responses_save_path, fmri_responses)

results_utils.py at commit af1f772, no license · at the source

Overview

Authors: Alex Mulrooney1, Zhi Li1, Austin J Brockmeier1,2
ORCID iDs: Alex Mulrooney
  1. Department of Electrical and Computer Engineering, University of Delaware, Newark, Delaware, United States of America
  2. Department of Computer and Information Sciences, University of Delaware, Newark, Delaware, United States of America
Institutions: University of Delaware (United States)
Journal: PLoS computational biology, volume 22, issue 8, article e1014656
Dates: received 28 August 2025; accepted 3 August 2026; published online 17 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014656 · PMID 42607104 · PMCID PMC13492995 · OpenAlex W7203628548
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), human (organism), systems (subfield)
Methods: Smoothing, state filtering, decompositions, Machine learning, Statistics, Connectivity, fMRI & imaging
MeSH: Models, Neurological*, Visual Cortex*, Algorithms, Brain Mapping, Computational Biology, Convolutional Neural Networks, Humans, Image Processing, Computer-Assisted, Magnetic Resonance Imaging (* major topic)
Topic: Face Recognition and Perception (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: University of Delaware Department of Electrical and Computer Engineering First Year Fellowship; University of Delaware College of Engineering (None); University of Delaware Research Foundation; University of Delaware Electrical and Computer Engineering Department (None); Office of Naval Research (N00014-24-1-2259); University of Delaware General University Research fund
Citations: cited by 1 paper (Europe PMC); 59 references in the paper

Abstract

Predicting the neural response to natural images in the visual cortex requires extracting relevant features from the images and relating those feature to the observed responses. In this work, we optimize the feature extraction in order to maximize the information shared between the image features and the neural response across voxels in a given region of interest (ROI) extracted from the BOLD signal measured by functional magnetic resonance imaging (fMRI). We adapt contrastive learning (CL) to fine-tune a convolutional neural network, which was pretrained for image classification, such that a mapping of a given image’s features are more similar to the corresponding fMRI response than to the responses to other images. We exploit the Natural Scenes Dataset as organized for the Algonauts Project, which contains the high-resolution fMRI responses of eight subjects to tens of thousands of naturalistic images. We show that CL fine-tuning creates feature extraction models that enable higher encoding accuracy in both early and higher visual ROIs as compared to the features from the pretrained network. Quantitatively, the performance is similar to a baseline approach that directly uses a regression loss at the output of the network to tune it for fMRI response encoding. We investigate inter-subject transfer of the CL fine-tuned models, including subjects from the Natural Object Dataset, another lower-resolution dataset with 9 subjects. We also pool subjects for fine-tuning, which further improves encoding performance in early ROIs. Finally, we examine the performance of the fine-tuned models on common image classification tasks, explore the landscape of ROI-specific models by applying dimensionality reduction on the Bhattacharya dissimilarity matrix created using the predictions on those tasks, and show that these landscapes match those based on representational similarity analysis. Finally, we generate images via Stable Diffusion based on vector-space prompts created by aligning the CL-tuned models embeddings for different ROIs, showing that generated images have similar embeddings to the original but that estimates of the intrinsic dimension are lower for generated versus original representations.

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.

alexmul1114/fmri_encoding_contrastive_learning

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: af1f772661f2c5789154d2f1078b077569fceba4, 18 August 2026
Languages: Python (13), Jupyter (1), Shell (1)
Size: 564 files, 15 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, environment (requirements.txt), 1 notebook
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: PyTorch (13 files), NumPy (8 files), scikit-learn (4 files), Pillow (3 files), SciPy (3 files), DataLad (1 file), Matplotlib (1 file), NiBabel (1 file), pandas (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
16 files
At the source: github.com/alexmul1114/fmri_encoding_contrastive_learning

The paper's code and data availability statement is in the Data section.

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;
  • 15 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

Datasets cited

Data Availability

All data we use in our paper is publicly available. The Natural Scenes Dataset is available at https://naturalscenesdataset.org/. The Natural Object Dataset is available at https://openneuro.org/datasets/ds004496/versions/2.1.2. The ImageNet dataset is available at https://www.image-net.org/. The Caltech256 dataset is available at https://data.caltech.edu/records/nyy15-4j048. The Places365 dataset is available at https://data.caltech.edu/records/nyy15-4j048. The MS-COCO dataset is available at https://cocodataset.org/#home. We provide code for reproducing our results as well as instructions for obtaining the data on Github at https://github.com/alexmul1114/fmri_encoding_contrastive_learning.

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

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 9 MeSH terms, 6 funders, 36 references.

Cite

This paper

Mulrooney, A., Li, Z., & Brockmeier, A. J. (2026). Contrastive learning to fine-tune feature extraction models for the visual cortex. PLoS computational biology, 22(8), e1014656. https://doi.org/10.1371/journal.pcbi.1014656

BibTeX

@article{mulrooney2026contrastive,
author = {Mulrooney, Alex and Li, Zhi and Brockmeier, Austin J},
title = {{Contrastive learning to fine-tune feature extraction models for the visual cortex}},
journal = {PLoS computational biology},
year = {2026},
month = aug,
volume = {22},
number = {8},
pages = {e1014656},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014656},
url = {https://doi.org/10.1371/journal.pcbi.1014656},
pmid = {42607104},
pmcid = {PMC13492995}
}

RIS

TY - JOUR
AU - Mulrooney, Alex
AU - Li, Zhi
AU - Brockmeier, Austin J
TI - Contrastive learning to fine-tune feature extraction models for the visual cortex
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/08/17
VL - 22
IS - 8
SP - e1014656
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014656
UR - https://doi.org/10.1371/journal.pcbi.1014656
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014656",
"type": "article-journal",
"title": "Contrastive learning to fine-tune feature extraction models for the visual cortex",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Mulrooney",
"given": "Alex"
},
{
"family": "Li",
"given": "Zhi"
},
{
"family": "Brockmeier",
"given": "Austin J"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "8",
"page": "e1014656",
"DOI": "10.1371/journal.pcbi.1014656",
"PMID": "42607104",
"PMCID": "PMC13492995",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014656",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
17
]
]
}
}

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/s41597-026-07248-6 [code]
A large-scale fMRI dataset for vision-language semantic association.
Journal: Scientific data
In common: Pillow, NiBabel, PyTorch, 5 other tools, fMRI, 8 references
[2] doi:10.1038/s42003-026-10169-0 [code]
Shared representations in brains and models reveal a two-route cortical organization during scene perception.
Journal: Communications biology
In common: Pillow, NiBabel, PyTorch, 5 other tools, 7 references
[3] doi:10.7554/elife.107933 [code]
Modality-agnostic decoding of vision and language from fMRI.
Journal: eLife
In common: Pillow, NiBabel, PyTorch, 5 other tools, fMRI, 5 references
[4] doi:10.1371/journal.pcbi.1014263 [code]
MIRAGE: Robust multi-modal architectures translate fMRI-to-image models from vision to mental imagery.
Journal: PLoS computational biology
In common: Pillow, NiBabel, PyTorch, 5 other tools, fMRI, 4 references
[5] doi:10.1038/s41597-025-05174-7 [code]
A large-scale MEG and EEG dataset for object recognition in naturalistic scenes
Journal: n/a
In common: Pillow, NiBabel, PyTorch, 5 other tools, fMRI, 4 references
[6] doi:10.1162/imag.a.1309 [code]
Probing the content of semantic representations in body-selective regions.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Pillow, NiBabel, PyTorch, 5 other tools, 4 references
[7] doi:10.1162/imag.a.1256 [code]
Gamer in the scanner: Event-related analysis of fMRI activity during retro videogame play guided by automated annotations of game content.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: DataLad, Pillow, NiBabel, 6 other tools, fMRI, 1 reference
[8] doi:10.1162/imag.a.1142 [code]
Encoding neural representations of time-continuous stimulus-response transformations in the human brain with advanced deep neural networks.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Pillow, PyTorch, NumPy, fMRI, systems, 5 references
[9] doi:10.1038/s41467-026-76098-y [code]
A single computational objective can produce specialization of streams in visual cortex.
Journal: Nature communications
In common: Pillow, NiBabel, PyTorch, 5 other tools, 3 references
[10] doi:10.1523/jneurosci.0038-26.2026 [code]
Multidimensional Feature Tuning in Category Selective Areas of Human Visual Cortex.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: Pillow, NiBabel, PyTorch, 5 other tools, fMRI, systems, 2 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.