From 3cf7daea8dcc667684b793d9bffff552627aa886 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 6 Jan 2017 16:37:49 -0800 Subject: [PATCH 01/10] add MVPA voxel selection using Searchlight --- brainiak/fcma/mvpa_voxelselector.py | 178 ++++++++++++++++++++++++++ examples/fcma/file_io.py | 36 +++++- examples/fcma/mvpa_voxel_selection.py | 74 +++++++++++ 3 files changed, 287 insertions(+), 1 deletion(-) create mode 100644 brainiak/fcma/mvpa_voxelselector.py create mode 100644 examples/fcma/mvpa_voxel_selection.py diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py new file mode 100644 index 000000000..cfff512ae --- /dev/null +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -0,0 +1,178 @@ +# Copyright 2016 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. +"""Full Correlation Matrix Analysis (FCMA) + +This implementation is based on the following publications: + +.. [Wang2015] Full correlation matrix analysis (FCMA): An unbiased method for + task-related functional connectivity", + Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. + Journal of Neuroscience Methods, 2015. +""" + +# Authors: Yida Wang +# (Intel Labs), 2017 + +import numpy as np +from sklearn import model_selection +from scipy.stats.mstats import zscore +import logging +from mpi4py import MPI +from brainiak.searchlight.searchlight import Searchlight + +logger = logging.getLogger(__name__) + +__all__ = [ + "MVPAVoxelSelector", +] + + +class MVPAVoxelSelector: + """Activity-based voxel selection component of FCMA + + Parameters + ---------- + + raw\_data: list of 4D array + Assumption: 1. all activity data contains the same number of voxels + 2. the activity data has been z-scored, + ready to compute correlation as matrix multiplication + 3. all subjects have the same number of epochs + 4. epochs belonging to the same subject are adjacent + in the list + 5. voxel selection is always done in the auto-correlation, + i.e. raw_data correlate with themselves + 6. if MPI jobs are running on multiple nodes, the user + home directory is shared by all nodes + + mask: 3D array + + num\_voxels: int + the number of processed voxels defined by mask + + epoch\_info: list of tuple (label, sid, start, end) + the condition labels of the epochs + len(labels) labels equals the number of epochs + assuming the epochs of the same sid are adjacent, + and sid is the index of raw\_data + + num\_folds: int + the number of folds to be conducted in the cross validation + + sl: Searchlight + the distributed Searchlight object + """ + def __init__(self, + raw_data, + mask, + epoch_info, + num_folds, + sl_rad=2, + max_blk_edge=10 + ): + self.raw_data = raw_data + self.mask = mask + self.num_voxels = np.sum(mask==1) + self.epoch_info = epoch_info + self.num_folds = num_folds + self.sl = Searchlight(sl_rad=sl_rad, max_blk_edge=max_blk_edge) + if self.num_voxels == 0: + raise ValueError('Zero processed voxels') + + def _preprocess_data(self): + """ process the raw data according to epoch info + + This is done in rank 0 which has the raw_data read in + Average the activity within epochs and z-scoring within subject. + Write the results to self.processed_data, + which is a 4D array of averaged epoch by epoch processed data + Also write the labels to self.label as a 1D numpy array + """ + num_epochs = len(self.epoch_info) + (d1, d2, d3, _) = self.raw_data[0].shape + self.processed_data = np.empty([d1, d2, d3, num_epochs]) + self.labels = np.empty(num_epochs) + subject_count = [0] # counting the epochs per subject for z-scoring + cur_sid = -1 + # averaging + for idx, epoch in enumerate(self.epoch_info): + self.labels[idx] = epoch[0] + if cur_sid != epoch[1]: + subject_count.append(0) + cur_sid = epoch[1] + subject_count[-1] += 1 + print(self.processed_data.shape, epoch) + self.processed_data[:, :, :, idx] = \ + np.mean(self.raw_data[cur_sid][:, :, :, epoch[2]:epoch[3]], axis=3) + # z-scoring + cur_epoch = 0 + for i in subject_count: + if i > 1: + self.processed_data[:, :, :, cur_epoch:cur_epoch + i] = \ + zscore(self.processed_data[:, :, :, cur_epoch:cur_epoch + i], + axis=3, ddof=0) + cur_epoch += i + # if zscore fails (standard deviation is zero), set all values to be zero + self.processed_data = np.nan_to_num(self.processed_data) + + def run(self, clf): + """ run activity-based voxel selection + + Sort the voxels based on the cross-validation accuracy + of their activity vectors within the searchlight + + Parameters + ---------- + clf: classification function + the classifier to be used in cross validation + + Returns + ------- + results: list of tuple (voxel_id, accuracy) + the accuracy numbers of all voxels, in accuracy descending order + the length of array equals the number of voxels + """ + logger.info( + 'running activity-based voxel selection via Searchlight' + ) + rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + self._preprocess_data() + self.sl.distribute([self.processed_data], self.mask) + self.sl.broadcast(self.labels) + logger.info( + 'data preparation done' + ) + # Searchlight kernel function + def _sfn(l, mask, myrad, bcast_var): + data = l[0][mask,:].T + skf = model_selection.StratifiedKFold(n_splits=self.num_folds, + shuffle=False) + accuracy = np.mean(model_selection.cross_val_score(clf, data, + y=bcast_var, + cv=skf, n_jobs=1)) + return accuracy + # obtain a 3D array with accuracy numbers + result_volume = self.sl.run_searchlight(_sfn) + # get result tuple list from the volume + result_list = result_volume[self.mask] + results = [] + for idx, value in enumerate(result_list): + results.append((idx, value)) + # Sort the voxels + results.sort(key=lambda tup: tup[1], reverse=True) + logger.info( + 'activity-based voxel selection via Searchlight is done' + ) + return results diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py index 7dfe25f00..5111a3a15 100644 --- a/examples/fcma/file_io.py +++ b/examples/fcma/file_io.py @@ -126,7 +126,6 @@ def separate_epochs(activity_data, epoch_list): ) return raw_data, labels - def prepare_data(data_dir, extension, mask_file, epoch_file): """ read the data in and generate epochs of interests, then broadcast to all workers @@ -178,3 +177,38 @@ def prepare_data(data_dir, extension, mask_file, epoch_file): (time2 - time1) ) return raw_data, labels + +def generate_epochs_info(epoch_list): + """ use epoch_list to generate epoch_info defined below + + Parameters + ---------- + epoch\_list: list of 3D (binary) array in shape [condition, nEpochs, nTRs] + specification of epochs and conditions + assuming all subjects have the same number of epochs + len(epoch_list) equals the number of subjects + also assuming an epoch is always a continuous time course + + Returns + ------- + epoch\_info: list of tuple (label, sid, start, end) + the condition labels of the epochs + len(labels) labels equals the number of epochs + assuming the epochs of the same sid are adjacent + """ + time1 = time.time() + epoch_info = [] + for sid, epoch in enumerate(epoch_list): + for cond in range(epoch.shape[0]): + sub_epoch = epoch[cond, :, :] + for eid in range(epoch.shape[1]): + r = np.sum(sub_epoch[eid, :]) + if r > 0: # there is an epoch in this condition + start = np.nonzero(sub_epoch[eid, :])[0][0] + epoch_info.append((cond, sid, start, start+r)) + time2 = time.time() + logger.info( + 'epoch separation done, takes %.2f s' % + (time2 - time1) + ) + return epoch_info diff --git a/examples/fcma/mvpa_voxel_selection.py b/examples/fcma/mvpa_voxel_selection.py new file mode 100644 index 000000000..62576ebc1 --- /dev/null +++ b/examples/fcma/mvpa_voxel_selection.py @@ -0,0 +1,74 @@ +# Copyright 2016 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. + +from brainiak.fcma.mvpa_voxelselector import MVPAVoxelSelector +from sklearn import svm +import sys +import os +from mpi4py import MPI +import logging +from file_io import generate_epochs_info +import nibabel as nib +import numpy as np + +format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s' +# if want to output log to a file instead of outputting log to the console, +# replace "stream=sys.stdout" with "filename='fcma.log'" +logging.basicConfig(level=logging.INFO, format=format, stream=sys.stdout) +logger = logging.getLogger(__name__) + +""" +example running command: +mpirun -np 2 python mvpa_voxel_selection.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz + data/fs_epoch_labels.npy 18 +""" +if __name__ == '__main__': + if MPI.COMM_WORLD.Get_rank()==0: + logger.info( + 'programming starts in %d process(es)' % + MPI.COMM_WORLD.Get_size() + ) + data_dir = sys.argv[1] + extension = sys.argv[2] + mask_file = sys.argv[3] + epoch_file = sys.argv[4] + + raw_data = [] + mask = None + epoch_info = None + if MPI.COMM_WORLD.Get_rank()==0: + mask_img = nib.load(mask_file) + mask = mask_img.get_data() + count = 0 + files = [f for f in sorted(os.listdir(data_dir)) + if os.path.isfile(os.path.join(data_dir, f)) + and f.endswith(extension)] + for f in files: + img = nib.load(os.path.join(data_dir, f)) + data = img.get_data() + raw_data.append(data) + logger.info( + 'file %s is loaded, with data shape %s' % (f, data.shape) + ) + epoch_list = np.load(epoch_file) + epoch_info = generate_epochs_info(epoch_list) + num_subjs = int(sys.argv[5]) + mvs = MVPAVoxelSelector(raw_data, mask, epoch_info, num_subjs) + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=10 + clf = svm.SVC(kernel='rbf', shrinking=False, C=10) + results = mvs.run(clf) + # this output is just for result checking + if MPI.COMM_WORLD.Get_rank()==0: + print(results[0:100]) From 956aa8cb191806b8e44ed078ab5e4cf60e90fe47 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 6 Jan 2017 18:52:06 -0800 Subject: [PATCH 02/10] make the mask binary --- brainiak/fcma/mvpa_voxelselector.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py index cfff512ae..59d02be13 100644 --- a/brainiak/fcma/mvpa_voxelselector.py +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -82,7 +82,7 @@ def __init__(self, max_blk_edge=10 ): self.raw_data = raw_data - self.mask = mask + self.mask = mask.astype(np.bool) self.num_voxels = np.sum(mask==1) self.epoch_info = epoch_info self.num_folds = num_folds @@ -112,7 +112,6 @@ def _preprocess_data(self): subject_count.append(0) cur_sid = epoch[1] subject_count[-1] += 1 - print(self.processed_data.shape, epoch) self.processed_data[:, :, :, idx] = \ np.mean(self.raw_data[cur_sid][:, :, :, epoch[2]:epoch[3]], axis=3) # z-scoring @@ -143,20 +142,23 @@ def run(self, clf): the accuracy numbers of all voxels, in accuracy descending order the length of array equals the number of voxels """ - logger.info( - 'running activity-based voxel selection via Searchlight' - ) rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + logger.info( + 'running activity-based voxel selection via Searchlight' + ) if rank == 0: self._preprocess_data() self.sl.distribute([self.processed_data], self.mask) self.sl.broadcast(self.labels) - logger.info( - 'data preparation done' - ) + if rank == 0: + logger.info( + 'data preparation done' + ) # Searchlight kernel function def _sfn(l, mask, myrad, bcast_var): - data = l[0][mask,:].T + data = l[0][mask, :].T + #print(l[0].shape, mask.shape, data.shape, bcast_var) skf = model_selection.StratifiedKFold(n_splits=self.num_folds, shuffle=False) accuracy = np.mean(model_selection.cross_val_score(clf, data, @@ -172,7 +174,8 @@ def _sfn(l, mask, myrad, bcast_var): results.append((idx, value)) # Sort the voxels results.sort(key=lambda tup: tup[1], reverse=True) - logger.info( - 'activity-based voxel selection via Searchlight is done' - ) + if rank == 0: + logger.info( + 'activity-based voxel selection via Searchlight is done' + ) return results From 01e5f1ca043de669e48756ecbe91df3208a0d2b7 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 9 Jan 2017 14:39:02 -0800 Subject: [PATCH 03/10] getting rid of using object properties in the closure; instead passing them in via bcast_var --- brainiak/fcma/mvpa_voxelselector.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py index 59d02be13..8a64c7ca0 100644 --- a/brainiak/fcma/mvpa_voxelselector.py +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -26,6 +26,7 @@ import numpy as np from sklearn import model_selection +from sklearn import svm from scipy.stats.mstats import zscore import logging from mpi4py import MPI @@ -150,7 +151,7 @@ def run(self, clf): if rank == 0: self._preprocess_data() self.sl.distribute([self.processed_data], self.mask) - self.sl.broadcast(self.labels) + self.sl.broadcast((self.labels, self.num_folds)) if rank == 0: logger.info( 'data preparation done' @@ -159,10 +160,10 @@ def run(self, clf): def _sfn(l, mask, myrad, bcast_var): data = l[0][mask, :].T #print(l[0].shape, mask.shape, data.shape, bcast_var) - skf = model_selection.StratifiedKFold(n_splits=self.num_folds, + skf = model_selection.StratifiedKFold(n_splits=bcast_var[1], shuffle=False) accuracy = np.mean(model_selection.cross_val_score(clf, data, - y=bcast_var, + y=bcast_var[0], cv=skf, n_jobs=1)) return accuracy # obtain a 3D array with accuracy numbers From b6468cda8821166aea872740720172c4405c7b2e Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 10 Jan 2017 12:45:15 -0800 Subject: [PATCH 04/10] read mask in all MPI processes; create the necessary class properties in all MPI processes; write out results in nifti format --- .gitignore | 1 + brainiak/fcma/mvpa_voxelselector.py | 38 ++++++++++++++++++--------- examples/fcma/file_io.py | 6 +++++ examples/fcma/mvpa_voxel_selection.py | 13 +++++---- 4 files changed, 41 insertions(+), 17 deletions(-) diff --git a/.gitignore b/.gitignore index e0035356c..c57a1ba9f 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ *.o *.npy *.log +*.nii.gz brainiak.egg-info build diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py index 8a64c7ca0..528804b37 100644 --- a/brainiak/fcma/mvpa_voxelselector.py +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -59,9 +59,6 @@ class MVPAVoxelSelector: mask: 3D array - num\_voxels: int - the number of processed voxels defined by mask - epoch\_info: list of tuple (label, sid, start, end) the condition labels of the epochs len(labels) labels equals the number of epochs @@ -73,6 +70,13 @@ class MVPAVoxelSelector: sl: Searchlight the distributed Searchlight object + + processed\_data: 4D array in shape [brain 3D + epoch] + contains the averaged and normalized brain data epoch by epoch + it is generated in _\preprocess\_data method + + labels: 1D array + contains the labels of the epochs, extracted from epoch\_info """ def __init__(self, raw_data, @@ -84,11 +88,12 @@ def __init__(self, ): self.raw_data = raw_data self.mask = mask.astype(np.bool) - self.num_voxels = np.sum(mask==1) self.epoch_info = epoch_info self.num_folds = num_folds self.sl = Searchlight(sl_rad=sl_rad, max_blk_edge=max_blk_edge) - if self.num_voxels == 0: + self.processed_data = None + self.labels = None + if np.sum(self.mask==True) == 0: raise ValueError('Zero processed voxels') def _preprocess_data(self): @@ -100,6 +105,10 @@ def _preprocess_data(self): which is a 4D array of averaged epoch by epoch processed data Also write the labels to self.label as a 1D numpy array """ + logger.info( + 'mask size: %d' % + np.sum(self.mask==True) + ) num_epochs = len(self.epoch_info) (d1, d2, d3, _) = self.raw_data[0].shape self.processed_data = np.empty([d1, d2, d3, num_epochs]) @@ -139,6 +148,8 @@ def run(self, clf): Returns ------- + result_volume: 3D array of accuracy numbers + contains the voxelwise accuracy numbers obtained via Searchlight results: list of tuple (voxel_id, accuracy) the accuracy numbers of all voxels, in accuracy descending order the length of array equals the number of voxels @@ -159,24 +170,27 @@ def run(self, clf): # Searchlight kernel function def _sfn(l, mask, myrad, bcast_var): data = l[0][mask, :].T - #print(l[0].shape, mask.shape, data.shape, bcast_var) + #print(l[0].shape, mask.shape, data.shape) skf = model_selection.StratifiedKFold(n_splits=bcast_var[1], shuffle=False) accuracy = np.mean(model_selection.cross_val_score(clf, data, y=bcast_var[0], - cv=skf, n_jobs=1)) + cv=skf, + n_jobs=1)) return accuracy # obtain a 3D array with accuracy numbers result_volume = self.sl.run_searchlight(_sfn) # get result tuple list from the volume result_list = result_volume[self.mask] results = [] - for idx, value in enumerate(result_list): - results.append((idx, value)) - # Sort the voxels - results.sort(key=lambda tup: tup[1], reverse=True) if rank == 0: + for idx, value in enumerate(result_list): + if value is None: + value = 0 + results.append((idx, value)) + # Sort the voxels + results.sort(key=lambda tup: tup[1], reverse=True) logger.info( 'activity-based voxel selection via Searchlight is done' ) - return results + return result_volume, results diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py index 5111a3a15..6af9d34ec 100644 --- a/examples/fcma/file_io.py +++ b/examples/fcma/file_io.py @@ -13,6 +13,7 @@ # limitations under the License. import nibabel as nib +from nibabel.nifti1 import Nifti1Pair import os import math import time @@ -212,3 +213,8 @@ def generate_epochs_info(epoch_list): (time2 - time1) ) return epoch_info + + +def write_nifti_file(data, affine, filename): + img = Nifti1Pair(data, affine) + nib.nifti1.save(img, filename) diff --git a/examples/fcma/mvpa_voxel_selection.py b/examples/fcma/mvpa_voxel_selection.py index 62576ebc1..59565187a 100644 --- a/examples/fcma/mvpa_voxel_selection.py +++ b/examples/fcma/mvpa_voxel_selection.py @@ -19,6 +19,7 @@ from mpi4py import MPI import logging from file_io import generate_epochs_info +from file_io import write_nifti_file import nibabel as nib import numpy as np @@ -45,11 +46,11 @@ epoch_file = sys.argv[4] raw_data = [] - mask = None + # all MPI processes read the mask; the mask file is small + mask_img = nib.load(mask_file) + mask = mask_img.get_data() epoch_info = None if MPI.COMM_WORLD.Get_rank()==0: - mask_img = nib.load(mask_file) - mask = mask_img.get_data() count = 0 files = [f for f in sorted(os.listdir(data_dir)) if os.path.isfile(os.path.join(data_dir, f)) @@ -68,7 +69,9 @@ # for cross validation, use SVM with precomputed kernel # no shrinking, set C=10 clf = svm.SVC(kernel='rbf', shrinking=False, C=10) - results = mvs.run(clf) + # only rank 0 has meaningful return values + result_volume, results = mvs.run(clf) # this output is just for result checking if MPI.COMM_WORLD.Get_rank()==0: - print(results[0:100]) + result_volume = np.nan_to_num(result_volume.astype(np.float)) + write_nifti_file(result_volume, mask_img.affine, 'result.nii.gz') From 7ec086ebf54e1124028aa0952de6ab9e71b270cc Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 10 Jan 2017 14:14:45 -0800 Subject: [PATCH 05/10] add mvpa searchlight test case --- brainiak/.DS_Store | Bin 0 -> 6148 bytes brainiak/fcma/.DS_Store | Bin 0 -> 6148 bytes examples/.DS_Store | Bin 0 -> 10244 bytes examples/fcma/.DS_Store | Bin 0 -> 6148 bytes examples/fcma/data/.DS_Store | Bin 0 -> 6148 bytes examples/fcma/model/logistic.pkl | Bin 0 -> 846 bytes examples/fcma/model/svm.pkl | Bin 0 -> 1409 bytes examples/fcma/mvpa_voxel_selection.py | 1 - tests/fcma/test_mvpa_voxel_selection.py | 61 ++++++++++++++++++++++++ 9 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 brainiak/.DS_Store create mode 100644 brainiak/fcma/.DS_Store create mode 100644 examples/.DS_Store create mode 100644 examples/fcma/.DS_Store create mode 100644 examples/fcma/data/.DS_Store create mode 100644 examples/fcma/model/logistic.pkl create mode 100644 examples/fcma/model/svm.pkl create mode 100644 tests/fcma/test_mvpa_voxel_selection.py diff --git a/brainiak/.DS_Store b/brainiak/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5a24c85f99298ad19b10ebb5b42c339b1d2faabf GIT binary patch literal 6148 zcmeHKO>fgc5S?w);-r9#0I3oeA3=&x359SBp&)TUiU?|jC<=CC8cVKi#ZICqij*@~ zE6)`ZA=?NjqYkd?l$`K36(#~E0{>eC_}g8R-)uk;Eq=e_hheraf|@J@UDo0_$6WWx;C*jj&Rr!2P_xcwF#5PcLQb8O79x6+uJvnc`j3 z0{44?TLO2F4#>^>$&nnwR~#Lo#KJfC^G0c+NBq-E!XV8yHhzeu<wDU6xmj;+KZ)b;*mk|`z#nz&OZVa^bA!n5FHm=Pf~Xsh zn(-j&rn_y`-cHcXdL8?6tM3m%ZMl)Z9i)L5gh4j8%6xU~^~t1Gt5u!ao%*!uOzQQG zsd3iyG(f6~=KuhV6b7#?e0N z*uRWG4Uhs(c!&h-Q&b2or~7bWlAOb|7-KW=k-P*fPT>JknNXHz%ET`+S7i?>-Ni}5 ziC*+Aqg@|u_T{uVz?+5n);ZCUEQ0VPg8oU*+L#5!E!v?USai$P85c4A-!(a*%fqxo- EU*VVM-T(jq literal 0 HcmV?d00001 diff --git a/brainiak/fcma/.DS_Store b/brainiak/fcma/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..7ace9f0ff99e00be5e351899d7212e167c115af2 GIT binary patch literal 6148 zcmeHK%}&BV5S|qhiGXtSSQ8IkAR)$kL(s&lJ^(_gfuyvgpq5)7!-w(#d>j2{b|w9x z(ReV#>@>SyyPa=m=WDy$B_cI=n)Zo0MASfGtRJD;VSJxU$qr1*F;Gb#b6V2~M~{{z zZ`)uMunO#(0{rbZ<(r*TNZIdqa~nk7B#0*P;OXH8<>EsOXJ4!(ZOEfB%8Y#Or@(qE zu-fnp;J?B#r$-c@Vwo=r@zh;`*#a{HcO)1KE5tlgVo95K)9Bi;_CWcZrNkeM6I zXfyCb#GROt_sH9kWtMqX&Eyep`|Q#sUD1NNG=(NHy+Dr;qgzpwk05GfBb-!5;Qo^W zn>ER@&!N1BxMO7HQ-bmMh|5oeGev&C8zn2K8W{x-CC0i!Z@e2h?>1sbP+{@yGse4- zW46GGpxzJOjU2Croa55d71%1k^MixJ=o(Bl%G7~MwgA97n#FL=KMTxp47vtWjkp6-nhMlZWwsbf z)6wo2Ki6QYQPWA8&4)6PmD!;viH`mqaVOUt~P7*|!tOq-(9%h)Bz>;1lu3l`!PP+45izN-+cIxY|fNw_8 z%vrMqaZWM)lgsDa7a6ifB97vu?KUN~o3(URi@EhK#Tm-hyDiJe)sTk}?g;9K=4VFCtKBbU~asK_`R6EP(b))XH~X2^m6XCDGw2EEVR$1iCUw zoXM8V@yGKVFQTr7dMW1W1iITOu%PbFwPxW5Ej-nv~=uC96CIDB=g`ixHU{U zCzbJbt~`OYX|`+F!(q;0-^J{9ujaNpr&ucSZ|mH){az_5S5QvvNIJS@+PzxJMj4gM zY$SI_+%jFx^~EhCoiOdFTpm}9c-l}Lt<^9h`e66fp<`2q z`ghq@$_cZ{`gNS2YsqO_+o|EC`molHQ`<6#XG#|?UB2>;VAEZ#o5#!K@(O95w4iU8 z65r5zM(CLsQ0%y7s=7J23r~va==*goC6@bwVz;RzwDF3{MT?hI`Td33*;)G2G5tUq z50P9g2;&uUd0AkYU)rPM7i_mwI-4rLBZzFIW|`l=n~sE28kVYSP`k=6b=$aAS#cFt zuaTtg2k;onw#ZEz*Gm2#Ly2jItcrY{<_(g+ZCF0%dCqdFq zI>|mfB%|aQIZj?AZ zLeL6Zp%?m~9~9URL!iU`V8cN;45M%a9)U;U7(5P7z;Sp1PQocT4X?vlcn{u(58xxX z1fRnfa2dXV@8Jjd5q^R{V8U1AE8(B%%U(=fyqzt5rC6Ju^X|Rvd;0pEuiA3YmX57A zG3T$8QlihDH-EvxK+UQ(Yd35R=Lk_ifqX790&<2i=4F}1TtA7<Z(?Y$bESRA8n8t1!TEAlaH>KnngSZ zauFZBQ)&^A?((ehC5skUZVE|~Oz3eHx%P^xV)5Tn_%m{ed`o^t@?Il~Q9 zum(0T+J>PWHp33+f}Ie79@q;}BybE=NPr1RIDjNR0vUJ+#^7Nj@nej{Pr?a!7M_FW zk<72b8F&NUgtuHWe+cK{A`<$m0_t{VNy{aCJw0^`eqkpYdMytD4*?GW4*?H>TO5HI z{QCp+{r~Ri|NnoB@6!8x4*?H>ze50u+9U04sF4lx3YFf$CXja literal 0 HcmV?d00001 diff --git a/examples/fcma/.DS_Store b/examples/fcma/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..fae1e99e9d2b2734783e89a01c43772fa54ab15b GIT binary patch literal 6148 zcmeHK&u`N(6n@@1$T|&W{E#Mb$q}S!+RC6DOBqNUkS2;YA%;rQhBZCzNf_4b~7*9a2kk~o_N{VQe{bfA-$at_|A*qFD6&4fv;z>q-{#bwg1dE>; ziyYNepWIAOcq~UqVlY6Ct#7F3b>gNT^KY*i263{s_Djsn=H@P3wC3{*{MbwR7nO$W zR1!C7$7Q$mT)hQ#-_zwSH)#yZPo)fx%dWQ__?@bK^}dV}H;DY!i0ST55Y=U;EZb2% z-mRqmc7l4+sM^=ITK)k{x7^6z4&uNI!XP=c=J~1G%iV6FP*`yacZR&zikA@7{m-`04YPuiw;=Nnp!Xaz*11zN2tV$U&tQg)-Vh9qU_| z=o;;V)KyHrq`eHWJjF*NL6TR25vl`$17k{9sA(RxnS+0c3Gpc$_bzyPM>Q7_HXK0v z2;WmfiGx>IBmKHG#qxxpB8+Lw_A(uy)z{oVv>Pk*(f%o+&tn|(fq|U8=m4-5&b}+A z2ZjN|z&T@p_Xh`wv81t|qS`u;sVe|5hiWA#^Ung?I2ubD`zhiJM5Ly`)KuszhDgm( z?`XY}#(s*aISGCF5V~cdZzw{w4*!l)Cs9h#8s7>|AUk7|KmpH%`jjX z_^%icxy{OE8AH-%>%`#rtaXrHBC)VuKSgDNOdrRpz(?^ik`jD#xBx6^?5Btti1`sv MG?>gV@J|`|1AIdJGynhq literal 0 HcmV?d00001 diff --git a/examples/fcma/data/.DS_Store b/examples/fcma/data/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..f9c0b52b2bd73f49a24bd25c9793fd8eab3fcb98 GIT binary patch literal 6148 zcmeHKu}Z{15PhRp9M~M#SZ--0_y@~z3f6vrnp3d|Mi2|Tbbj{vW@Z&b6s#=d4a~gk z?7VF9icAK8tzT|Wffaxy-4Smd=H~nEGuvw!Bhq<~Ck)u)0axf})!!$adx(2HlXExW z-=Jfr=gtWpx$}C_58JLEIx_rm&*{f>+wB--Qa}nw0VyB_q`>bK@ZL*X9waJC0VyB_ zJ{9oqL!mp?hGSxUIv8REATF2=<8{mu#Nr8JZ8#<}L$jn3lWJ|mu%t6zT3u~8CMF#g zC-Xk-WNQxTn#-TFwAi3Fkz;LKn}g2n?`ozSdA=sCFOPS4BXoP z#;@SYg`dD$VUkOHuvfD7d+e@$)Z)l6p#&R*Nu04Tmdsj_2yTM*B81A}PpB5K(jr@x zN<|GbYFU;`ZlTusMQXSg%LGWjX4ybIBpb^HoVZ06j6xm`IP1&My{ay<+Cwt`rWWRN6oN_a6tiE zGq^-Zu2{yQLG}$6GH`1N7nflZcPAK#g^W0Im>GN+DaC>?Fd!K)2 zkzd99Hrzcs{fhxUWEEe%A(OCwMYb_cgrjo;xOYh!xar-|F&F*R!9ka-1CKp{VXCM7dB@>Nm+T_o zTm`i;!;$F2Q^YQ*6%|+k3eUXz=LNi2?*4;!|8n8J)M>Ru>ChsZxWI{+4I3Q#PG)T4 zvEfyx=c#g>0O_Xinyi$^PvDJL?zDuXv(x>LzyBPxkFcdp{L6D4reC)`4BmEqeac%f(i1K%9Dn;%A4cInMe75P literal 0 HcmV?d00001 diff --git a/examples/fcma/model/svm.pkl b/examples/fcma/model/svm.pkl new file mode 100644 index 0000000000000000000000000000000000000000..7ff52141daffef3b2a9cfd860150100f5ad81ef8 GIT binary patch literal 1409 zcmZWpS##7z5ME;&Y+-H^ha)7MSs+ntZbQJ>1`HP7kZ@RUQ={40)v!9IHPT{Jr7C!# z@K52N@Zcx@3Gb;QJ-rKZA-(KO?|fhP*Vmg-5oBB?iuc@(h`B2wZj99RJ@IgAR^izLqySFYzWOL%15muV2Gz)gx+Z!oO{K5$z17P5?Q ztYutlq=zx&94}#K)1N@jj*4W$V+p(by?8bK*wwv=&CI%q-hc_;K_ic%^c9O^u3_>$ zp-d4_#NY_>emRh^PZo-y*cM1dzsg`eMm}Js8B^681I8fl5$-kUD3u zKeX{3DB)l;{4S(R6w;w^@_p!p*3S;~lN4}x8Gq>`GQmYCnNhDL{NR1`r|_SQC+#$5 zCg-_?A6x!zw6S)j=*H_AHp`3pnfV!gF6{ zL24wN+oUTqnrC?ZZE%&^s3=XKS>@MlfnH4->e3AEig>=_2L z1iVzjm(tk!}ZEQwYXtld)EJBGpEgnm-OKQ;7Vp8u`TcC@SjJEsMv>8aNk zMx6#cP?5?v;Mpb)zI}*$hZ&CjFAY5Rw)}l)lR6clVYue7Siwtp;Y|%{hj07ox{M60 Kt 0: + output.append(int(8*tuple[1])) + expected_output = [6, 6, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2] + assert np.allclose(output, expected_output, atol=1), \ + 'voxel selection via SVM does not provide correct results' + +if __name__ == '__main__': + test_mvpa_voxel_selection() From c107e8e7fe568a487f45e7ad89061ca6c835874f Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 10 Jan 2017 14:16:23 -0800 Subject: [PATCH 06/10] remove unncessary files --- brainiak/.DS_Store | Bin 6148 -> 0 bytes brainiak/fcma/.DS_Store | Bin 6148 -> 0 bytes examples/.DS_Store | Bin 10244 -> 0 bytes examples/fcma/.DS_Store | Bin 6148 -> 0 bytes examples/fcma/data/.DS_Store | Bin 6148 -> 0 bytes examples/fcma/model/logistic.pkl | Bin 846 -> 0 bytes examples/fcma/model/svm.pkl | Bin 1409 -> 0 bytes 7 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 brainiak/.DS_Store delete mode 100644 brainiak/fcma/.DS_Store delete mode 100644 examples/.DS_Store delete mode 100644 examples/fcma/.DS_Store delete mode 100644 examples/fcma/data/.DS_Store delete mode 100644 examples/fcma/model/logistic.pkl delete mode 100644 examples/fcma/model/svm.pkl diff --git a/brainiak/.DS_Store b/brainiak/.DS_Store deleted file mode 100644 index 5a24c85f99298ad19b10ebb5b42c339b1d2faabf..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKO>fgc5S?w);-r9#0I3oeA3=&x359SBp&)TUiU?|jC<=CC8cVKi#ZICqij*@~ zE6)`ZA=?NjqYkd?l$`K36(#~E0{>eC_}g8R-)uk;Eq=e_hheraf|@J@UDo0_$6WWx;C*jj&Rr!2P_xcwF#5PcLQb8O79x6+uJvnc`j3 z0{44?TLO2F4#>^>$&nnwR~#Lo#KJfC^G0c+NBq-E!XV8yHhzeu<wDU6xmj;+KZ)b;*mk|`z#nz&OZVa^bA!n5FHm=Pf~Xsh zn(-j&rn_y`-cHcXdL8?6tM3m%ZMl)Z9i)L5gh4j8%6xU~^~t1Gt5u!ao%*!uOzQQG zsd3iyG(f6~=KuhV6b7#?e0N z*uRWG4Uhs(c!&h-Q&b2or~7bWlAOb|7-KW=k-P*fPT>JknNXHz%ET`+S7i?>-Ni}5 ziC*+Aqg@|u_T{uVz?+5n);ZCUEQ0VPg8oU*+L#5!E!v?USai$P85c4A-!(a*%fqxo- EU*VVM-T(jq diff --git a/brainiak/fcma/.DS_Store b/brainiak/fcma/.DS_Store deleted file mode 100644 index 7ace9f0ff99e00be5e351899d7212e167c115af2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHK%}&BV5S|qhiGXtSSQ8IkAR)$kL(s&lJ^(_gfuyvgpq5)7!-w(#d>j2{b|w9x z(ReV#>@>SyyPa=m=WDy$B_cI=n)Zo0MASfGtRJD;VSJxU$qr1*F;Gb#b6V2~M~{{z zZ`)uMunO#(0{rbZ<(r*TNZIdqa~nk7B#0*P;OXH8<>EsOXJ4!(ZOEfB%8Y#Or@(qE zu-fnp;J?B#r$-c@Vwo=r@zh;`*#a{HcO)1KE5tlgVo95K)9Bi;_CWcZrNkeM6I zXfyCb#GROt_sH9kWtMqX&Eyep`|Q#sUD1NNG=(NHy+Dr;qgzpwk05GfBb-!5;Qo^W zn>ER@&!N1BxMO7HQ-bmMh|5oeGev&C8zn2K8W{x-CC0i!Z@e2h?>1sbP+{@yGse4- zW46GGpxzJOjU2Croa55d71%1k^MixJ=o(Bl%G7~MwgA97n#FL=KMTxp47vtWjkp6-nhMlZWwsbf z)6wo2Ki6QYQPWA8&4)6PmD!;viH`mqaVOUt~P7*|!tOq-(9%h)Bz>;1lu3l`!PP+45izN-+cIxY|fNw_8 z%vrMqaZWM)lgsDa7a6ifB97vu?KUN~o3(URi@EhK#Tm-hyDiJe)sTk}?g;9K=4VFCtKBbU~asK_`R6EP(b))XH~X2^m6XCDGw2EEVR$1iCUw zoXM8V@yGKVFQTr7dMW1W1iITOu%PbFwPxW5Ej-nv~=uC96CIDB=g`ixHU{U zCzbJbt~`OYX|`+F!(q;0-^J{9ujaNpr&ucSZ|mH){az_5S5QvvNIJS@+PzxJMj4gM zY$SI_+%jFx^~EhCoiOdFTpm}9c-l}Lt<^9h`e66fp<`2q z`ghq@$_cZ{`gNS2YsqO_+o|EC`molHQ`<6#XG#|?UB2>;VAEZ#o5#!K@(O95w4iU8 z65r5zM(CLsQ0%y7s=7J23r~va==*goC6@bwVz;RzwDF3{MT?hI`Td33*;)G2G5tUq z50P9g2;&uUd0AkYU)rPM7i_mwI-4rLBZzFIW|`l=n~sE28kVYSP`k=6b=$aAS#cFt zuaTtg2k;onw#ZEz*Gm2#Ly2jItcrY{<_(g+ZCF0%dCqdFq zI>|mfB%|aQIZj?AZ zLeL6Zp%?m~9~9URL!iU`V8cN;45M%a9)U;U7(5P7z;Sp1PQocT4X?vlcn{u(58xxX z1fRnfa2dXV@8Jjd5q^R{V8U1AE8(B%%U(=fyqzt5rC6Ju^X|Rvd;0pEuiA3YmX57A zG3T$8QlihDH-EvxK+UQ(Yd35R=Lk_ifqX790&<2i=4F}1TtA7<Z(?Y$bESRA8n8t1!TEAlaH>KnngSZ zauFZBQ)&^A?((ehC5skUZVE|~Oz3eHx%P^xV)5Tn_%m{ed`o^t@?Il~Q9 zum(0T+J>PWHp33+f}Ie79@q;}BybE=NPr1RIDjNR0vUJ+#^7Nj@nej{Pr?a!7M_FW zk<72b8F&NUgtuHWe+cK{A`<$m0_t{VNy{aCJw0^`eqkpYdMytD4*?GW4*?H>TO5HI z{QCp+{r~Ri|NnoB@6!8x4*?H>ze50u+9U04sF4lx3YFf$CXja diff --git a/examples/fcma/.DS_Store b/examples/fcma/.DS_Store deleted file mode 100644 index fae1e99e9d2b2734783e89a01c43772fa54ab15b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHK&u`N(6n@@1$T|&W{E#Mb$q}S!+RC6DOBqNUkS2;YA%;rQhBZCzNf_4b~7*9a2kk~o_N{VQe{bfA-$at_|A*qFD6&4fv;z>q-{#bwg1dE>; ziyYNepWIAOcq~UqVlY6Ct#7F3b>gNT^KY*i263{s_Djsn=H@P3wC3{*{MbwR7nO$W zR1!C7$7Q$mT)hQ#-_zwSH)#yZPo)fx%dWQ__?@bK^}dV}H;DY!i0ST55Y=U;EZb2% z-mRqmc7l4+sM^=ITK)k{x7^6z4&uNI!XP=c=J~1G%iV6FP*`yacZR&zikA@7{m-`04YPuiw;=Nnp!Xaz*11zN2tV$U&tQg)-Vh9qU_| z=o;;V)KyHrq`eHWJjF*NL6TR25vl`$17k{9sA(RxnS+0c3Gpc$_bzyPM>Q7_HXK0v z2;WmfiGx>IBmKHG#qxxpB8+Lw_A(uy)z{oVv>Pk*(f%o+&tn|(fq|U8=m4-5&b}+A z2ZjN|z&T@p_Xh`wv81t|qS`u;sVe|5hiWA#^Ung?I2ubD`zhiJM5Ly`)KuszhDgm( z?`XY}#(s*aISGCF5V~cdZzw{w4*!l)Cs9h#8s7>|AUk7|KmpH%`jjX z_^%icxy{OE8AH-%>%`#rtaXrHBC)VuKSgDNOdrRpz(?^ik`jD#xBx6^?5Btti1`sv MG?>gV@J|`|1AIdJGynhq diff --git a/examples/fcma/data/.DS_Store b/examples/fcma/data/.DS_Store deleted file mode 100644 index f9c0b52b2bd73f49a24bd25c9793fd8eab3fcb98..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKu}Z{15PhRp9M~M#SZ--0_y@~z3f6vrnp3d|Mi2|Tbbj{vW@Z&b6s#=d4a~gk z?7VF9icAK8tzT|Wffaxy-4Smd=H~nEGuvw!Bhq<~Ck)u)0axf})!!$adx(2HlXExW z-=Jfr=gtWpx$}C_58JLEIx_rm&*{f>+wB--Qa}nw0VyB_q`>bK@ZL*X9waJC0VyB_ zJ{9oqL!mp?hGSxUIv8REATF2=<8{mu#Nr8JZ8#<}L$jn3lWJ|mu%t6zT3u~8CMF#g zC-Xk-WNQxTn#-TFwAi3Fkz;LKn}g2n?`ozSdA=sCFOPS4BXoP z#;@SYg`dD$VUkOHuvfD7d+e@$)Z)l6p#&R*Nu04Tmdsj_2yTM*B81A}PpB5K(jr@x zN<|GbYFU;`ZlTusMQXSg%LGWjX4ybIBpb^HoVZ06j6xm`IP1&My{ay<+Cwt`rWWRN6oN_a6tiE zGq^-Zu2{yQLG}$6GH`1N7nflZcPAK#g^W0Im>GN+DaC>?Fd!K)2 zkzd99Hrzcs{fhxUWEEe%A(OCwMYb_cgrjo;xOYh!xar-|F&F*R!9ka-1CKp{VXCM7dB@>Nm+T_o zTm`i;!;$F2Q^YQ*6%|+k3eUXz=LNi2?*4;!|8n8J)M>Ru>ChsZxWI{+4I3Q#PG)T4 zvEfyx=c#g>0O_Xinyi$^PvDJL?zDuXv(x>LzyBPxkFcdp{L6D4reC)`4BmEqeac%f(i1K%9Dn;%A4cInMe75P diff --git a/examples/fcma/model/svm.pkl b/examples/fcma/model/svm.pkl deleted file mode 100644 index 7ff52141daffef3b2a9cfd860150100f5ad81ef8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1409 zcmZWpS##7z5ME;&Y+-H^ha)7MSs+ntZbQJ>1`HP7kZ@RUQ={40)v!9IHPT{Jr7C!# z@K52N@Zcx@3Gb;QJ-rKZA-(KO?|fhP*Vmg-5oBB?iuc@(h`B2wZj99RJ@IgAR^izLqySFYzWOL%15muV2Gz)gx+Z!oO{K5$z17P5?Q ztYutlq=zx&94}#K)1N@jj*4W$V+p(by?8bK*wwv=&CI%q-hc_;K_ic%^c9O^u3_>$ zp-d4_#NY_>emRh^PZo-y*cM1dzsg`eMm}Js8B^681I8fl5$-kUD3u zKeX{3DB)l;{4S(R6w;w^@_p!p*3S;~lN4}x8Gq>`GQmYCnNhDL{NR1`r|_SQC+#$5 zCg-_?A6x!zw6S)j=*H_AHp`3pnfV!gF6{ zL24wN+oUTqnrC?ZZE%&^s3=XKS>@MlfnH4->e3AEig>=_2L z1iVzjm(tk!}ZEQwYXtld)EJBGpEgnm-OKQ;7Vp8u`TcC@SjJEsMv>8aNk zMx6#cP?5?v;Mpb)zI}*$hZ&CjFAY5Rw)}l)lR6clVYue7Siwtp;Y|%{hj07ox{M60 Kt Date: Tue, 10 Jan 2017 14:16:55 -0800 Subject: [PATCH 07/10] more gitignore patterns --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index c57a1ba9f..9746561ff 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ *.npy *.log *.nii.gz +*.pkl brainiak.egg-info build @@ -22,6 +23,7 @@ htmlcov .ipynb_checkpoints junit-*.xml __pycache__ +.DS_Store brainiak/fcma/cython_blas.c brainiak/eventseg/_utils.c From b32887c984af527d87dbab2808ee1f9216f48382 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 10 Jan 2017 16:55:52 -0800 Subject: [PATCH 08/10] fix format issues revealed by pr-check --- brainiak/fcma/__init__.py | 19 ++++++++++++++++++- brainiak/fcma/classifier.py | 7 +------ brainiak/fcma/mvpa_voxelselector.py | 27 +++++++++++++-------------- brainiak/fcma/voxelselector.py | 16 +--------------- examples/fcma/file_io.py | 11 +++++++++++ 5 files changed, 44 insertions(+), 36 deletions(-) diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py index 5f9a9900f..a3485d6d8 100644 --- a/brainiak/fcma/__init__.py +++ b/brainiak/fcma/__init__.py @@ -11,4 +11,21 @@ # 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. -"""Full correlation matrix analysis""" +"""Full correlation matrix analysis + +The implementation is based on the following publications: + +.. [Wang2015-1] Full correlation matrix analysis (FCMA): An unbiased method for + task-related functional connectivity", + Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. + Journal of Neuroscience Methods, 2015. + +.. [Wang2015-2] "Full correlation matrix analysis of fMRI data on Intel® Xeon + Phi™ coprocessors", + Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, + Kai Li, Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, + Theodore L. Willke. + In Proceedings of the International Conference for + High Performance Computing, + Networking, Storage and Analysis. 2015. +""" diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index c2781ba30..1dc489b9c 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -13,12 +13,7 @@ # limitations under the License. """Full Correlation Matrix Analysis (FCMA) -This implementation is based on the following publications: - -.. [Wang2015] Full correlation matrix analysis (FCMA): An unbiased method for - task-related functional connectivity", - Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. - Journal of Neuroscience Methods, 2015. +Corrlation-based training and predicrtion """ # Authors: Yida Wang diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py index 528804b37..a9b8ca54d 100644 --- a/brainiak/fcma/mvpa_voxelselector.py +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -13,12 +13,7 @@ # limitations under the License. """Full Correlation Matrix Analysis (FCMA) -This implementation is based on the following publications: - -.. [Wang2015] Full correlation matrix analysis (FCMA): An unbiased method for - task-related functional connectivity", - Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. - Journal of Neuroscience Methods, 2015. +Activity-based voxel selction """ # Authors: Yida Wang @@ -26,7 +21,6 @@ import numpy as np from sklearn import model_selection -from sklearn import svm from scipy.stats.mstats import zscore import logging from mpi4py import MPI @@ -93,7 +87,8 @@ def __init__(self, self.sl = Searchlight(sl_rad=sl_rad, max_blk_edge=max_blk_edge) self.processed_data = None self.labels = None - if np.sum(self.mask==True) == 0: + num_voxels = np.sum(self.mask) + if num_voxels == 0: raise ValueError('Zero processed voxels') def _preprocess_data(self): @@ -107,13 +102,13 @@ def _preprocess_data(self): """ logger.info( 'mask size: %d' % - np.sum(self.mask==True) + np.sum(self.mask) ) num_epochs = len(self.epoch_info) (d1, d2, d3, _) = self.raw_data[0].shape self.processed_data = np.empty([d1, d2, d3, num_epochs]) self.labels = np.empty(num_epochs) - subject_count = [0] # counting the epochs per subject for z-scoring + subject_count = [0] # counting the epochs per subject for z-scoring cur_sid = -1 # averaging for idx, epoch in enumerate(self.epoch_info): @@ -123,16 +118,19 @@ def _preprocess_data(self): cur_sid = epoch[1] subject_count[-1] += 1 self.processed_data[:, :, :, idx] = \ - np.mean(self.raw_data[cur_sid][:, :, :, epoch[2]:epoch[3]], axis=3) + np.mean(self.raw_data[cur_sid][:, :, :, epoch[2]:epoch[3]], + axis=3) # z-scoring cur_epoch = 0 for i in subject_count: if i > 1: self.processed_data[:, :, :, cur_epoch:cur_epoch + i] = \ - zscore(self.processed_data[:, :, :, cur_epoch:cur_epoch + i], + zscore(self.processed_data[:, :, :, + cur_epoch:cur_epoch + i], axis=3, ddof=0) cur_epoch += i - # if zscore fails (standard deviation is zero), set all values to be zero + # if zscore fails (standard deviation is zero), + # set all values to be zero self.processed_data = np.nan_to_num(self.processed_data) def run(self, clf): @@ -167,10 +165,11 @@ def run(self, clf): logger.info( 'data preparation done' ) + # Searchlight kernel function def _sfn(l, mask, myrad, bcast_var): data = l[0][mask, :].T - #print(l[0].shape, mask.shape, data.shape) + # print(l[0].shape, mask.shape, data.shape) skf = model_selection.StratifiedKFold(n_splits=bcast_var[1], shuffle=False) accuracy = np.mean(model_selection.cross_val_score(clf, data, diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index 31ad4532d..e3891e78e 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -13,21 +13,7 @@ # limitations under the License. """Full Correlation Matrix Analysis (FCMA) -This implementation is based on the following publications: - -.. [Wang2015-1] Full correlation matrix analysis (FCMA): An unbiased method for - task-related functional connectivity", - Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. - Journal of Neuroscience Methods, 2015. - -.. [Wang2015-2] "Full correlation matrix analysis of fMRI data on Intel® Xeon - Phi™ coprocessors", - Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, - Kai Li, Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, - Theodore L. Willke. - In Proceedings of the International Conference for - High Performance Computing, - Networking, Storage and Analysis. 2015. +Correlation-based voxel selection """ # Authors: Yida Wang diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py index 6af9d34ec..d73acc281 100644 --- a/examples/fcma/file_io.py +++ b/examples/fcma/file_io.py @@ -216,5 +216,16 @@ def generate_epochs_info(epoch_list): def write_nifti_file(data, affine, filename): + """ write a nifti file given data and affine + + Parameters + ---------- + data: 3D/4D numpy array + the brain data with/without time dimension + affine: 2D numpy array + affine of the image, usually inherited from an existing image + filename: string + the output filename + """ img = Nifti1Pair(data, affine) nib.nifti1.save(img, filename) From 2798a08d10d94aa63a36d43b729d9bf1dadf7f11 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 10 Jan 2017 21:56:48 -0800 Subject: [PATCH 09/10] install scikit-learn[alldeps] so that it works when no wheel is available --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 33552e9c9..adc9aedc3 100644 --- a/setup.py +++ b/setup.py @@ -124,7 +124,7 @@ def finalize_options(self): 'cython', 'mpi4py', 'numpy', - 'scikit-learn>=0.18', + 'scikit-learn[alldeps]>=0.18', 'scipy', 'pymanopt', 'theano', From db1a1b4b3404987513c692cc210e31fed43b5f5b Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 11 Jan 2017 16:39:50 -0800 Subject: [PATCH 10/10] address review comments --- .gitignore | 3 ++ brainiak/fcma/classifier.py | 2 +- brainiak/fcma/mvpa_voxelselector.py | 60 ++++++++++------------- examples/fcma/download_data.sh | 4 ++ examples/fcma/file_io.py | 21 ++++---- examples/fcma/mvpa_voxel_selection.py | 9 ++-- examples/fcma/run_mvpa_voxel_selection.sh | 2 + examples/fcma/run_voxel_selection.sh | 3 ++ examples/fcma/voxel_selection.py | 4 +- tests/fcma/test_mvpa_voxel_selection.py | 4 +- 10 files changed, 60 insertions(+), 52 deletions(-) create mode 100644 examples/fcma/download_data.sh create mode 100644 examples/fcma/run_mvpa_voxel_selection.sh create mode 100644 examples/fcma/run_voxel_selection.sh diff --git a/.gitignore b/.gitignore index 9746561ff..611eed748 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ *.npy *.log *.nii.gz +*.zip *.pkl brainiak.egg-info @@ -24,6 +25,8 @@ htmlcov junit-*.xml __pycache__ .DS_Store +__MACOSX brainiak/fcma/cython_blas.c brainiak/eventseg/_utils.c +brainiak/examples/fcma/face_scene/ diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 1dc489b9c..672e5a294 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -13,7 +13,7 @@ # limitations under the License. """Full Correlation Matrix Analysis (FCMA) -Corrlation-based training and predicrtion +Correlation-based training and prediction """ # Authors: Yida Wang diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py index a9b8ca54d..d1f2ecdb2 100644 --- a/brainiak/fcma/mvpa_voxelselector.py +++ b/brainiak/fcma/mvpa_voxelselector.py @@ -13,7 +13,7 @@ # limitations under the License. """Full Correlation Matrix Analysis (FCMA) -Activity-based voxel selction +Activity-based voxel selection """ # Authors: Yida Wang @@ -24,7 +24,6 @@ from scipy.stats.mstats import zscore import logging from mpi4py import MPI -from brainiak.searchlight.searchlight import Searchlight logger = logging.getLogger(__name__) @@ -40,24 +39,18 @@ class MVPAVoxelSelector: ---------- raw\_data: list of 4D array - Assumption: 1. all activity data contains the same number of voxels - 2. the activity data has been z-scored, - ready to compute correlation as matrix multiplication - 3. all subjects have the same number of epochs - 4. epochs belonging to the same subject are adjacent - in the list - 5. voxel selection is always done in the auto-correlation, - i.e. raw_data correlate with themselves - 6. if MPI jobs are running on multiple nodes, the user - home directory is shared by all nodes + each element of the list is the raw data of a subject, + in the shape of [x, y, z, t] mask: 3D array - epoch\_info: list of tuple (label, sid, start, end) - the condition labels of the epochs - len(labels) labels equals the number of epochs - assuming the epochs of the same sid are adjacent, - and sid is the index of raw\_data + epoch\_info: list of tuple (label, sid, start, end). + label is the condition labels of the epochs; + sid is the subject id, corresponding to the index of raw_data; + start is the start TR of an epoch (inclusive); + end is the end TR of an epoch(exclusive). + Assuming len(labels) labels equals the number of epochs and + the epochs of the same sid are adjacent in epoch_info num\_folds: int the number of folds to be conducted in the cross validation @@ -65,11 +58,11 @@ class MVPAVoxelSelector: sl: Searchlight the distributed Searchlight object - processed\_data: 4D array in shape [brain 3D + epoch] + processed\_data\_: 4D array in shape [brain 3D + epoch] contains the averaged and normalized brain data epoch by epoch - it is generated in _\preprocess\_data method + it is generated from raw\_data and epoch\_info - labels: 1D array + labels\_: 1D array contains the labels of the epochs, extracted from epoch\_info """ def __init__(self, @@ -77,16 +70,15 @@ def __init__(self, mask, epoch_info, num_folds, - sl_rad=2, - max_blk_edge=10 + sl ): self.raw_data = raw_data self.mask = mask.astype(np.bool) self.epoch_info = epoch_info self.num_folds = num_folds - self.sl = Searchlight(sl_rad=sl_rad, max_blk_edge=max_blk_edge) - self.processed_data = None - self.labels = None + self.sl = sl + self.processed_data_ = None + self.labels_ = None num_voxels = np.sum(self.mask) if num_voxels == 0: raise ValueError('Zero processed voxels') @@ -106,32 +98,32 @@ def _preprocess_data(self): ) num_epochs = len(self.epoch_info) (d1, d2, d3, _) = self.raw_data[0].shape - self.processed_data = np.empty([d1, d2, d3, num_epochs]) - self.labels = np.empty(num_epochs) + self.processed_data_ = np.empty([d1, d2, d3, num_epochs]) + self.labels_ = np.empty(num_epochs) subject_count = [0] # counting the epochs per subject for z-scoring cur_sid = -1 # averaging for idx, epoch in enumerate(self.epoch_info): - self.labels[idx] = epoch[0] + self.labels_[idx] = epoch[0] if cur_sid != epoch[1]: subject_count.append(0) cur_sid = epoch[1] subject_count[-1] += 1 - self.processed_data[:, :, :, idx] = \ + self.processed_data_[:, :, :, idx] = \ np.mean(self.raw_data[cur_sid][:, :, :, epoch[2]:epoch[3]], axis=3) # z-scoring cur_epoch = 0 for i in subject_count: if i > 1: - self.processed_data[:, :, :, cur_epoch:cur_epoch + i] = \ - zscore(self.processed_data[:, :, :, + self.processed_data_[:, :, :, cur_epoch:cur_epoch + i] = \ + zscore(self.processed_data_[:, :, :, cur_epoch:cur_epoch + i], axis=3, ddof=0) cur_epoch += i # if zscore fails (standard deviation is zero), # set all values to be zero - self.processed_data = np.nan_to_num(self.processed_data) + self.processed_data_ = np.nan_to_num(self.processed_data_) def run(self, clf): """ run activity-based voxel selection @@ -159,8 +151,8 @@ def run(self, clf): ) if rank == 0: self._preprocess_data() - self.sl.distribute([self.processed_data], self.mask) - self.sl.broadcast((self.labels, self.num_folds)) + self.sl.distribute([self.processed_data_], self.mask) + self.sl.broadcast((self.labels_, self.num_folds)) if rank == 0: logger.info( 'data preparation done' diff --git a/examples/fcma/download_data.sh b/examples/fcma/download_data.sh new file mode 100644 index 000000000..3edd8ffbf --- /dev/null +++ b/examples/fcma/download_data.sh @@ -0,0 +1,4 @@ +#!/bin/sh +curl --location -o face_scene.zip https://www.dropbox.com/s/11adrjdkt0w1tr3/face_scene.zip?dl=0 +unzip -qo face_scene.zip +rm -f face_scene.zip diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py index d73acc281..929bbe07e 100644 --- a/examples/fcma/file_io.py +++ b/examples/fcma/file_io.py @@ -185,17 +185,20 @@ def generate_epochs_info(epoch_list): Parameters ---------- epoch\_list: list of 3D (binary) array in shape [condition, nEpochs, nTRs] - specification of epochs and conditions - assuming all subjects have the same number of epochs - len(epoch_list) equals the number of subjects - also assuming an epoch is always a continuous time course + Contains specification of epochs and conditions, + Assumption: 1. all subjects have the same number of epochs; + 2. len(epoch_list) equals the number of subjects; + 3. an epoch is always a continuous time course. Returns ------- - epoch\_info: list of tuple (label, sid, start, end) - the condition labels of the epochs - len(labels) labels equals the number of epochs - assuming the epochs of the same sid are adjacent + epoch\_info: list of tuple (label, sid, start, end). + label is the condition labels of the epochs; + sid is the subject id, corresponding to the index of raw_data; + start is the start TR of an epoch (inclusive); + end is the end TR of an epoch(exclusive). + Assuming len(labels) labels equals the number of epochs and + the epochs of the same sid are adjacent in epoch_info """ time1 = time.time() epoch_info = [] @@ -208,7 +211,7 @@ def generate_epochs_info(epoch_list): start = np.nonzero(sub_epoch[eid, :])[0][0] epoch_info.append((cond, sid, start, start+r)) time2 = time.time() - logger.info( + logger.debug( 'epoch separation done, takes %.2f s' % (time2 - time1) ) diff --git a/examples/fcma/mvpa_voxel_selection.py b/examples/fcma/mvpa_voxel_selection.py index 8798314c3..38b41b699 100644 --- a/examples/fcma/mvpa_voxel_selection.py +++ b/examples/fcma/mvpa_voxel_selection.py @@ -22,6 +22,7 @@ from file_io import write_nifti_file import nibabel as nib import numpy as np +from brainiak.searchlight.searchlight import Searchlight format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s' # if want to output log to a file instead of outputting log to the console, @@ -30,9 +31,7 @@ logger = logging.getLogger(__name__) """ -example running command: -mpirun -np 2 python mvpa_voxel_selection.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz - data/fs_epoch_labels.npy 18 +example running command in run_mvpa_voxel_selection.sh """ if __name__ == '__main__': if MPI.COMM_WORLD.Get_rank()==0: @@ -64,7 +63,9 @@ epoch_list = np.load(epoch_file) epoch_info = generate_epochs_info(epoch_list) num_subjs = int(sys.argv[5]) - mvs = MVPAVoxelSelector(raw_data, mask, epoch_info, num_subjs) + # create a Searchlight object + sl = Searchlight(sl_rad=2) + mvs = MVPAVoxelSelector(raw_data, mask, epoch_info, num_subjs, sl) # for cross validation, use SVM with precomputed kernel # no shrinking, set C=10 clf = svm.SVC(kernel='rbf', shrinking=False, C=10) diff --git a/examples/fcma/run_mvpa_voxel_selection.sh b/examples/fcma/run_mvpa_voxel_selection.sh new file mode 100644 index 000000000..73521190c --- /dev/null +++ b/examples/fcma/run_mvpa_voxel_selection.sh @@ -0,0 +1,2 @@ +#!/bin/sh +mpirun -np 1 python mvpa_voxel_selection.py face_scene bet.nii.gz face_scene/visual_top_mask.nii.gz face_scene/fs_epoch_labels.npy 18 diff --git a/examples/fcma/run_voxel_selection.sh b/examples/fcma/run_voxel_selection.sh new file mode 100644 index 000000000..a9bc330b2 --- /dev/null +++ b/examples/fcma/run_voxel_selection.sh @@ -0,0 +1,3 @@ +#!/bin/sh + +mpirun -np 2 python voxel_selection.py face_scene bet.nii.gz face_scene/mask.nii.gz face_scene/fs_epoch_labels.npy 12 18 diff --git a/examples/fcma/voxel_selection.py b/examples/fcma/voxel_selection.py index 88a1e9dfc..33ea937b3 100644 --- a/examples/fcma/voxel_selection.py +++ b/examples/fcma/voxel_selection.py @@ -26,9 +26,7 @@ logger = logging.getLogger(__name__) """ -example running command: -mpirun -np 2 python voxel_selection.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz - data/fs_epoch_labels.npy 12 18 +example running command in run_voxel_selection.sh """ if __name__ == '__main__': if MPI.COMM_WORLD.Get_rank()==0: diff --git a/tests/fcma/test_mvpa_voxel_selection.py b/tests/fcma/test_mvpa_voxel_selection.py index ecc07c2fe..68ccf6d92 100644 --- a/tests/fcma/test_mvpa_voxel_selection.py +++ b/tests/fcma/test_mvpa_voxel_selection.py @@ -13,6 +13,7 @@ # limitations under the License. from brainiak.fcma.mvpa_voxelselector import MVPAVoxelSelector +from brainiak.searchlight.searchlight import Searchlight from sklearn import svm import numpy as np from mpi4py import MPI @@ -42,7 +43,8 @@ def test_mvpa_voxel_selection(): epoch_info.append((1, 1, 70, 80)) # 2 subjects, 4 epochs per subject - mvs = MVPAVoxelSelector(raw_data, mask, epoch_info, 2, sl_rad=1) + sl = Searchlight(sl_rad=1) + mvs = MVPAVoxelSelector(raw_data, mask, epoch_info, 2, sl) # for cross validation, use SVM with precomputed kernel clf = svm.SVC(kernel='rbf', C=10)