From 003cd8fc7ca49b37651c6649490fb6f210f39fb1 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Mon, 20 Feb 2017 16:58:03 -0500 Subject: [PATCH 01/11] Adding isfc functions --- brainiak/__init__.py | 1 + brainiak/isfc.py | 150 +++++++++++++++++++++++++++++++++ examples/isfc/isfc.py | 62 ++++++++++++++ examples/isfc/requirements.txt | 2 + tests/test_isfc.py | 38 +++++++++ 5 files changed, 253 insertions(+) create mode 100644 brainiak/isfc.py create mode 100644 examples/isfc/isfc.py create mode 100644 examples/isfc/requirements.txt create mode 100644 tests/test_isfc.py diff --git a/brainiak/__init__.py b/brainiak/__init__.py index 10e76d6a9..e25d80792 100644 --- a/brainiak/__init__.py +++ b/brainiak/__init__.py @@ -14,6 +14,7 @@ """Brain Imaging Analysis Kit.""" import sys +from . import isfc if sys.version_info < (3, 4): raise Exception( diff --git a/brainiak/isfc.py b/brainiak/isfc.py new file mode 100644 index 000000000..33e5bf29b --- /dev/null +++ b/brainiak/isfc.py @@ -0,0 +1,150 @@ +# Copyright 2017 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Intersubject analyses (ISC/ISFC) + +Functions for computing intersubject correlation (ISC) and intersubject +functional correlation (ISFC), and utility functions for loading subject data +files and masks +""" + +# Authors: Christopher Baldassano and Mor Regev +# Princeton University, 2017 + +from brainiak.fcma.util import compute_correlation +from scipy import stats +import nibabel as nib +import numpy as np +from matplotlib import pyplot as plt +from scipy.cluster.hierarchy import fcluster, linkage + + +def isc(D, collapse_subj=True): + """Intersubject correlation + + For each voxel, computes the correlation of each subject's timecourse with + the mean of all other subjects' timecourses. By default the result is + averaged across subjects, unless collapse_subj is set to False. + + Parameters + ---------- + D : voxel by time by subject ndarray + fMRI data for which to compute ISC + + collapse_subj : bool, default:True + Whether to average across subjects before returning result + + Returns + ------- + ISC : voxel ndarray (or voxel by subject ndarray, if collapse_subj=False) + pearson correlation for each voxel, across subjects + """ + + n_vox = D.shape[0] + n_subj = D.shape[2] + ISC = np.zeros((n_vox, n_subj)) + + for loo_subj in range(n_subj): + group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) + subj = D[:, :, loo_subj] + for v in range(n_vox): + ISC[v, loo_subj] = stats.pearsonr(group[v, :], subj[v, :])[0] + + if collapse_subj: + ISC = np.mean(ISC, axis=1) + return ISC + + +def isfc(D, collapse_subj=True): + """Intersubject functional correlation + + Computes the correlation between the timecoure of each voxel in each subject + with the average of all other subjects' timecourses in *all* voxels. By + default the result is averaged across subjects, unless collapse_subj is set + to False. + + Uses the high performance compute_correlation routine from fcma.util + + Parameters + ---------- + D : voxel by time by subject ndarray + fMRI data for which to compute ISFC + + collapse_subj : bool, default:True + Whether to average across subjects before returning result + + Returns + ------- + ISFC : voxel by voxel ndarray + (or voxel by voxel by subject ndarray, if collapse_subj=False) + pearson correlation between all pairs of voxels, across subjects + """ + + n_vox = D.shape[0] + n_subj = D.shape[2] + ISFC = np.zeros((n_vox, n_vox, n_subj)) + for loo_subj in range(D.shape[2]): + group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) + subj = D[:, :, loo_subj] + ISFC[:, :, loo_subj] = compute_correlation(group, subj) + ISFC[:, :, loo_subj] = (ISFC[:, :, loo_subj] + + ISFC[:, :, loo_subj].T) / 2 + + if collapse_subj: + ISFC = np.mean(ISFC, axis=2) + return ISFC + + +def load_subjects_nii(data_files, mask_file, mask_func=None): + """Loading masked nifti data into matrix + + Given a list of subject data files and a mask, loads voxel timecourses + inside mask into a matrix, which is returned along with the voxel + coordinates of the mask. If mask_func is given, it specifies which mask + values should be included (e.g. for thresholding a continuous-valued mask) + + Parameters + ---------- + data_files : list of filenames of subject nii files + + mask_file : filename of mask nii + + mask_func : Callable[[ndarray], bool] : default x>0 + + Returns + ------- + (D, coords) : tuple of voxel by time by subject ndarray, and + x,y,x (as provided by nibabel) coordinates of voxel locations (as a + tuple of ndarrays + """ + + mask_nii = nib.load(mask_file) + mask = mask_nii.get_data() + if mask_func is not None: + mask = mask_func(mask) + else: + mask = mask > 0 + mask_shape = mask_nii.shape + coords = np.where(mask) + + data_shape = nib.load(data_files[0]).shape + D = np.zeros((np.sum(mask), data_shape[3], len(data_files))) + for s in range(len(data_files)): + nii = nib.load(data_files[s]) + if nii.shape != data_shape: + raise ValueError("Data has different shapes across subjects") + if nii.shape[:3] != mask_shape: + raise ValueError("Data and mask have different shapes") + D[:, :, s] = nii.get_data()[mask] + + return (D, coords) diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py new file mode 100644 index 000000000..ceab42690 --- /dev/null +++ b/examples/isfc/isfc.py @@ -0,0 +1,62 @@ +# Copyright 2017 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Intersubject analyses (ISC/ISFC) + +Functions for computing intersubject correlation (ISC) and intersubject +functional correlation (ISFC), and utility functions for loading subject data +files and masks +""" + +# Authors: Christopher Baldassano and Mor Regev +# Princeton University, 2017 + +import brainiak +import nibabel as nib +import numpy as np +from matplotlib import pyplot as plt +from scipy.cluster.hierarchy import fcluster, linkage +import sys, os +dir = os.path.dirname(__file__) + +brain_fname = os.path.join(dir,'avg152T1_gray_3mm.nii.gz') +fnames = [os.path.join(dir,'sub-0' + format(subj, '02') + '-task-intact1.nii.gz') for + subj in np.arange(1, 5)] + +print('Loading data from ', len(fnames), ' subjects...') +D, coords = brainiak.isfc.load_subjects_nii(fnames, brain_fname, lambda x: x > 50) + +print('Calculating ISC on ', D.shape[0], ' voxels') +ISC = brainiak.isfc.isc(D) +ISC[np.isnan(ISC)] = 0 + +brain_nii = nib.load(brain_fname) +ISC_vol = np.zeros(brain_nii.shape) +ISC_vol[coords] = ISC +ISC_nifti = nib.Nifti1Image(ISC_vol, + brain_nii.get_affine(), brain_nii.header) +nib.save(ISC_nifti, 'ISC.nii.gz') + +ISC_mask = ISC > 0.2 +print('Calculating ISFC on ', np.sum(ISC_mask), ' voxels...') +D_masked = D[ISC_mask, :, :] +ISFC = brainiak.isfc.isfc(D_masked) +coords_ISFC = [x[ISC_mask] for x in coords] + +print('Clustering ISFC...') +Z = linkage(ISFC, 'ward') +z = fcluster(Z, 2, criterion='maxclust') +clust_inds = np.argsort(z) + +plt.imshow(ISFC[np.ix_(clust_inds,clust_inds)]) +plt.show() diff --git a/examples/isfc/requirements.txt b/examples/isfc/requirements.txt new file mode 100644 index 000000000..f0b259309 --- /dev/null +++ b/examples/isfc/requirements.txt @@ -0,0 +1,2 @@ +matplotlib + diff --git a/tests/test_isfc.py b/tests/test_isfc.py new file mode 100644 index 000000000..561acba83 --- /dev/null +++ b/tests/test_isfc.py @@ -0,0 +1,38 @@ +import brainiak.isfc +import numpy as np + +def test_create_event_segmentation(): + es = brainiak.eventseg.event.EventSegment(5) + assert es, "Invalid EventSegment instance" + + +def test_fit_shapes(): + K = 5 + V = 3 + T = 10 + es = brainiak.eventseg.event.EventSegment(K, n_iter=2) + sample_data = np.random.rand(V, T) + es.fit(sample_data.T) + + assert es.segments_[0].shape == (T,K), "Segmentation from fit " \ + "has incorrect shape" + assert np.isclose(np.sum(es.segments_[0], axis=1), np.ones(T)).all(), \ + "Segmentation from learn_events not correctly normalized" + + T2 = 15 + sample_data2 = np.random.rand(V, T2) + test_segments, test_ll = es.find_events(sample_data2.T) + + assert test_segments.shape == (T2, K), "Segmentation from find_events " \ + "has incorrect shape" + assert np.isclose(np.sum(test_segments, axis=1), np.ones(T2)).all(), \ + "Segmentation from find_events not correctly normalized" + +def test_simple_boundary(): + es = brainiak.eventseg.event.EventSegment(2, n_iter=10) + sample_data = np.asarray([[1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 1]]) + es.fit(sample_data.T) + + events = np.argmax(es.segments_[0], axis=1) + assert np.array_equal(events, [0, 0, 0, 1, 1, 1, 1]),\ + "Failed to correctly segment two events" From d0bd0c4205cf061ee31a46b70d051ee7676a4342 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 21 Feb 2017 11:35:08 -0500 Subject: [PATCH 02/11] Full ISFC code prior to pr checks --- examples/isfc/isfc.py | 17 +++++++++-------- tests/isfc/mask.nii.gz | Bin 0 -> 87 bytes tests/isfc/subj1.nii.gz | Bin 0 -> 124 bytes tests/isfc/subj2.nii.gz | Bin 0 -> 124 bytes tests/isfc/test_isfc.py | 41 ++++++++++++++++++++++++++++++++++++++++ tests/test_isfc.py | 38 ------------------------------------- 6 files changed, 50 insertions(+), 46 deletions(-) create mode 100644 tests/isfc/mask.nii.gz create mode 100644 tests/isfc/subj1.nii.gz create mode 100644 tests/isfc/subj2.nii.gz create mode 100644 tests/isfc/test_isfc.py delete mode 100644 tests/test_isfc.py diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py index ceab42690..23336ca22 100644 --- a/examples/isfc/isfc.py +++ b/examples/isfc/isfc.py @@ -11,11 +11,10 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -"""Intersubject analyses (ISC/ISFC) +"""Example of intersubject analyses (ISC/ISFC) -Functions for computing intersubject correlation (ISC) and intersubject -functional correlation (ISFC), and utility functions for loading subject data -files and masks +Computes ISC for all voxels within a brain mask, and +computes ISFC for voxels with high ISC """ # Authors: Christopher Baldassano and Mor Regev @@ -27,14 +26,16 @@ from matplotlib import pyplot as plt from scipy.cluster.hierarchy import fcluster, linkage import sys, os -dir = os.path.dirname(__file__) +curr_dir = os.path.dirname(__file__) -brain_fname = os.path.join(dir,'avg152T1_gray_3mm.nii.gz') -fnames = [os.path.join(dir,'sub-0' + format(subj, '02') + '-task-intact1.nii.gz') for +brain_fname = os.path.join(curr_dir,'avg152T1_gray_3mm.nii.gz') +fnames = [os.path.join(curr_dir, + 'sub-0' + format(subj, '02') + '-task-intact1.nii.gz') for subj in np.arange(1, 5)] print('Loading data from ', len(fnames), ' subjects...') -D, coords = brainiak.isfc.load_subjects_nii(fnames, brain_fname, lambda x: x > 50) +D, coords = brainiak.isfc.load_subjects_nii(fnames, brain_fname, + lambda x: x > 50) print('Calculating ISC on ', D.shape[0], ' voxels') ISC = brainiak.isfc.isc(D) diff --git a/tests/isfc/mask.nii.gz b/tests/isfc/mask.nii.gz new file mode 100644 index 0000000000000000000000000000000000000000..83b31ffbdea8d48fd5c9c745520306fde30fac26 GIT binary patch literal 87 zcmV-d0I2^TiwFp;X{=ZR|7~G&Yc6hSX#k61WFQJKGq8XmBLf7YYGPp!01GtOGce%8 t8=N8XsJaJ^$ArgfkU!7?v0@AidD@0x%|r$-P%Uj#0sybv>&4&!004TWAtnF- literal 0 HcmV?d00001 diff --git a/tests/isfc/subj1.nii.gz b/tests/isfc/subj1.nii.gz new file mode 100644 index 0000000000000000000000000000000000000000..9cf26af98981dfc65b707115007cfe881f6cd8a0 GIT binary patch literal 124 zcmV-?0E7P@iwFqsSgcqA|8sR>YB4TuX=wn9VPqf*urNR%D+41KqH1Dc5C97_*fTKT z!yB9-@~FB8j>m+@YLGwB0kL8X40+myV9k)w8Yq&Lfe9L4gk3-jl`IU*MB7grl`KTt ePrOP}LJF2&2!}uMDhZo~Pyqm%9mAC10{{SDpDw)s literal 0 HcmV?d00001 diff --git a/tests/isfc/subj2.nii.gz b/tests/isfc/subj2.nii.gz new file mode 100644 index 0000000000000000000000000000000000000000..ec6bb9a44eb799b2f30af4babaeceb26f91dfa43 GIT binary patch literal 124 zcmV-?0E7P@iwFqvSgcqA|8sR>YBDZvX=wn9VPqf*urNR%D+41KqH1Dc5C97_*fTKT z!yB9-@~FB8j>m+@YLGwB0kL8X40+myV9k)w8Yq&Lfe9L4gk3-jl`IU*MB7grl`KTt ePrOP}LJF2&2!}uMDhZo~Pyqm%9mAC10{{SFbuPgG literal 0 HcmV?d00001 diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py new file mode 100644 index 000000000..053d0c41e --- /dev/null +++ b/tests/isfc/test_isfc.py @@ -0,0 +1,41 @@ +import brainiak.isfc +import numpy as np +import sys, os + +def test_ISC(): + D = np.zeros((2,5,3)) + D[:,:,0] = \ + [[-0.36225433, -0.43482456, 0.26723158, 0.16461712, -0.37991465], \ + [-0.62305959, -0.46660116, -0.50037994, 1.81083754, 0.23499509]] + D[:,:,1] = \ + [[-0.31484153, -0.47486988, 0.18966625, -0.12568572, -0.40535156], \ + [-0.78433298, -0.35875085, -0.6121344 , 1.68267639, 0.28603493]] + D[:,:,2] = \ + [[-0.26593192, -0.56914734, 0.28397317, 0.30276589, -0.42637472], \ + [-0.67598379, -0.51549055, -0.64196342, 1.60686666, 0.04127293]] + + ISC = brainiak.isfc.isc(D) + assert np.isclose(ISC, [0.9540602, 0.99585304]).all(), \ + "Calculated ISC does not match ground truth" + +def test_loading_and_ISFC(): + curr_dir = os.path.dirname(__file__) + + mask_fname = os.path.join(curr_dir,'mask.nii.gz') + fnames = [os.path.join(curr_dir,'subj1.nii.gz'), + os.path.join(curr_dir,'subj2.nii.gz')] + + D, coords = brainiak.isfc.load_subjects_nii(fnames, mask_fname) + + assert D.shape == (4,5,2), "Loaded data has incorrect shape" + + ISFC = brainiak.isfc.isfc(D) + + gt = \ + [[ 1, 1, 0, -1], + [ 1, 1, 0, -1], + [ 0, 0, 1, 0], + [-1, -1, 0, 1]] + + assert np.isclose(ISFC, gt).all(), \ + "Calculated ISFC does not match ground truth" diff --git a/tests/test_isfc.py b/tests/test_isfc.py deleted file mode 100644 index 561acba83..000000000 --- a/tests/test_isfc.py +++ /dev/null @@ -1,38 +0,0 @@ -import brainiak.isfc -import numpy as np - -def test_create_event_segmentation(): - es = brainiak.eventseg.event.EventSegment(5) - assert es, "Invalid EventSegment instance" - - -def test_fit_shapes(): - K = 5 - V = 3 - T = 10 - es = brainiak.eventseg.event.EventSegment(K, n_iter=2) - sample_data = np.random.rand(V, T) - es.fit(sample_data.T) - - assert es.segments_[0].shape == (T,K), "Segmentation from fit " \ - "has incorrect shape" - assert np.isclose(np.sum(es.segments_[0], axis=1), np.ones(T)).all(), \ - "Segmentation from learn_events not correctly normalized" - - T2 = 15 - sample_data2 = np.random.rand(V, T2) - test_segments, test_ll = es.find_events(sample_data2.T) - - assert test_segments.shape == (T2, K), "Segmentation from find_events " \ - "has incorrect shape" - assert np.isclose(np.sum(test_segments, axis=1), np.ones(T2)).all(), \ - "Segmentation from find_events not correctly normalized" - -def test_simple_boundary(): - es = brainiak.eventseg.event.EventSegment(2, n_iter=10) - sample_data = np.asarray([[1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 1]]) - es.fit(sample_data.T) - - events = np.argmax(es.segments_[0], axis=1) - assert np.array_equal(events, [0, 0, 0, 1, 1, 1, 1]),\ - "Failed to correctly segment two events" From 10ec5729ecbdb5c59f811a2d38e7f4c6ef23a0ac Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 21 Feb 2017 12:19:30 -0500 Subject: [PATCH 03/11] Full isfc commit --- brainiak/__init__.py | 1 - brainiak/isfc.py | 10 ++++------ examples/isfc/isfc.py | 4 +++- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/brainiak/__init__.py b/brainiak/__init__.py index e25d80792..10e76d6a9 100644 --- a/brainiak/__init__.py +++ b/brainiak/__init__.py @@ -14,7 +14,6 @@ """Brain Imaging Analysis Kit.""" import sys -from . import isfc if sys.version_info < (3, 4): raise Exception( diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 33e5bf29b..97246dbc4 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -25,8 +25,6 @@ from scipy import stats import nibabel as nib import numpy as np -from matplotlib import pyplot as plt -from scipy.cluster.hierarchy import fcluster, linkage def isc(D, collapse_subj=True): @@ -68,10 +66,10 @@ def isc(D, collapse_subj=True): def isfc(D, collapse_subj=True): """Intersubject functional correlation - Computes the correlation between the timecoure of each voxel in each subject - with the average of all other subjects' timecourses in *all* voxels. By - default the result is averaged across subjects, unless collapse_subj is set - to False. + Computes the correlation between the timecoure of each voxel in each + subject with the average of all other subjects' timecourses in *all* + voxels. By default the result is averaged across subjects, unless + collapse_subj is set to False. Uses the high performance compute_correlation routine from fcma.util diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py index 23336ca22..b4bf5c7a3 100644 --- a/examples/isfc/isfc.py +++ b/examples/isfc/isfc.py @@ -20,12 +20,13 @@ # Authors: Christopher Baldassano and Mor Regev # Princeton University, 2017 -import brainiak +import brainiak.isfc import nibabel as nib import numpy as np from matplotlib import pyplot as plt from scipy.cluster.hierarchy import fcluster, linkage import sys, os + curr_dir = os.path.dirname(__file__) brain_fname = os.path.join(curr_dir,'avg152T1_gray_3mm.nii.gz') @@ -41,6 +42,7 @@ ISC = brainiak.isfc.isc(D) ISC[np.isnan(ISC)] = 0 +print('Writing ISC map to file...') brain_nii = nib.load(brain_fname) ISC_vol = np.zeros(brain_nii.shape) ISC_vol[coords] = ISC From 298001e0cbf57f584b6555ee0d129961f61be21a Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 21 Feb 2017 13:01:37 -0500 Subject: [PATCH 04/11] Adding script to download example data --- examples/isfc/download_data.sh | 4 ++++ 1 file changed, 4 insertions(+) create mode 100755 examples/isfc/download_data.sh diff --git a/examples/isfc/download_data.sh b/examples/isfc/download_data.sh new file mode 100755 index 000000000..19d0f05be --- /dev/null +++ b/examples/isfc/download_data.sh @@ -0,0 +1,4 @@ +#!/bin/sh +curl --location -o isfc.zip https://api.onedrive.com/v1.0/shares/s!Aobi2ryypFQCgqQOm2Zhs-TSLu9leQ/root/content +unzip -qo isfc.zip +rm -f isfc.zip From 931b0a07acd8e4e065a7dd8fb6c8e5370e55fe05 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:26:32 -0500 Subject: [PATCH 05/11] Minor changes in response to comments --- brainiak/isfc.py | 24 +++++++++++++++++++----- examples/isfc/isfc.py | 8 ++++---- tests/isfc/test_isfc.py | 4 ++-- 3 files changed, 25 insertions(+), 11 deletions(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 97246dbc4..9730e622f 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -16,13 +16,21 @@ Functions for computing intersubject correlation (ISC) and intersubject functional correlation (ISFC), and utility functions for loading subject data files and masks + +Paper references: +ISC: Hasson, U., Nir, Y., Levy, I., Fuhrmann, G. & Malach, R. Intersubject +synchronization of cortical activity during natural vision. Science 303, +1634–1640 (2004). + +ISFC: Simony E, Honey CJ, Chen J, Lositsky O, Yeshurun Y, Wiesel A, Hasson U +(2016) Dynamic reconfiguration of the default mode network during narrative +comprehension. Nat Commun 7. """ # Authors: Christopher Baldassano and Mor Regev # Princeton University, 2017 from brainiak.fcma.util import compute_correlation -from scipy import stats import nibabel as nib import numpy as np @@ -52,11 +60,12 @@ def isc(D, collapse_subj=True): n_subj = D.shape[2] ISC = np.zeros((n_vox, n_subj)) + # Loop across choice of leave-one-out subject for loo_subj in range(n_subj): group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) subj = D[:, :, loo_subj] for v in range(n_vox): - ISC[v, loo_subj] = stats.pearsonr(group[v, :], subj[v, :])[0] + ISC[v, loo_subj] = compute_correlation(group[v, :], subj[v, :]) if collapse_subj: ISC = np.mean(ISC, axis=1) @@ -91,10 +100,14 @@ def isfc(D, collapse_subj=True): n_vox = D.shape[0] n_subj = D.shape[2] ISFC = np.zeros((n_vox, n_vox, n_subj)) + + # Loop across choice of leave-one-out subject for loo_subj in range(D.shape[2]): group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) subj = D[:, :, loo_subj] ISFC[:, :, loo_subj] = compute_correlation(group, subj) + + # Symmetrize matrix ISFC[:, :, loo_subj] = (ISFC[:, :, loo_subj] + ISFC[:, :, loo_subj].T) / 2 @@ -121,9 +134,10 @@ def load_subjects_nii(data_files, mask_file, mask_func=None): Returns ------- - (D, coords) : tuple of voxel by time by subject ndarray, and - x,y,x (as provided by nibabel) coordinates of voxel locations (as a - tuple of ndarrays + D : voxel by time by subject ndarray + all data within mask, from all subjects + coords : tuple of 3 ndarrays + x,y,z (as provided by nibabel) coordinates of mask voxel locations """ mask_nii = nib.load(mask_file) diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py index b4bf5c7a3..cee151f89 100644 --- a/examples/isfc/isfc.py +++ b/examples/isfc/isfc.py @@ -25,7 +25,8 @@ import numpy as np from matplotlib import pyplot as plt from scipy.cluster.hierarchy import fcluster, linkage -import sys, os +import sys +import os curr_dir = os.path.dirname(__file__) @@ -46,20 +47,19 @@ brain_nii = nib.load(brain_fname) ISC_vol = np.zeros(brain_nii.shape) ISC_vol[coords] = ISC -ISC_nifti = nib.Nifti1Image(ISC_vol, - brain_nii.get_affine(), brain_nii.header) +ISC_nifti = nib.Nifti1Image(ISC_vol, brain_nii.affine(), brain_nii.header) nib.save(ISC_nifti, 'ISC.nii.gz') ISC_mask = ISC > 0.2 print('Calculating ISFC on ', np.sum(ISC_mask), ' voxels...') D_masked = D[ISC_mask, :, :] ISFC = brainiak.isfc.isfc(D_masked) -coords_ISFC = [x[ISC_mask] for x in coords] print('Clustering ISFC...') Z = linkage(ISFC, 'ward') z = fcluster(Z, 2, criterion='maxclust') clust_inds = np.argsort(z) +# Show the ISFC matrix, sorted to show the two main clusters plt.imshow(ISFC[np.ix_(clust_inds,clust_inds)]) plt.show() diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py index 053d0c41e..c16d66b8d 100644 --- a/tests/isfc/test_isfc.py +++ b/tests/isfc/test_isfc.py @@ -31,11 +31,11 @@ def test_loading_and_ISFC(): ISFC = brainiak.isfc.isfc(D) - gt = \ + ground_truth = \ [[ 1, 1, 0, -1], [ 1, 1, 0, -1], [ 0, 0, 1, 0], [-1, -1, 0, 1]] - assert np.isclose(ISFC, gt).all(), \ + assert np.isclose(ISFC, ground_truth).all(), \ "Calculated ISFC does not match ground truth" From a04f2f7ce55e208b3b45e54a7fd4dabef3c5cfd8 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:29:22 -0500 Subject: [PATCH 06/11] Updated description of nii files --- brainiak/isfc.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 9730e622f..7b720b6d4 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -127,8 +127,12 @@ def load_subjects_nii(data_files, mask_file, mask_func=None): Parameters ---------- data_files : list of filenames of subject nii files + each nii file should be 4D (space+time), with all dimensions identical + for all subjects mask_file : filename of mask nii + mask should be 3D (space), with dimensions identical to the first + three dimensions of the data files mask_func : Callable[[ndarray], bool] : default x>0 From 3d62874c5cf1498b012af1ce0a64dabdcc120caf Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:35:43 -0500 Subject: [PATCH 07/11] Fixed bug in example --- examples/isfc/isfc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py index cee151f89..9d0bfd2f1 100644 --- a/examples/isfc/isfc.py +++ b/examples/isfc/isfc.py @@ -47,7 +47,7 @@ brain_nii = nib.load(brain_fname) ISC_vol = np.zeros(brain_nii.shape) ISC_vol[coords] = ISC -ISC_nifti = nib.Nifti1Image(ISC_vol, brain_nii.affine(), brain_nii.header) +ISC_nifti = nib.Nifti1Image(ISC_vol, brain_nii.affine, brain_nii.header) nib.save(ISC_nifti, 'ISC.nii.gz') ISC_mask = ISC > 0.2 From c811e7c1f146692480e9ea0756eaf10e21fb040c Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:39:59 -0500 Subject: [PATCH 08/11] Fixed bug in calling compute_correlation --- brainiak/isfc.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 7b720b6d4..80f23bf66 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -65,7 +65,9 @@ def isc(D, collapse_subj=True): group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) subj = D[:, :, loo_subj] for v in range(n_vox): - ISC[v, loo_subj] = compute_correlation(group[v, :], subj[v, :]) + ISC[v, loo_subj] = compute_correlation( + np.expand_dims(group[v, :], 0), + np.expand_dims(subj[v, :],0)) if collapse_subj: ISC = np.mean(ISC, axis=1) From eec71b2fa7a1b01d7de747bc927b18f1c99d24ed Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:50:13 -0500 Subject: [PATCH 09/11] Fixed whitespace --- brainiak/isfc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 80f23bf66..103fb8383 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -67,7 +67,7 @@ def isc(D, collapse_subj=True): for v in range(n_vox): ISC[v, loo_subj] = compute_correlation( np.expand_dims(group[v, :], 0), - np.expand_dims(subj[v, :],0)) + np.expand_dims(subj[v, :], 0)) if collapse_subj: ISC = np.mean(ISC, axis=1) From 3d4608f38c68198ace0fb957b6b2c1d5935fa829 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 22 Feb 2017 14:52:23 -0500 Subject: [PATCH 10/11] Fixed whitespace --- brainiak/isfc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 80f23bf66..103fb8383 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -67,7 +67,7 @@ def isc(D, collapse_subj=True): for v in range(n_vox): ISC[v, loo_subj] = compute_correlation( np.expand_dims(group[v, :], 0), - np.expand_dims(subj[v, :],0)) + np.expand_dims(subj[v, :], 0)) if collapse_subj: ISC = np.mean(ISC, axis=1) From 7ed3f65dba5876f499fa40aef92fb3ad3beedc40 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Fri, 24 Feb 2017 14:42:18 -0500 Subject: [PATCH 11/11] Reverting correlation computation to use scipy --- brainiak/isfc.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 103fb8383..dca2a8670 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -33,6 +33,7 @@ from brainiak.fcma.util import compute_correlation import nibabel as nib import numpy as np +from scipy import stats def isc(D, collapse_subj=True): @@ -65,9 +66,7 @@ def isc(D, collapse_subj=True): group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2) subj = D[:, :, loo_subj] for v in range(n_vox): - ISC[v, loo_subj] = compute_correlation( - np.expand_dims(group[v, :], 0), - np.expand_dims(subj[v, :], 0)) + ISC[v, loo_subj] = stats.pearsonr(group[v, :], subj[v, :])[0] if collapse_subj: ISC = np.mean(ISC, axis=1)