diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py new file mode 100644 index 000000000..56f79a62e --- /dev/null +++ b/brainiak/fcma/classifier.py @@ -0,0 +1,280 @@ +# 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), 2016 + +import numpy as np +import time +from sklearn.base import BaseEstimator +import sklearn +from . import fcma_extension +from . import cython_blas as blas +import logging + +logger = logging.getLogger(__name__) + +__all__ = [ + "Classifier", +] + + +class Classifier(BaseEstimator): + """Correlation-based classification component of FCMA + + Parameters + ---------- + + clf: class + The classifier used, normally a classifier class of sklearn + + epochs_per_subj: int, default 0 + The number of epochs of each subject + within-subject normalization will be performed during + classifier training if epochs_per_subj is specified + default 0 means no within-subject normalization + + + Attributes + ---------- + + training_data_: 2D numpy array in shape [num_samples, num_features] + default None + training_data is None except clf is SVM.SVC with precomputed kernel, + in which case training data is needed to compute + the similarity vector for each sample to be classified + + num_voxels_: int + The number of voxels per brain used in this classifier + this is defined by the applied mask, normally the top voxels + selected by FCMA voxel selection + num_voxels must be consistent in both training and classification + + num_samples_: int + The number of samples of the training set + """ + def __init__(self, + clf, + epochs_per_subj=0): + self.clf = clf + self.epochs_per_subj = epochs_per_subj + return + + def fit(self, X, y): + """ use correlation data to train a model + + first compute the correlation of the input data, + and then normalize within subject + if more than one sample in one subject, + and then fit to a model defined by self.clf. + + Parameters + ---------- + X: a list of numpy array in shape [num_TRs, num_voxels] + X contains the activity data filtered by top voxels + and prepared for correlation computation. + assuming all elements of X has the same num_voxels value + y: labels, len(X) equals len(Y) + + Returns + ------- + self: return the object itself + """ + time1 = time.time() + assert len(X) == len(y), \ + 'the number of samples does not match the number labels' + num_samples = len(X) + num_voxels = X[0].shape[1] # see assumption above + self.num_voxels_ = num_voxels + self.num_samples_ = num_samples + corr_data = np.zeros((num_samples, num_voxels, num_voxels), + np.float32, order='C') + # compute correlation + count = 0 + for data in X: + num_TRs = data.shape[0] + # syrk performs slower in this case + # blas.compute_single_self_correlation_syrk('L', 'N', + # num_voxels, + # num_TRs, + # 1.0, data, + # num_voxels, 0.0, + # corr_data, + # num_voxels, count) + blas.compute_single_self_correlation_gemm('N', 'T', + num_voxels, + num_voxels, + num_TRs, + 1.0, data, + num_voxels, num_voxels, + 0.0, corr_data, + num_voxels, count) + count += 1 + logger.debug( + 'correlation computation done' + ) + # normalize if necessary + if self.epochs_per_subj > 0: + corr_data = corr_data.reshape(1, + num_samples, + num_voxels * num_voxels) + fcma_extension.normalization(corr_data, self.epochs_per_subj) + corr_data = corr_data.reshape(num_samples, num_voxels, num_voxels) + logger.debug( + 'normalization done' + ) + # training + if isinstance(self.clf, sklearn.svm.SVC) \ + and self.clf.kernel == 'precomputed': + kernel_matrix = np.zeros((num_samples, num_samples), + np.float32, + order='C') + # for using kernel matrix computation from voxel selection + corr_data = corr_data.reshape(1, + num_samples, + num_voxels * num_voxels) + blas.compute_kernel_matrix('L', 'T', + num_samples, num_voxels * num_voxels, + 1.0, corr_data, + 0, num_voxels * num_voxels, + 0.0, kernel_matrix, num_samples) + data = kernel_matrix + # training data is in shape [num_samples, num_voxels * num_voxels] + self.training_data_ = corr_data.reshape(num_samples, + num_voxels * num_voxels) + logger.debug( + 'kernel computation done' + ) + else: + data = corr_data.reshape(num_samples, num_voxels * num_voxels) + self.training_data_ = None + + self.clf = self.clf.fit(data, y) + time2 = time.time() + logger.info( + 'training done, takes %.2f s' % + (time2 - time1) + ) + return self + + def predict(self, X): + """ use a trained model to predict correlation data + + first compute the correlation of the input data, + and then normalize across all samples in the list + if len(X) > 1, + and then predict via self.clf. + + Parameters + ---------- + X: a list of numpy array in shape [num_TRs, self.num_voxels\_] + X contains the activity data filtered by top voxels + and prepared for correlation computation. + len(X) is the number of test samples + if len(X) > 0: normalization is done + on all test samples + + Returns + ------- + y_pred: the predicted label of X, in shape [len(X),] + """ + time1 = time.time() + num_test_samples = len(X) + assert num_test_samples > 0, \ + 'at least one sample is needed' + corr_data = np.zeros((num_test_samples, + self.num_voxels_, + self.num_voxels_), + np.float32, + order='C') + # compute correlation + count = 0 + for data in X: + num_TRs = data.shape[0] + num_voxels = data.shape[1] + assert self.num_voxels_ == num_voxels, \ + 'the number of voxels provided by X does not match ' \ + 'the number of voxels defined in the model' + blas.compute_single_self_correlation_gemm('N', 'T', + num_voxels, + num_voxels, + num_TRs, + 1.0, data, + num_voxels, num_voxels, + 0.0, corr_data, + num_voxels, count) + count += 1 + logger.debug( + 'correlation computation done' + ) + # normalize if necessary + if num_test_samples > 1: + corr_data = corr_data.reshape(1, + num_test_samples, + num_voxels * num_voxels) + fcma_extension.normalization(corr_data, + num_test_samples) + corr_data = corr_data.reshape(num_test_samples, + num_voxels, + num_voxels) + logger.debug( + 'normalization done' + ) + # predict + if isinstance(self.clf, sklearn.svm.SVC) \ + and self.clf.kernel == 'precomputed': + assert self.training_data_ is not None, \ + 'when using precomputed kernel of SVM, ' \ + 'all training data must be provided' + num_training_samples = self.training_data_.shape[0] + data = np.zeros((num_test_samples, num_training_samples), + np.float32, + order='C') + corr_data = corr_data.reshape(num_test_samples, + num_voxels * num_voxels) + # compute the similarity matrix using corr_data and training_data + blas.compute_single_matrix_multiplication('T', 'N', + num_training_samples, + num_test_samples, + num_voxels * num_voxels, + 1.0, + self.training_data_, + num_voxels * num_voxels, + corr_data, + num_voxels * num_voxels, + 0.0, + data, + num_training_samples) + logger.debug( + 'similarity matrix computation done' + ) + else: + data = corr_data.reshape(num_test_samples, + num_voxels * num_voxels) + y_pred = self.clf.predict(data) + time2 = time.time() + logger.info( + 'prediction done, takes %.2f s' % + (time2 - time1) + ) + return y_pred diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 6f4e2ad3f..ffc1617cb 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -12,11 +12,15 @@ # See the License for the specific language governing permissions and # limitations under the License. +# Authors: Yida Wang +# (Intel Labs), 2016 + cimport scipy.linalg.cython_blas as blas -def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a, py_lda, - int py_start_voxel, py_ldb, py_beta, py_c, py_ldc, int py_start_epoch): - """ use blas API wrapped by scipy.linalg.cython_blas to compute correlation +def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, + py_alpha, py_a, py_lda, int py_start_voxel, + py_ldb, py_beta, py_c, py_ldc, int py_start_epoch): + """ use blas API sgemm wrapped by scipy to compute correlation The blas APIs process matrices in column-major, but our matrices are in row-major, @@ -100,14 +104,15 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a cdef float[:, :, ::1] C C = py_c blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda, - &A[0, py_start_voxel], &ldb, &beta, &C[0, py_start_epoch,0], &ldc) + &A[0, py_start_voxel], &ldb, &beta, &C[0, py_start_epoch, 0], &ldc) -def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_start_voxel, py_lda, +def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, + int py_start_voxel, py_lda, py_beta, py_c, py_ldc): - """ use blas API wrapped by scipy.linalg.cython_blas to compute kernel matrix of SVM + """ use blas API syrk wrapped by scipy to compute kernel matrix of SVM - The blas APIs process matrices in column-major, but our matrices are in row-major, - so we play the transpose trick here, i.e. A*B=(B^T*A^T)^T + The blas APIs process matrices in column-major, but our matrices are + in row-major, so we play the transpose trick here, i.e. A*B=(B^T*A^T)^T In SVM with linear kernel, the distance of two samples is essentially the dot product of them. @@ -115,8 +120,8 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ Since the kernel matrix is symmetric, ssyrk is used, the other half of the matrix is assigned later. In our case, the dimension of samples is much larger than - the number samples, so we proportionally shrink the values of the kernel matrix - for getting more robust alpha values in SVM iteration. + the number samples, so we proportionally shrink the values of + the kernel matrix for getting more robust alpha values in SVM iteration. Parameters ---------- @@ -139,8 +144,8 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ py_alpha: float the weight applied to the input matrix A - py_a: 3D array in shape [num_assigned_voxels, num_epochs, num_voxels] in our case - the normalized correlation values of a voxel + py_a: 3D array in shape [num_assigned_voxels, num_epochs, num_voxels] + in our case the normalized correlation values of a voxel py_start_voxel: int the processed voxel @@ -182,13 +187,14 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ C = py_c blas.ssyrk(uplo, trans, &N, &K, &alpha, &A[py_start_voxel, 0, 0], &lda, &beta, &C[0, 0], &ldc) - # shrink the values for getting more stable alpha values in SVM training iteration + # shrink the values for getting more stable alpha values + # in SVM training iteration num_digits = len(str(int(py_c[0, 0]))) - if (num_digits > 2): + if num_digits > 2: proportion = 10**(2-num_digits) py_c *= proportion # complete the other half of the kernel matrix - if (py_uplo=='L'): + if py_uplo == 'L': for j in range(py_c.shape[0]): for k in range(j): py_c[j, k] = py_c[k, j] @@ -197,7 +203,275 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ for k in range(j): py_c[k, j] = py_c[j, k] -def installed(): +def compute_single_self_correlation_syrk(py_uplo, py_trans, py_n, py_k, + py_alpha, py_a, py_lda, + py_beta, py_c, py_ldc, + int py_start_sample): + """ use blas API syrk wrapped by scipy to compute correlation matrix + + This is to compute the correlation between selected voxels for + final training and classification. Since the resulting correlation + matrix is symmetric, syrk is used. However, it looks like that in most + cases, syrk performs much worse than gemm (the next function). + Here we assume that the resulting matrix is stored in a compact way, + i.e. py_ldc == py_n. + + Parameters + ---------- + py_uplo: str + getting the upper or lower triangle of the matrix + + py_trans: str + do transpose or not for the input matrix A + + py_n: int + the row and column of the resulting matrix C + in our case, is num_selected_voxels + + py_k: int + the collapsed dimension of the multiplying matrices + i.e. the column of the first matrix after transpose if necessary + the row of the second matrix after transpose if necessary + in our case, is num_TRs + + py_alpha: float + the weight applied to the input matrix A + + py_a: 2D array in shape [num_TRs, num_selected_voxels] + in our case the normalized activity values + + py_lda: int + the stride of the input matrix A + + py_beta: float + the weight applied to the resulting matrix C + + py_c: 3D array + in shape [num_samples, num_selected_voxels, num_selected_voxels] + place to store the resulting kernel matrix + + py_ldc: int + the stride of the resulting matrix + + py_start_sample: int + the processed sample + used to locate the resulting matrix C + + Returns + ------- + py_c: 3D array + in shape [num_samples, num_selected_voxels, num_selected_voxels] + write the resulting correlation matrices + for the processed sample + """ + cdef bytes by_uplo=py_uplo.encode() + cdef bytes by_trans=py_trans.encode() + cdef char* uplo = by_uplo + cdef char* trans = by_trans + cdef int N, K, lda, ldc + N = py_n + K = py_k + lda = py_lda + ldc = py_ldc + cdef float alpha, beta + alpha = py_alpha + beta = py_beta + cdef float[:, ::1] A + A = py_a + cdef float[:, :, ::1] C + C = py_c + blas.ssyrk(uplo, trans, &N, &K, &alpha, &A[0, 0], &lda, + &beta, &C[py_start_sample, 0, 0], &ldc) + # complete the other half of the kernel matrix + if py_uplo == 'L': + for j in range(py_c.shape[1]): + for k in range(j): + py_c[py_start_sample, j, k] = py_c[py_start_sample, k, j] + else: + for j in range(py_c.shape[1]): + for k in range(j): + py_c[py_start_sample, k, j] = py_c[py_start_sample, j, k] + +def compute_single_self_correlation_gemm(py_trans_a, py_trans_b, py_m, py_n, + py_k, py_alpha, py_a, py_lda, + py_ldb, py_beta, py_c, py_ldc, + int py_start_sample): + """ use blas API gemm wrapped by scipy to compute correlation matrix + + This is to compute the correlation between selected voxels for + final training and classification. Although the resulting correlation + matrix is symmetric, in most cases, gemm performs better than syrk. + Here we assume that the resulting matrix is stored in a compact way, + i.e. py_ldc == py_n. + + Parameters + ---------- + py_trans_a: str + do transpose or not for the first matrix A + + py_trans_b: str + do transpose or not for the first matrix B + + py_m: int + the row of the resulting matrix C + in our case, is num_selected_voxels + + py_n: int + the column of the resulting matrix C + in our case, is num_selected_voxels + + py_k: int + the collapsed dimension of the multiplying matrices + i.e. the column of the first matrix after transpose if necessary + the row of the second matrix after transpose if necessary + in our case, is num_TRs + + py_alpha: float + the weight applied to the input matrix A + + py_a: 2D array in shape [num_TRs, num_selected_voxels] + in our case the normalized activity values + both multipliers are specified here as the same one + + py_lda: int + the stride of the input matrix A + + py_ldb: int + the stride of the input matrix B + in our case, the same as py_lda + + py_beta: float + the weight applied to the resulting matrix C + + py_c: 3D array + in shape [num_samples, num_selected_voxels, num_selected_voxels] + place to store the resulting kernel matrix + + py_ldc: int + the stride of the resulting matrix + + py_start_sample: int + the processed sample + used to locate the resulting matrix C + + Returns + ------- + py_c: 3D array + in shape [num_samples, num_selected_voxels, num_selected_voxels] + write the resulting correlation matrices + for the processed sample """ - This is an empty method for installing cython_blas library + cdef bytes by_trans_a=py_trans_a.encode() + cdef bytes by_trans_b=py_trans_b.encode() + cdef char* trans_a = by_trans_a + cdef char* trans_b = by_trans_b + cdef int M, N, K, lda, ldb, ldc + M = py_m + N = py_n + K = py_k + lda = py_lda + ldb = py_ldb + ldc = py_ldc + cdef float alpha, beta + alpha = py_alpha + beta = py_beta + cdef float[:, ::1] A + A = py_a + cdef float[:, :, ::1] C + C = py_c + blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda, + &A[0, 0], &ldb, &beta, &C[py_start_sample, 0, 0], &ldc) + +def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, + py_k, py_alpha, py_a, py_lda, + py_b, py_ldb, py_beta, py_c, py_ldc): + """ use blas API gemm wrapped by scipy to compute similarity matrix of SVM + + This is to compute the similarity matrix between test sample(s) and + training samples. The blas APIs process matrices in column-major, + but our matrices are in row-major, so we play the transpose trick here, + i.e. A*B=(B^T*A^T)^T + + Parameters + ---------- + py_trans_a: str + do transpose or not for the first matrix A + + py_trans_b: str + do transpose or not for the first matrix B + + py_m: int + the row of the resulting matrix C + in our case, is num_training_samples + + py_n: int + the column of the resulting matrix C + in our case, is num_test_samples + + py_k: int + the collapsed dimension of the multiplying matrices + i.e. the column of the first matrix after transpose if necessary + the row of the second matrix after transpose if necessary + in our case, is num_selected_voxels*num_selected_voxels + + py_alpha: float + the weight applied to the input matrix A + + py_a: 2D array + in shape [num_training_samples, num_selected_voxels*num_selected_voxels] + + py_lda: int + the stride of the input matrix A + + py_b: 2D array + in shape [num_test_samples, num_selected_voxels*num_selected_voxels] + + py_ldb: int + the stride of the input matrix B + + py_beta: float + the weight applied to the resulting matrix C + + py_c: 2D array + in shape [num_training_samples, num_test_samples] of column-major + in fact it is + in shape [num_test_samples, num_training_samples] of row-major + place to store the resulting similarity matrix + + py_ldc: int + the stride of the resulting matrix + + Returns + ------- + py_c: 3D array + in shape [num_training_samples, num_test_samples] of column-major + write the resulting similarity matrix """ + cdef bytes by_trans_a=py_trans_a.encode() + cdef bytes by_trans_b=py_trans_b.encode() + cdef char* trans_a = by_trans_a + cdef char* trans_b = by_trans_b + cdef int M, N, K, lda, ldb, ldc + M = py_m + N = py_n + K = py_k + lda = py_lda + ldb = py_ldb + ldc = py_ldc + cdef float alpha, beta + alpha = py_alpha + beta = py_beta + cdef float[:, ::1] A + A = py_a + cdef float[:, ::1] B + B = py_b + cdef float[:, ::1] C + C = py_c + blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda, + &B[0, 0], &ldb, &beta, &C[0, 0], &ldc) + # shrink the values for getting more stable alpha values + # in SVM training iteration + num_digits = len(str(int(py_c[0, 0]))) + if num_digits > 2: + proportion = 10**(2-num_digits) + py_c *= proportion diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index 432b2343f..1826a5265 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -79,7 +79,7 @@ class VoxelSelector: The number of folds to be conducted in the cross validation voxel_unit: int, default 100 - The number of voxel assigned to a worker each time + The number of voxels assigned to a worker each time master_rank: int, default 0 The process which serves as the master diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py new file mode 100644 index 000000000..7e6896d53 --- /dev/null +++ b/examples/fcma/classification.py @@ -0,0 +1,51 @@ +# 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.classifier import Classifier +from sklearn import svm +#from sklearn.linear_model import LogisticRegression +import sys +import logging +from file_io import prepare_data +import numpy as np +#from sklearn.externals import joblib + +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__) + + +# python classification.py /Users/yidawang/data/face_scene/raw nii.gz +# /Users/yidawang/data/face_scene/prefrontal_top_mask.nii.gz data/fs_epoch_labels.npy 12 +if __name__ == '__main__': + data_dir = sys.argv[1] + extension = sys.argv[2] + mask_file = sys.argv[3] + epoch_file = sys.argv[4] + raw_data, labels = prepare_data(data_dir, extension, mask_file, epoch_file) + epochs_per_subj = int(sys.argv[5]) + # no shrinking, set C=1 + use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + #use_clf = LogisticRegression() + clf = Classifier(use_clf, epochs_per_subj=epochs_per_subj) + training_data = raw_data[0:204] + test_data = raw_data[204:] + clf.fit(training_data, labels[0:204]) + # joblib can be used for saving and loading models + #joblib.dump(clf, 'model/logistic.pkl') + #clf = joblib.load('model/svm.pkl') + print(clf.predict(test_data)) + print(np.asanyarray(labels[204:])) diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py new file mode 100644 index 000000000..7dfe25f00 --- /dev/null +++ b/examples/fcma/file_io.py @@ -0,0 +1,180 @@ +# 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. + +import nibabel as nib +import os +import math +import time +import numpy as np +import logging +from scipy.stats.mstats import zscore +from mpi4py import MPI + +logger = logging.getLogger(__name__) + +def read_activity_data(dir, file_extension, mask_file): + """ read data in NIfTI format and apply the spatial mask to them + + Parameters + ---------- + dir: str + the path to all subject files + file_extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, we apply the mask right after + reading a file for saving memory + + Returns + ------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats + len(activity_data) equals the number of subjects + """ + time1 = time.time() + mask_img = nib.load(mask_file) + mask = mask_img.get_data() + count = 0 + for index in np.ndindex(mask.shape): + if mask[index] != 0: + count += 1 + files = [f for f in sorted(os.listdir(dir)) + if os.path.isfile(os.path.join(dir, f)) + and f.endswith(file_extension)] + activity_data = [] + for f in files: + img = nib.load(os.path.join(dir, f)) + data = img.get_data() + (d1, d2, d3, d4) = data.shape + masked_data = np.zeros([d4, count], np.float32, order='C') + count1 = 0 + for index in np.ndindex(mask.shape): + if mask[index] != 0: + masked_data[:, count1] = np.copy(data[index]) + count1 += 1 + activity_data.append(masked_data) + logger.info( + 'file %s is loaded and masked, with data shape %s' % + (f, masked_data.shape) + ) + time2 = time.time() + logger.info( + 'data reading done, takes %.2f s' % + (time2 - time1) + ) + return activity_data + + +def separate_epochs(activity_data, epoch_list): + """ separate data into epochs of interest specified in epoch_list + and z-score them for computing correlation + + Parameters + ---------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats of all subjects + epoch_list: list of 3D 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 + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + and z-scored in preparation of correlation computation + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs + """ + time1 = time.time() + raw_data = [] + labels = [] + for sid in range(len(epoch_list)): + epoch = epoch_list[sid] + 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 + # mat is row-major + # regardless of the order of acitvity_data[sid] + mat = activity_data[sid][sub_epoch[eid, :] == 1, :] + mat = zscore(mat, axis=0, ddof=0) + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) + mat = mat / math.sqrt(r) + raw_data.append(mat) + labels.append(cond) + time2 = time.time() + logger.info( + 'epoch separation done, takes %.2f s' % + (time2 - time1) + ) + 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 + + Parameters + ---------- + data_dir: str + the path to all subject files + extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, + we apply the mask right after reading a file for saving memory + epoch_file: str + the absolute path of the epoch file + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs + """ + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + labels = [] + raw_data = [] + if rank == 0: + activity_data = read_activity_data(data_dir, extension, mask_file) + # a list of numpy array in shape [condition, nEpochs, nTRs] + epoch_list = np.load(epoch_file) + raw_data, labels = separate_epochs(activity_data, epoch_list) + time1 = time.time() + raw_data_length = len(raw_data) + raw_data_length = comm.bcast(raw_data_length, root=0) + # broadcast the data subject by subject to prevent size overflow + for i in range(raw_data_length): + if rank != 0: + raw_data.append(None) + raw_data[i] = comm.bcast(raw_data[i], root=0) + if comm.Get_size() > 1: + labels = comm.bcast(labels, root=0) + if comm.Get_size() > 1 and rank == 0: + time2 = time.time() + logger.info( + 'data broadcasting done, takes %.2f s' % + (time2 - time1) + ) + return raw_data, labels diff --git a/examples/fcma/voxel_selection.py b/examples/fcma/voxel_selection.py index afcf29963..88a1e9dfc 100644 --- a/examples/fcma/voxel_selection.py +++ b/examples/fcma/voxel_selection.py @@ -13,16 +13,11 @@ # limitations under the License. from brainiak.fcma.voxelselector import VoxelSelector -from scipy.stats.mstats import zscore from sklearn import svm import sys from mpi4py import MPI -import nibabel as nib -import os -import math -import time -import numpy as np import logging +from file_io import prepare_data 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,162 +25,6 @@ logging.basicConfig(level=logging.INFO, format=format, stream=sys.stdout) logger = logging.getLogger(__name__) -def readActivityData(dir, file_extension, mask_file): - """ read data in NIfTI format and apply the spatial mask to them - - Parameters - ---------- - dir: str - the path to all subject files - file_extension: str - the file extension, usually nii.gz or nii - mask_file: str - the absolute path of the mask file, we apply the mask right after - reading a file for saving memory - - Returns - ------- - activity_data: list of 2D array in shape [nTRs, nVoxels] - the masked activity data organized in TR*voxel formats - len(activity_data) equals the number of subjects - """ - time1 = time.time() - mask_img = nib.load(mask_file) - mask = mask_img.get_data() - count = 0 - for index in np.ndindex(mask.shape): - if mask[index] != 0: - count += 1 - files = [f for f in sorted(os.listdir(dir)) - if os.path.isfile(os.path.join(dir, f)) - and f.endswith(file_extension)] - activity_data = [] - for f in files: - img = nib.load(os.path.join(dir, f)) - data = img.get_data() - (d1, d2, d3, d4) = data.shape - masked_data = np.zeros([d4, count], np.float32, order='C') - count1 = 0 - for index in np.ndindex(mask.shape): - if mask[index] != 0: - masked_data[:, count1] = np.copy(data[index]) - count1 += 1 - activity_data.append(masked_data) - logger.info( - 'file %s is loaded and masked, with data shape %s' % - (f, masked_data.shape) - ) - time2 = time.time() - logger.info( - 'data reading done, takes %.2f s' % - (time2 - time1) - ) - return activity_data - - -def separateEpochs(activity_data, epoch_list): - """ separate data into epochs of interest specified in epoch_list - and z-score them for computing correlation - - Parameters - ---------- - activity_data: list of 2D array in shape [nTRs, nVoxels] - the masked activity data organized in TR*voxel formats of all subjects - epoch_list: list of 3D 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 - - Returns - ------- - raw_data: list of 2D array in shape [epoch length, nVoxels] - the data organized in epochs - and z-scored in preparation of correlation computation - len(raw_data) equals the number of epochs - labels: list of 1D array - the condition labels of the epochs - len(labels) labels equals the number of epochs - """ - time1 = time.time() - raw_data = [] - labels = [] - for sid in range(len(epoch_list)): - epoch = epoch_list[sid] - 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 - # mat is row-major - # regardless of the order of acitvity_data[sid] - mat = activity_data[sid][sub_epoch[eid, :] == 1, :] - mat = zscore(mat, axis=0, ddof=0) - # if zscore fails (standard deviation is zero), - # set all values to be zero - mat = np.nan_to_num(mat) - mat = mat / math.sqrt(r) - raw_data.append(mat) - labels.append(cond) - time2 = time.time() - logger.info( - 'epoch separation done, takes %.2f s' % - (time2 - time1) - ) - return raw_data, labels - - -def prepareData(data_dir, extension, mask_file, epoch_file): - """ read the data in and generate epochs of interests, - then broadcast to all workers - - Parameters - ---------- - data_dir: str - the path to all subject files - extension: str - the file extension, usually nii.gz or nii - mask_file: str - the absolute path of the mask file, - we apply the mask right after reading a file for saving memory - epoch_file: str - the absolute path of the epoch file - - Returns - ------- - raw_data: list of 2D array in shape [epoch length, nVoxels] - the data organized in epochs - len(raw_data) equals the number of epochs - labels: list of 1D array - the condition labels of the epochs - len(labels) labels equals the number of epochs - """ - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - labels = [] - raw_data = [] - if rank == 0: - activity_data = readActivityData(data_dir, extension, mask_file) - # a list of numpy array in shape [condition, nEpochs, nTRs] - epoch_list = np.load(epoch_file) - raw_data, labels = separateEpochs(activity_data, epoch_list) - time1 = time.time() - raw_data_length = len(raw_data) - raw_data_length = comm.bcast(raw_data_length, root=0) - # broadcast the data subject by subject to prevent size overflow - for i in range(raw_data_length): - if rank != 0: - raw_data.append(None) - raw_data[i] = comm.bcast(raw_data[i], root=0) - labels = comm.bcast(labels, root=0) - if rank == 0: - time2 = time.time() - logger.info( - 'data broadcasting done, takes %.2f s' % - (time2 - time1) - ) - return raw_data, labels - - """ 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 @@ -201,7 +40,7 @@ def prepareData(data_dir, extension, mask_file, epoch_file): extension = sys.argv[2] mask_file = sys.argv[3] epoch_file = sys.argv[4] - raw_data, labels = prepareData(data_dir, extension, mask_file, epoch_file) + raw_data, labels = prepare_data(data_dir, extension, mask_file, epoch_file) epochs_per_subj = int(sys.argv[5]) num_subjs = int(sys.argv[6]) vs = VoxelSelector(raw_data, epochs_per_subj, labels, num_subjs) diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py new file mode 100644 index 000000000..4b3d85070 --- /dev/null +++ b/tests/fcma/test_classification.py @@ -0,0 +1,67 @@ +# 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.classifier import Classifier +from scipy.stats.mstats import zscore +from sklearn import svm +from sklearn.linear_model import LogisticRegression +import numpy as np +import math +from numpy.random import RandomState +from scipy.spatial.distance import hamming + +# specify the random state to fix the random numbers +prng = RandomState(1234567890) + +def create_epoch(idx): + row = 12 + col = 5 + mat = prng.rand(row, col).astype(np.float32) + # impose a pattern to even epochs + if idx % 2 == 0: + mat = np.sort(mat, axis=0) + mat = zscore(mat, axis=0, ddof=0) + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) + mat = mat / math.sqrt(mat.shape[0]) + return mat + +def test_classification(): + fake_raw_data = [create_epoch(i) for i in range(20)] + labels = [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1] + # 4 subjects, 4 epochs per subject + epochs_per_subj = 4 + # svm + svm_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + training_data = fake_raw_data[0: 12] + clf = Classifier(svm_clf, epochs_per_subj=epochs_per_subj) + clf.fit(training_data, labels) + y_pred = clf.predict(fake_raw_data[12:]) + expected_output = [0, 0, 0, 1, 0, 1, 0, 1] + hamming_distance = hamming(y_pred, expected_output) * len(y_pred) + assert hamming_distance <= 1, \ + 'classification via SVM does not provide correct results' + # logistic regression + lr_clf = LogisticRegression() + clf = Classifier(lr_clf, epochs_per_subj=epochs_per_subj) + clf.fit(training_data, labels[0:12]) + y_pred = clf.predict(fake_raw_data[12:]) + hamming_distance = hamming(y_pred, expected_output) * len(y_pred) + assert hamming_distance <= 1, \ + 'classification via logistic regression ' \ + 'does not provide correct results' + +if __name__ == '__main__': + test_classification() diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index cd9d1c6c6..c4d7f2a75 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -36,10 +36,7 @@ def create_epoch(): return mat def test_voxel_selection(): - fake_raw_data = [create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch()] + fake_raw_data = [create_epoch() for i in range(8)] labels = [0, 1, 0, 1, 0, 1, 0, 1] # 2 subjects, 4 epochs per subject vs = VoxelSelector(fake_raw_data, 4, labels, 2, voxel_unit=1) @@ -77,7 +74,7 @@ def test_voxel_selection(): output[tuple[0]] = int(8*tuple[1]) expected_output = [6, 3, 6, 4, 4] assert np.allclose(output, expected_output, atol=1), \ - 'voxel selection via SVM does not provide correct results' + 'voxel selection via logistic regression does not provide correct results' if __name__ == '__main__': test_voxel_selection()