diff --git a/.gitignore b/.gitignore index 892caef59..afd089b49 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,8 @@ *.swp *.lib *.o +*.npy +*.log brainiak.egg-info build @@ -18,3 +20,6 @@ htmlcov .idea .ipynb_checkpoints junit-*.xml +__pycache__ + +brainiak/fcma/cython_blas.c diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py new file mode 100644 index 000000000..5f9a9900f --- /dev/null +++ b/brainiak/fcma/__init__.py @@ -0,0 +1,14 @@ +# 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""" diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx new file mode 100644 index 000000000..6f4e2ad3f --- /dev/null +++ b/brainiak/fcma/cython_blas.pyx @@ -0,0 +1,203 @@ +# 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. + +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 + + 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 resulting matrix in shape [num_assigned_voxels, num_voxels] + is stored in an alternate way to make sure that + the correlation vectors of the same voxel stored continuously + + 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_voxels + + py_n: int + the column of the resulting matrix C + in our case, is num_assigned_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 + + py_alpha: float + the weight applied to the first matrix A + + py_a: 2D array in shape [epoch_length, num_voxels] in our case + the activity data of an epoch + + py_lda: int + the stride of the first matrix A + + py_start_voxel: int + the starting voxel of assigned voxels + used to locate the second matrix B + + py_ldb: int + the stride of the second matrix B + + py_beta: float + the weight applied to the resulting matrix C + + py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + place to store the resulting correlation values + + py_ldc: int + the stride of the resulting matrix + in our case, num_voxels*num_epochs + + py_start_epoch: int + the epoch over which the correlation is computed + + Returns + ------- + py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + write the resulting correlation values in an alternate way + for the processing epoch + """ + 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, 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, + py_beta, py_c, py_ldc): + """ use blas API wrapped by scipy.linalg.cython_blas 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 + + In SVM with linear kernel, the distance of two samples + is essentially the dot product of them. + Therefore, the kernel matrix can be obtained by matrix multiplication. + 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. + + 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_epochs + + 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_voxels + + 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_start_voxel: int + the processed voxel + used to locate the input matrix A + + py_lda: int + the stride of the input matrix A + + py_beta: float + the weight applied to the resulting matrix C + + py_c: 2D array in shape [num_epochs, num_epochs] + place to store the resulting kernel matrix + + py_ldc: int + the stride of the resulting matrix + + Returns + ------- + py_c: 2D array in shape [num_epochs, num_epochs] + write the resulting kernel_matrix + for the processing voxel + """ + 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[py_start_voxel, 0, 0], &lda, + &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 + # complete the other half of the kernel matrix + 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] + else: + for j in range(py_c.shape[0]): + for k in range(j): + py_c[k, j] = py_c[j, k] + +def installed(): + """ + This is an empty method for installing cython_blas library + """ diff --git a/brainiak/fcma/src/fcma_extension.cc b/brainiak/fcma/src/fcma_extension.cc new file mode 100644 index 000000000..260b71e78 --- /dev/null +++ b/brainiak/fcma/src/fcma_extension.cc @@ -0,0 +1,97 @@ +/* +// 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. +// */ +// +#include "Python.h" +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +void within_subject_norm_native(py::array_t py_data, int epochsPerSubj) +{ +//This function calculates the within-subject normalization of the input correlation values +// +//Input parameters are: +//py_data: the correlation value array, in shape [num_voxels, num_epochs, num_selected_voxels] +// containing the raw correlation values +//epochsPerSubj: the number pf epochs per subject +// assuming all subjects have the same number of epochs +// +//Ouput parameters are: +//py_data: the correlation value array, in shape [num_voxels, num_epochs, num_selected_voxels] +// containing the normalized correlation values +// + py::buffer_info buf = py_data.request(); + float* data = (float*)buf.ptr; + // sanity check + if (buf.ndim != 3) + throw std::runtime_error("The multi-subject correlation data structure must be 3D"); + size_t nSelectedVoxels = buf.shape[0]; + size_t nEpochs = buf.shape[1]; + size_t nVoxels = buf.shape[2]; + size_t nSubjs = nEpochs / epochsPerSubj; + + #pragma omp parallel for + for(int v = 0 ; v < nSelectedVoxels*nSubjs ; v++) + { + int s = v % nSubjs; // subject id + int i = v / nSubjs; // voxel id + + float *mat = data+i*nEpochs*nVoxels; + #pragma simd + for(int j = 0 ; j < nVoxels ; j++) + { + float mean = 0.0f; + float std_dev = 0.0f; + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + float num = 1.0f + mat[b*nVoxels+j]; + float den = 1.0f - mat[b*nVoxels+j]; + num = (num <= 0.0f) ? 1e-4 : num; + den = (den <= 0.0f) ? 1e-4 : den; + mat[b*nVoxels+j] = 0.5f * logf(num/den); + mean += mat[b*nVoxels+j]; + std_dev += mat[b*nVoxels+j] * mat[b*nVoxels+j]; + } + mean = mean / (float)epochsPerSubj; + std_dev = std_dev / (float)epochsPerSubj - mean*mean; + float inv_std_dev = (std_dev <= 0.0f) ? 0.0f : 1.0f / sqrt(std_dev); + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + mat[b*nVoxels+j] = (mat[b*nVoxels+j] - mean) * inv_std_dev; + } + } + } + return; +} + +PYBIND11_PLUGIN(fcma_extension) { + py::module m("fcma_extension"); + m.def("normalization", &within_subject_norm_native, "Within-subejct correlation normalization"); + +#ifdef VERSION_INFO + m.attr("__version__") = py::str(VERSION_INFO); +#else + m.attr("__version__") = py::str("dev"); +#endif + return m.ptr(); +} diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py new file mode 100644 index 000000000..fc0a747e6 --- /dev/null +++ b/brainiak/fcma/voxelselector.py @@ -0,0 +1,411 @@ +# 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-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. +""" + +# Authors: Yida Wang +# (Intel Labs), 2016 + +import numpy as np +import time +from mpi4py import MPI +from scipy.stats.mstats import zscore +from sklearn import cross_validation +import sklearn +from . import fcma_extension +from . import cython_blas as blas +import logging + +logger = logging.getLogger(__name__) + +__all__ = [ + "VoxelSelector", +] + + +class VoxelSelector: + """Correlation-based voxel selection component of FCMA + + Parameters + ---------- + + raw_data: list of 2D array in shape [epoch length, nVoxels] + 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 + + epochs_per_subj: int + The number of epochs of each subject + + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs + + num_folds: int + 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 + + master_rank: int, default 0 + The process which serves as the master + """ + def __init__(self, + raw_data, + epochs_per_subj, + labels, + num_folds, + voxel_unit=100, + master_rank=0): + self.raw_data = raw_data + self.epochs_per_subj = epochs_per_subj + self.num_voxels = raw_data[0].shape[1] + self.labels = labels + self.num_folds = num_folds + self.voxel_unit = voxel_unit + self.master_rank = master_rank + if self.num_voxels == 0: + raise ValueError('Zero processed voxels') + if MPI.COMM_WORLD.Get_size() == 1: + raise RuntimeError('one process cannot run the ' + 'master-worker model') + if self.master_rank >= MPI.COMM_WORLD.Get_size(): + logger.warn('Master rank exceeds the number of ' + 'launched processes, set to 0') + self.master_rank = 0 + + # tags for MPI messages + _WORKTAG = 0 + _TERMINATETAG = 1 + + def run(self, clf): + """ run correlation-based voxel selection in master-worker model + + Sort the voxels based on the cross-validation accuracy + of their correlation vectors + + 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 + """ + rank = MPI.COMM_WORLD.Get_rank() + if rank == self.master_rank: + results = self._master() + # Sort the voxels + results.sort(key=lambda tup: tup[1], reverse=True) + else: + self._worker(clf) + results = [] + return results + + def _master(self): + """ master node's operation + + Assigning tasks to workers and collecting results from them + + Parameters + ---------- + None + + 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( + 'Master starts to allocate tasks' + ) + results = [] + comm = MPI.COMM_WORLD + size = comm.Get_size() + sending_voxels = self.voxel_unit if self.voxel_unit < self.num_voxels \ + else self.num_voxels + current_task = (0, sending_voxels) + status = MPI.Status() + # using_size is used when the number of tasks + # is smaller than the number of workers + using_size = size + for i in range(0, size): + if i == self.master_rank: + continue + if current_task[1] == 0: + using_size = i + break + comm.send(current_task, + dest=i, + tag=self._WORKTAG) + next_start = current_task[0] + current_task[1] + sending_voxels = self.voxel_unit \ + if self.voxel_unit < self.num_voxels - next_start \ + else self.num_voxels - next_start + current_task = (next_start, sending_voxels) + + while using_size == size: + if current_task[1] == 0: + break + result = comm.recv(source=MPI.ANY_SOURCE, + tag=MPI.ANY_TAG, + status=status) + results += result + comm.send(current_task, + dest=status.Get_source(), + tag=self._WORKTAG) + next_start = current_task[0] + current_task[1] + sending_voxels = self.voxel_unit \ + if self.voxel_unit < self.num_voxels - next_start \ + else self.num_voxels - next_start + current_task = (next_start, sending_voxels) + + for i in range(0, using_size): + if i == self.master_rank: + continue + result = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG) + results += result + + for i in range(0, size): + if i == self.master_rank: + continue + comm.send(None, + dest=i, + tag=self._TERMINATETAG) + + return results + + def _worker(self, clf): + """ worker node's operation + + Receiving tasks from the master to process and sending the result back + + Parameters + ---------- + clf: classification function + the classifier to be used in cross validation + + Returns + ------- + None + """ + comm = MPI.COMM_WORLD + status = MPI.Status() + while 1: + task = comm.recv(source=self.master_rank, + tag=MPI.ANY_TAG, + status=status) + if status.Get_tag(): + break + comm.send(self._voxelScoring(task, clf), + dest=self.master_rank) + + def _correlationComputation(self, task): + """ use BLAS API to do correlation computation (matrix multiplication) + + Parameters + ---------- + task: tuple (start_voxel_id, num_assigned_voxels) + depicting the voxels assigned to compute + + Returns + ------- + corr: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + the correlation values of all subjects in all epochs + for the assigned values, in row-major + corr[i, e, s + j] = corr[j, e, s + i] + """ + time1 = time.time() + s = task[0] + nEpochs = len(self.raw_data) + corr = np.zeros((task[1], nEpochs, self.num_voxels), + np.float32, order='C') + count = 0 + for mat in self.raw_data: + no_trans = 'N' + trans = 'T' + blas.compute_correlation(no_trans, trans, + self.num_voxels, task[1], + mat.shape[0], 1.0, + mat, self.num_voxels, + s, self.num_voxels, + 0.0, corr, + self.num_voxels * nEpochs, count) + count += 1 + time2 = time.time() + logger.debug( + 'correlation computation for %d voxels, takes %.2f s' % + ((task[1] - task[0]), (time2 - time1)) + ) + return corr + + def _correlationNormalization(self, corr): + """ within-subject normalization + + This method uses scipy.zscore to normalize the data, + but is much slower than its C++ counterpart。 + It is doing in-place z-score. + + Parameters + ---------- + corr: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + the correlation values of all subjects in all epochs + for the assigned values, in row-major + + Returns + ------- + corr: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + the normalized correlation values of all subjects in all epochs + for the assigned values, in row-major + """ + time1 = time.time() + (sv, e, av) = corr.shape + for i in range(sv): + start = 0 + while start < e: + cur_val = corr[i, start: start + self.epochs_per_subj, :] + cur_val = .5 * np.log(cur_val + 1) / (1 - cur_val) + corr[i, start: start + self.epochs_per_subj, :] = \ + zscore(cur_val, axis=0, ddof=0) + start += self.epochs_per_subj + # if zscore fails (standard deviation is zero), + # set all values to be zero + corr = np.nan_to_num(corr) + time2 = time.time() + logger.debug( + 'within-subject normalization for %d voxels ' + 'using numpy zscore function, takes %.2f s' % + (sv, (time2 - time1)) + ) + return corr + + def _crossValidation(self, task, corr, clf): + """ voxelwise cross validation based on correlation vectors + + Parameters + ---------- + task: tuple (start_voxel_id, num_assigned_voxels) + depicting the voxels assigned to compute + corr: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + the normalized correlation values of all subjects in all epochs + for the assigned values, in row-major + + 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 assigned voxels + """ + time1 = time.time() + (sv, e, av) = corr.shape + kernel_matrix = np.zeros((e, e), np.float32, order='C') + results = [] + for i in range(sv): + if isinstance(clf, sklearn.svm.SVC) \ + and clf.kernel == 'precomputed': + blas.compute_kernel_matrix('L', 'T', + e, self.num_voxels, + 1.0, corr, + i, self.num_voxels, + 0.0, kernel_matrix, e) + data = kernel_matrix + else: + data = corr[i, :, :] + # no shuffling in cv + skf = cross_validation.StratifiedKFold(self.labels, + n_folds=self.num_folds, + shuffle=False) + scores = cross_validation.cross_val_score(clf, data, self.labels, + cv=skf, n_jobs=1) + results.append((i + task[0], scores.mean())) + time2 = time.time() + logger.debug( + 'cross validation for %d voxels, takes %.2f s' % + ((task[1] - task[0]), (time2 - time1)) + ) + return results + + def _voxelScoring(self, task, clf): + """ voxel selection processing done in the worker node + + Take the task in, + do analysis on voxels specified by the task (voxel id, num_voxels) + It is a three-stage pipeline consisting of: + 1. correlation computation + 2. within-subject normalization + 3. voxelwise cross validaion + + Parameters + ---------- + task: tuple (start_voxel_id, num_assigned_voxels), + depicting the voxels assigned to compute + 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 assigned voxels + """ + time1 = time.time() + # correlation computation + corr = self._correlationComputation(task) + # normalization + # corr = self._correlationNormalization(corr) + time3 = time.time() + fcma_extension.normalization(corr, self.epochs_per_subj) + time4 = time.time() + logger.debug( + 'within-subject normalization for %d voxels ' + 'using C++, takes %.2f s' % + ((task[1] - task[0]), (time4 - time3)) + ) + + # cross validation + results = self._crossValidation(task, corr, clf) + time2 = time.time() + logger.info( + 'in rank %d, task %d takes %.2f s' % + (MPI.COMM_WORLD.Get_rank(), + (int(task[0] / self.voxel_unit)), (time2 - time1)) + ) + return results diff --git a/examples/fcma/voxel_selection.py b/examples/fcma/voxel_selection.py new file mode 100644 index 000000000..afcf29963 --- /dev/null +++ b/examples/fcma/voxel_selection.py @@ -0,0 +1,214 @@ +# 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.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 + +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__) + +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 + data/fs_epoch_labels.npy 12 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, labels = prepareData(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) + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=10 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) + results = vs.run(clf) + # this output is just for result checking + if MPI.COMM_WORLD.Get_rank()==0: + print(results[0:100]) diff --git a/requirements-dev.txt b/requirements-dev.txt index b7cd60cf5..37097bf57 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,7 +2,6 @@ coverage flake8 flake8-print pytest -pytest-cython restructuredtext-lint sphinx sphinx_rtd_theme diff --git a/setup.cfg b/setup.cfg index eebc7c44b..e095400f4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,3 @@ -[tool:pytest] -addopts = - --doctest-cython - [coverage:run] source = brainiak branch = True diff --git a/setup.py b/setup.py index b2757fad0..a0eb1b366 100644 --- a/setup.py +++ b/setup.py @@ -1,8 +1,11 @@ +from distutils import sysconfig + from setuptools import setup, Extension, find_packages from setuptools.command.build_ext import build_ext import os import sys import setuptools +from copy import deepcopy __version__ = '0.2' @@ -28,10 +31,22 @@ class get_pybind_include(object): def __init__(self, user=False): self.user = user + # Required by Cython + def __add__(self, x): + import pybind11 + return pybind11.get_include(self.user) + x + def __str__(self): import pybind11 return pybind11.get_include(self.user) + def endswith(self, x): + import pybind11 + return pybind11.get_include(self.user).endswith(x) + + def startswith(self, x): + import pybind11 + return pybind11.get_include(self.user).startswith(x) ext_modules = [ Extension( @@ -43,6 +58,19 @@ def __str__(self): get_pybind_include(user=True) ], ), + Extension( + 'brainiak.fcma.fcma_extension', + ['brainiak/fcma/src/fcma_extension.cc'], + include_dirs=[ + # Path to pybind11 headers + get_pybind_include(), + get_pybind_include(user=True) + ], + ), + Extension( + 'brainiak.fcma.cython_blas', + ['brainiak/fcma/cython_blas.pyx'], + ), ] @@ -82,6 +110,12 @@ class BuildExt(build_ext): 'unix': ['-g0', '-fopenmp'], } + # FIXME Workaround for using the Intel compiler by setting the CC env var + # Other uses of ICC (e.g., cc binary linked to icc) are not supported + if (('CC' in os.environ and 'icc' in os.environ['CC']) + or 'icc' in sysconfig.get_config_var('CC')): + c_opts['unix'] += ['-lirc', '-lintlc'] + if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', '-ftemplate-depth-1024'] @@ -92,12 +126,13 @@ def build_extensions(self): if ct == 'unix': opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version()) - opts.append(cpp_flag(self.compiler)) - if has_flag(self.compiler, '-fvisibility=hidden'): - opts.append('-fvisibility=hidden') for ext in self.extensions: - ext.extra_compile_args = opts - ext.extra_link_args = opts + ext.extra_compile_args = deepcopy(opts) + ext.extra_link_args = deepcopy(opts) + lang = ext.language or self.compiler.detect_language(ext.sources) + if lang == 'c++': + ext.extra_compile_args.append(cpp_flag(self.compiler)) + ext.extra_link_args.append(cpp_flag(self.compiler)) build_ext.build_extensions(self) setup( diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py new file mode 100644 index 000000000..cd9d1c6c6 --- /dev/null +++ b/tests/fcma/test_voxel_selection.py @@ -0,0 +1,83 @@ +# 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.voxelselector import VoxelSelector +from scipy.stats.mstats import zscore +from sklearn import svm +from sklearn.linear_model import LogisticRegression +import numpy as np +import math +from mpi4py import MPI +from numpy.random import RandomState + +# specify the random state to fix the random numbers +prng = RandomState(1234567890) + +def create_epoch(): + row = 12 + col = 5 + mat = prng.rand(row, col).astype(np.float32) + 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_voxel_selection(): + fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch()] + 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) + # test scipy normalization + fake_corr = prng.rand(1, 4, 5).astype(np.float32) + fake_corr = vs._correlationNormalization(fake_corr) + if MPI.COMM_WORLD.Get_rank() == 0: + expected_fake_corr = [[[1.19203866, 0.18862808, -0.54350245, + -1.18334889, -0.16860008], + [-1.06594729, -1.08742261, -1.19447124, + 1.14114654, -0.67860204], + [0.7839641, 1.53981364, 0.24948341, + 0.82626557, 1.67902875], + [-0.91005552, -0.64101928, 1.48848987, + -0.78406328, -0.83182675]]] + assert np.allclose(fake_corr, expected_fake_corr), \ + 'within-subject normalization does not provide correct results' + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=1 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + results = vs.run(clf) + if MPI.COMM_WORLD.Get_rank() == 0: + output = [None] * len(results) + for tuple in results: + output[tuple[0]] = int(8*tuple[1]) + expected_output = [7, 4, 6, 4, 4] + assert np.allclose(output, expected_output, atol=1), \ + 'voxel selection via SVM does not provide correct results' + # for cross validation, use logistic regression + clf = LogisticRegression() + results = vs.run(clf) + if MPI.COMM_WORLD.Get_rank() == 0: + output = [None] * len(results) + for tuple in results: + 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' + +if __name__ == '__main__': + test_voxel_selection()