From 550825d7a0e827bd00c09f4467c4282d8c85ab00 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 26 Sep 2016 21:37:20 -0700 Subject: [PATCH 01/17] start working on fcma classification --- brainiak/fcma/classifier.py | 125 +++++++++++++++++++++ brainiak/fcma/cython_blas.pyx | 32 +++++- examples/fcma/classification.py | 37 +++++++ examples/fcma/file_io.py | 180 +++++++++++++++++++++++++++++++ examples/fcma/voxel_selection.py | 163 +--------------------------- 5 files changed, 373 insertions(+), 164 deletions(-) create mode 100644 brainiak/fcma/classifier.py create mode 100644 examples/fcma/classification.py create mode 100644 examples/fcma/file_io.py diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py new file mode 100644 index 000000000..819db67f8 --- /dev/null +++ b/brainiak/fcma/classifier.py @@ -0,0 +1,125 @@ +# 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.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): + """ + the data has been processed by top voxels and prapared for correlation computation + """ + def __init__(self, + epochs_per_subj=0, + clf=None): + self.epochs_per_subj = epochs_per_subj + self.clf = clf + return + + def fit(self, X, y): + """ + Parameters: + ---------- + X: a list of numpy array in shape [nun_TRs, num_voxels] + 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_TRs = X[0].shape[0] + num_voxels = X[0].shape[1] + corr_data = np.zeros((num_samples, num_voxels, num_voxels), + np.float32, order='C') + # compute correlation + count = 0 + for data in X: + blas.compute_single_self_correlation('L', 'N', + num_voxels, + num_TRs, + 1.0, data, + num_voxels, 0.0, + corr_data, + num_voxels, count) + count += 1 + logger.info( + 'correlation computation done' + ) + # normalization 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.info( + '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 + logger.info( + 'kernel computation done' + ) + else: + data = corr_data + self.clf = self.clf.fit(data, y) + time2 = time.time() + logger.info( + 'training done, takes %.2f s' % + (time2 - time1) + ) diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 6f4e2ad3f..19fc81b61 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -197,7 +197,35 @@ 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(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_lda, + py_beta, py_c, py_ldc, int py_start_sample): """ - This is an empty method for installing cython_blas library + assumption: py_ldc == py_n """ + 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] diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py new file mode 100644 index 000000000..5830c959c --- /dev/null +++ b/examples/fcma/classification.py @@ -0,0 +1,37 @@ +# 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 +import sys +import logging +from file_io import prepareData + +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__) + +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 = prepareData(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) + clf = Classifier(epochs_per_subj, use_clf) + clf.fit(raw_data, labels) diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py new file mode 100644 index 000000000..600486468 --- /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 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) + 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..c2fa0fc8f 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 prepareData 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 From abcf911a5436caefaed182a2e78e1146583ce305 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 26 Sep 2016 22:06:58 -0700 Subject: [PATCH 02/17] use gemm to compute kernel matrix for training --- brainiak/fcma/classifier.py | 14 +++++++++++--- brainiak/fcma/cython_blas.pyx | 30 ++++++++++++++++++++++++++++-- examples/fcma/classification.py | 3 +++ 3 files changed, 42 insertions(+), 5 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 819db67f8..9b57ce464 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -81,12 +81,20 @@ def fit(self, X, y): # compute correlation count = 0 for data in X: - blas.compute_single_self_correlation('L', 'N', + #blas.compute_single_self_correlation('L', 'N', + # num_voxels, + # num_TRs, + # 1.0, data, + # num_voxels, 0.0, + # corr_data, + # num_voxels, count) + blas.compute_single_self_correlation2('N', 'T', + num_voxels, num_voxels, num_TRs, 1.0, data, - num_voxels, 0.0, - corr_data, + num_voxels, num_voxels, + 0.0, corr_data, num_voxels, count) count += 1 logger.info( diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 19fc81b61..9b6ef565c 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -100,7 +100,7 @@ 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, py_beta, py_c, py_ldc): @@ -221,7 +221,7 @@ def compute_single_self_correlation(py_uplo, py_trans, py_n, py_k, py_alpha, py_ 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'): + 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] @@ -229,3 +229,29 @@ def compute_single_self_correlation(py_uplo, py_trans, py_n, py_k, py_alpha, py_ 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_correlation2(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): + """ + assumption: py_ldc == py_n + """ + 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) diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index 5830c959c..1c202bc7f 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -24,6 +24,9 @@ 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] From 3ce46f3fd95165432e3aedde268d3e12442d739e Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 26 Sep 2016 23:24:00 -0700 Subject: [PATCH 03/17] prediction is in --- brainiak/fcma/classifier.py | 93 ++++++++++++++++++++++++++++++--- examples/fcma/classification.py | 7 ++- 2 files changed, 92 insertions(+), 8 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 9b57ce464..d893ef66a 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -58,6 +58,8 @@ def __init__(self, clf=None): self.epochs_per_subj = epochs_per_subj self.clf = clf + self.training_data = None + self.num_voxels = -1 return def fit(self, X, y): @@ -65,22 +67,24 @@ def fit(self, X, y): Parameters: ---------- X: a list of numpy array in shape [nun_TRs, num_voxels] + 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_TRs = X[0].shape[0] - num_voxels = X[0].shape[1] + num_voxels = X[0].shape[1] # see assumption above 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] #blas.compute_single_self_correlation('L', 'N', # num_voxels, # num_TRs, @@ -97,15 +101,15 @@ def fit(self, X, y): 0.0, corr_data, num_voxels, count) count += 1 - logger.info( + logger.debug( 'correlation computation done' ) - # normalization if necessary + # 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.info( + logger.debug( 'normalization done' ) # training @@ -120,14 +124,89 @@ def fit(self, X, y): 0, num_voxels * num_voxels, 0.0, kernel_matrix, num_samples) data = kernel_matrix - logger.info( + # 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 + data = corr_data.reshape(num_samples, num_voxels * num_voxels) + self.num_voxels = num_voxels self.clf = self.clf.fit(data, y) time2 = time.time() logger.info( 'training done, takes %.2f s' % (time2 - time1) ) + + def predict(self, X): + """ + Parameters: + ---------- + X: a list of numpy array in shape [nun_TRs, num_voxels] + len(X) equals num_samples + if num_samples > 0: normalization is done on all subjects + num_voxels equals the one used in the model + + Returns: + ------- + y_pred: the predicted label of X, in shape [num_samples,] + """ + time1 = time.time() + num_samples = len(X) + assert num_samples > 0, \ + 'at least one sample is needed' + corr_data = np.zeros((num_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_correlation2('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_samples > 1: + corr_data = corr_data.reshape(1, num_samples, num_voxels*num_voxels) + fcma_extension.normalization(corr_data, num_samples) + corr_data = corr_data.reshape(num_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[1] + #data = np.zeros((num_samples, num_training_samples), np.float32, order='C') + corr_data = corr_data.reshape(num_samples, num_voxels * num_voxels) + # compute the similarity matrix using corr_data and training_data + data = np.dot(corr_data, self.training_data.transpose()) + print(data.shape, corr_data.shape, self.training_data.shape) + logger.info( + 'similarity matrix computation done' + ) + else: + data = corr_data.reshape(num_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/examples/fcma/classification.py b/examples/fcma/classification.py index 1c202bc7f..43329a7a9 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -17,6 +17,7 @@ import sys import logging from file_io import prepareData +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, @@ -37,4 +38,8 @@ # no shrinking, set C=1 use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) clf = Classifier(epochs_per_subj, use_clf) - clf.fit(raw_data, labels) + training_data = raw_data[0:204] + test_data = raw_data[204:] + clf.fit(training_data, labels[0:204]) + print(clf.predict(test_data)) + print(np.asanyarray(labels[204:])) From f676368e5954b290f2a528118ba4cc2d7c97acbb Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 26 Sep 2016 23:54:38 -0700 Subject: [PATCH 04/17] use blas to compute similarity matrix, test --- brainiak/fcma/classifier.py | 21 +++++++++++++++++---- brainiak/fcma/cython_blas.pyx | 25 +++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 4 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index d893ef66a..5dd4cc9a3 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -192,12 +192,25 @@ def predict(self, X): 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[1] - #data = np.zeros((num_samples, num_training_samples), np.float32, order='C') + num_training_samples = self.training_data.shape[0] + data = np.zeros((num_samples, num_training_samples), np.float32, order='C') corr_data = corr_data.reshape(num_samples, num_voxels * num_voxels) # compute the similarity matrix using corr_data and training_data - data = np.dot(corr_data, self.training_data.transpose()) - print(data.shape, corr_data.shape, self.training_data.shape) + data2 = np.dot(corr_data, self.training_data.transpose()) + print('shapes:', data.shape, corr_data.shape, self.training_data.shape) + blas.compute_single_matrix_multiplication('T', 'N', + num_training_samples, + num_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) + print(self.training_data.dtype, corr_data.dtype, data.dtype, data2.dtype) + print(data[0,0], data2[0,0]) + assert np.allclose(data, data2, atol=1e-3), \ + 'error!!!' logger.info( 'similarity matrix computation done' ) diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 9b6ef565c..774ba53f0 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -255,3 +255,28 @@ def compute_single_self_correlation2(py_trans_a, py_trans_b, py_m, py_n, py_k, p 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): + 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, + &A[0, 0], &ldb, &beta, &C[0, 0], &ldc) From 7e772752bad7c0aabcb518c8d9f4a03f5c386148 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 09:44:39 -0700 Subject: [PATCH 05/17] use gemm to compute similarity matrix in prediction --- brainiak/fcma/classifier.py | 8 +------- brainiak/fcma/cython_blas.pyx | 7 ++++++- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 5dd4cc9a3..49291177f 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -196,8 +196,6 @@ def predict(self, X): data = np.zeros((num_samples, num_training_samples), np.float32, order='C') corr_data = corr_data.reshape(num_samples, num_voxels * num_voxels) # compute the similarity matrix using corr_data and training_data - data2 = np.dot(corr_data, self.training_data.transpose()) - print('shapes:', data.shape, corr_data.shape, self.training_data.shape) blas.compute_single_matrix_multiplication('T', 'N', num_training_samples, num_samples, @@ -207,11 +205,7 @@ def predict(self, X): corr_data, num_voxels * num_voxels, 0.0, data, num_training_samples) - print(self.training_data.dtype, corr_data.dtype, data.dtype, data2.dtype) - print(data[0,0], data2[0,0]) - assert np.allclose(data, data2, atol=1e-3), \ - 'error!!!' - logger.info( + logger.debug( 'similarity matrix computation done' ) else: diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 774ba53f0..14703d4e6 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -279,4 +279,9 @@ def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, py_ 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[0, 0], &ldc) + &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 From 4dac89d0589829ff13dc4f8fc1e8840a7fb3c0a1 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 18:13:21 -0700 Subject: [PATCH 06/17] add model dump and load; try logestic regression --- brainiak/fcma/classifier.py | 3 +-- examples/fcma/classification.py | 5 +++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 49291177f..1acc5c674 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -35,8 +35,6 @@ import numpy as np import time -from mpi4py import MPI -from scipy.stats.mstats import zscore from sklearn.base import BaseEstimator import sklearn from . import fcma_extension @@ -138,6 +136,7 @@ def fit(self, X, y): 'training done, takes %.2f s' % (time2 - time1) ) + return self def predict(self, X): """ diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index 43329a7a9..f68a0eb48 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -14,10 +14,12 @@ from brainiak.fcma.classifier import Classifier from sklearn import svm +from sklearn.linear_model import LogisticRegression import sys import logging from file_io import prepareData 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, @@ -37,9 +39,12 @@ 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(epochs_per_subj, use_clf) training_data = raw_data[0:204] test_data = raw_data[204:] clf.fit(training_data, labels[0:204]) + #joblib.dump(clf, 'model/logestic.pkl') + #clf = joblib.load('model/svm.pkl') print(clf.predict(test_data)) print(np.asanyarray(labels[204:])) From af6026bdd065891c82d5dae7695ee866f7600f01 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 21:57:03 -0700 Subject: [PATCH 07/17] formatting; add test code for fcma classification --- brainiak/fcma/classifier.py | 88 +++++++++++++++++------------- tests/fcma/test_classification.py | 65 ++++++++++++++++++++++ tests/fcma/test_voxel_selection.py | 2 +- 3 files changed, 117 insertions(+), 38 deletions(-) create mode 100644 tests/fcma/test_classification.py diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 1acc5c674..5099acde8 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -15,19 +15,10 @@ This implementation is based on the following publications: -.. [Wang2015-1] Full correlation matrix analysis (FCMA): An unbiased method for +.. [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. - -.. [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 @@ -47,9 +38,11 @@ "Classifier", ] + class Classifier(BaseEstimator): """ - the data has been processed by top voxels and prapared for correlation computation + the data has been processed by top voxels + and prapared for correlation computation """ def __init__(self, epochs_per_subj=0, @@ -61,14 +54,20 @@ def __init__(self, return def fit(self, X, y): - """ - Parameters: + """ use correlation data to train a model + + the input data X is activity data, which needs to be first + converted to correlation, and then normalized 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 [nun_TRs, num_voxels] assuming all elements of X has the same num_voxels value - Y: labels, len(X) equals len(Y) + y: labels, len(X) equals len(Y) - Returns: + Returns ------- self: return the object itself """ @@ -83,7 +82,7 @@ def fit(self, X, y): count = 0 for data in X: num_TRs = data.shape[0] - #blas.compute_single_self_correlation('L', 'N', + # blas.compute_single_self_correlation('L', 'N', # num_voxels, # num_TRs, # 1.0, data, @@ -91,20 +90,22 @@ def fit(self, X, y): # corr_data, # num_voxels, count) blas.compute_single_self_correlation2('N', 'T', - num_voxels, - num_voxels, - num_TRs, - 1.0, data, - num_voxels, num_voxels, - 0.0, corr_data, - num_voxels, count) + 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) + 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( @@ -113,9 +114,13 @@ def fit(self, X, y): # 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') + 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) + 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, @@ -123,7 +128,8 @@ def fit(self, X, y): 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) + self.training_data = corr_data.reshape(num_samples, + num_voxels * num_voxels) logger.debug( 'kernel computation done' ) @@ -140,14 +146,14 @@ def fit(self, X, y): def predict(self, X): """ - Parameters: + Parameters ---------- X: a list of numpy array in shape [nun_TRs, num_voxels] len(X) equals num_samples if num_samples > 0: normalization is done on all subjects num_voxels equals the one used in the model - Returns: + Returns ------- y_pred: the predicted label of X, in shape [num_samples,] """ @@ -156,15 +162,16 @@ def predict(self, X): assert num_samples > 0, \ 'at least one sample is needed' corr_data = np.zeros((num_samples, self.num_voxels, self.num_voxels), - np.float32, order='C') + 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' + 'the number of voxels provided by X does not match ' \ + 'the number of voxels defined in the model' blas.compute_single_self_correlation2('N', 'T', num_voxels, num_voxels, @@ -179,7 +186,9 @@ def predict(self, X): ) # normalize if necessary if num_samples > 1: - corr_data = corr_data.reshape(1, num_samples, num_voxels*num_voxels) + corr_data = corr_data.reshape(1, + num_samples, + num_voxels * num_voxels) fcma_extension.normalization(corr_data, num_samples) corr_data = corr_data.reshape(num_samples, num_voxels, num_voxels) logger.debug( @@ -192,7 +201,9 @@ def predict(self, X): 'when using precomputed kernel of SVM, ' \ 'all training data must be provided' num_training_samples = self.training_data.shape[0] - data = np.zeros((num_samples, num_training_samples), np.float32, order='C') + data = np.zeros((num_samples, num_training_samples), + np.float32, + order='C') corr_data = corr_data.reshape(num_samples, num_voxels * num_voxels) # compute the similarity matrix using corr_data and training_data blas.compute_single_matrix_multiplication('T', 'N', @@ -200,10 +211,13 @@ def predict(self, X): num_samples, num_voxels * num_voxels, 1.0, - self.training_data, num_voxels * num_voxels, - corr_data, num_voxels * num_voxels, + self.training_data, + num_voxels * num_voxels, + corr_data, + num_voxels * num_voxels, 0.0, - data, num_training_samples) + data, + num_training_samples) logger.debug( 'similarity matrix computation done' ) diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py new file mode 100644 index 000000000..ed6823ec2 --- /dev/null +++ b/tests/fcma/test_classification.py @@ -0,0 +1,65 @@ +# 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 + +# 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_classification(): + fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + 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, 0, 1, 0, 1] + # 3 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: 8] + clf = Classifier(epochs_per_subj, svm_clf) + clf.fit(training_data, labels[0:8]) + y_pred = clf.predict(fake_raw_data[8:]) + expected_output = [0, 0, 1, 0] + assert np.array_equal(y_pred, expected_output), \ + 'classification via SVM does not provide correct results' + # logistic regression + lr_clf = LogisticRegression() + clf = Classifier(epochs_per_subj, lr_clf) + clf.fit(training_data, labels[0:8]) + y_pred = clf.predict(fake_raw_data[8:]) + assert np.array_equal(y_pred, expected_output), \ + '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..f93cb798d 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -77,7 +77,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() From ab5449fb93cc2ccc00c04db510a4469e3e2aff77 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 22:35:11 -0700 Subject: [PATCH 08/17] use Hamming distance to measure the results and the expected results --- tests/fcma/test_classification.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py index ed6823ec2..cf6ec7be2 100644 --- a/tests/fcma/test_classification.py +++ b/tests/fcma/test_classification.py @@ -19,6 +19,7 @@ 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) @@ -51,14 +52,16 @@ def test_classification(): clf.fit(training_data, labels[0:8]) y_pred = clf.predict(fake_raw_data[8:]) expected_output = [0, 0, 1, 0] - assert np.array_equal(y_pred, expected_output), \ + 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(epochs_per_subj, lr_clf) clf.fit(training_data, labels[0:8]) y_pred = clf.predict(fake_raw_data[8:]) - assert np.array_equal(y_pred, expected_output), \ + 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__': From 05b8efa2ef53e7bf0e94ca9200b3db9cada0bd1a Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 22:42:45 -0700 Subject: [PATCH 09/17] more samples for fcma classification testing --- tests/fcma/test_classification.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py index cf6ec7be2..0d1ada27c 100644 --- a/tests/fcma/test_classification.py +++ b/tests/fcma/test_classification.py @@ -37,32 +37,38 @@ def create_epoch(): def test_classification(): fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), create_epoch(), create_epoch(), 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, 0, 1, 0, 1] - # 3 subjects, 4 epochs per subject + # 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: 8] + training_data = fake_raw_data[0: 12] clf = Classifier(epochs_per_subj, svm_clf) - clf.fit(training_data, labels[0:8]) - y_pred = clf.predict(fake_raw_data[8:]) - expected_output = [0, 0, 1, 0] + clf.fit(training_data, labels) + y_pred = clf.predict(fake_raw_data[12:]) + expected_output = [0, 1, 1, 0, 0, 0, 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(epochs_per_subj, lr_clf) - clf.fit(training_data, labels[0:8]) - y_pred = clf.predict(fake_raw_data[8:]) + clf.fit(training_data, labels[0:12]) + y_pred = clf.predict(fake_raw_data[12:]) + expected_output = [0, 1, 1, 0, 0, 1, 0, 1] hamming_distance = hamming(y_pred, expected_output) * len(y_pred) assert hamming_distance <= 1, \ - 'classification via logistic regression does not provide correct results' + 'classification via logistic regression ' \ + 'does not provide correct results' if __name__ == '__main__': test_classification() From 927df34c424da6c9d4f1b473092731f908bca848 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 27 Sep 2016 23:46:07 -0700 Subject: [PATCH 10/17] docstrings to cython_blas.pyx --- brainiak/fcma/classifier.py | 35 ++--- brainiak/fcma/cython_blas.pyx | 244 ++++++++++++++++++++++++++++++---- 2 files changed, 235 insertions(+), 44 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 5099acde8..b8ca08344 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -82,21 +82,22 @@ def fit(self, X, y): count = 0 for data in X: num_TRs = data.shape[0] - # blas.compute_single_self_correlation('L', 'N', + # 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_correlation2('N', 'T', - num_voxels, - num_voxels, - num_TRs, - 1.0, data, - num_voxels, 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' @@ -172,14 +173,14 @@ def predict(self, X): 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_correlation2('N', 'T', - num_voxels, - num_voxels, - num_TRs, - 1.0, data, - num_voxels, 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' diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 14703d4e6..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, @@ -102,12 +106,13 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a 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, +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,10 +203,66 @@ 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 compute_single_self_correlation(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_lda, - py_beta, py_c, py_ldc, int py_start_sample): - """ - assumption: py_ldc == py_n +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() @@ -221,7 +283,7 @@ def compute_single_self_correlation(py_uplo, py_trans, py_n, py_k, py_alpha, py_ 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'): + 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] @@ -230,10 +292,74 @@ def compute_single_self_correlation(py_uplo, py_trans, py_n, py_k, py_alpha, py_ for k in range(j): py_c[py_start_sample, k, j] = py_c[py_start_sample, j, k] -def compute_single_self_correlation2(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): - """ - assumption: py_ldc == py_n +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 """ cdef bytes by_trans_a=py_trans_a.encode() cdef bytes by_trans_b=py_trans_b.encode() @@ -256,8 +382,71 @@ def compute_single_self_correlation2(py_trans_a, py_trans_b, py_m, py_n, py_k, p 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): +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 @@ -280,8 +469,9 @@ def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, py_ 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 + # 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 From 932c5141c7d6c89aa1d1a55ccb52b504f2db1302 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 28 Sep 2016 15:27:21 -0700 Subject: [PATCH 11/17] add number of training samples as an element of Classifier of FCMA --- brainiak/fcma/classifier.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index b8ca08344..b87ec95d6 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -51,6 +51,7 @@ def __init__(self, self.clf = clf self.training_data = None self.num_voxels = -1 + self.num_samples = -1 return def fit(self, X, y): @@ -137,6 +138,7 @@ def fit(self, X, y): else: data = corr_data.reshape(num_samples, num_voxels * num_voxels) self.num_voxels = num_voxels + self.num_samples = num_samples self.clf = self.clf.fit(data, y) time2 = time.time() logger.info( From ee6e36cce3b5951bac891c82abeb4380d0fa2644 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 28 Sep 2016 21:17:40 -0700 Subject: [PATCH 12/17] typo --- examples/fcma/classification.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index f68a0eb48..4a0f0c2f2 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -44,7 +44,7 @@ training_data = raw_data[0:204] test_data = raw_data[204:] clf.fit(training_data, labels[0:204]) - #joblib.dump(clf, 'model/logestic.pkl') + #joblib.dump(clf, 'model/logistic.pkl') #clf = joblib.load('model/svm.pkl') print(clf.predict(test_data)) print(np.asanyarray(labels[204:])) From f843f720933bc96663f6e11975ff8901203ffbde Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 29 Sep 2016 12:42:24 -0700 Subject: [PATCH 13/17] docstrings of classifier.py --- brainiak/fcma/classifier.py | 101 +++++++++++++++++++++--------- examples/fcma/classification.py | 7 ++- tests/fcma/test_classification.py | 4 +- 3 files changed, 78 insertions(+), 34 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index b87ec95d6..f4b953149 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -40,15 +40,40 @@ class Classifier(BaseEstimator): - """ - the data has been processed by top voxels - and prapared for correlation computation + """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 + + 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, - epochs_per_subj=0, - clf=None): - self.epochs_per_subj = epochs_per_subj + clf, + epochs_per_subj=0): self.clf = clf + self.epochs_per_subj = epochs_per_subj self.training_data = None self.num_voxels = -1 self.num_samples = -1 @@ -57,15 +82,17 @@ def __init__(self, def fit(self, X, y): """ use correlation data to train a model - the input data X is activity data, which needs to be first - converted to correlation, and then normalized within subject - if more than one sample in one subject, and then fit to a model - defined by self.clf + the input data X is activity data filtered by top voxels + and prepared for correlation computation. + X needs to be first converted to correlation, + and then normalized 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 [nun_TRs, num_voxels] - assuming all elements of X has the same num_voxels value + assuming all elements of X has the same num_voxels value y: labels, len(X) equals len(Y) Returns @@ -77,6 +104,8 @@ def fit(self, X, 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 @@ -107,7 +136,7 @@ def fit(self, X, y): if self.epochs_per_subj > 0: corr_data = corr_data.reshape(1, num_samples, - num_voxels*num_voxels) + 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( @@ -137,8 +166,6 @@ def fit(self, X, y): ) else: data = corr_data.reshape(num_samples, num_voxels * num_voxels) - self.num_voxels = num_voxels - self.num_samples = num_samples self.clf = self.clf.fit(data, y) time2 = time.time() logger.info( @@ -148,23 +175,34 @@ def fit(self, X, y): return self def predict(self, X): - """ + """ use a trained model to predict correlation data + + the input data X is activity data filtered by top voxels + and prepared for correlation computation. + X needs to be first converted to correlation, + and then normalized across all samples in the list + if len(X) > 1, + and then predicted via self.clf. + Parameters ---------- X: a list of numpy array in shape [nun_TRs, num_voxels] - len(X) equals num_samples - if num_samples > 0: normalization is done on all subjects + len(X) equals num_test_samples + if num_test_samples > 0: normalization is done + on all test samples num_voxels equals the one used in the model Returns ------- - y_pred: the predicted label of X, in shape [num_samples,] + y_pred: the predicted label of X, in shape [num_test_samples,] """ time1 = time.time() - num_samples = len(X) - assert num_samples > 0, \ + num_test_samples = len(X) + assert num_test_samples > 0, \ 'at least one sample is needed' - corr_data = np.zeros((num_samples, self.num_voxels, self.num_voxels), + corr_data = np.zeros((num_test_samples, + self.num_voxels, + self.num_voxels), np.float32, order='C') # compute correlation @@ -188,12 +226,15 @@ def predict(self, X): 'correlation computation done' ) # normalize if necessary - if num_samples > 1: + if num_test_samples > 1: corr_data = corr_data.reshape(1, - num_samples, + num_test_samples, num_voxels * num_voxels) - fcma_extension.normalization(corr_data, num_samples) - corr_data = corr_data.reshape(num_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' ) @@ -204,14 +245,15 @@ def predict(self, X): 'when using precomputed kernel of SVM, ' \ 'all training data must be provided' num_training_samples = self.training_data.shape[0] - data = np.zeros((num_samples, num_training_samples), + data = np.zeros((num_test_samples, num_training_samples), np.float32, order='C') - corr_data = corr_data.reshape(num_samples, num_voxels * num_voxels) + 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_samples, + num_test_samples, num_voxels * num_voxels, 1.0, self.training_data, @@ -225,7 +267,8 @@ def predict(self, X): 'similarity matrix computation done' ) else: - data = corr_data.reshape(num_samples, num_voxels*num_voxels) + data = corr_data.reshape(num_test_samples, + num_voxels * num_voxels) y_pred = self.clf.predict(data) time2 = time.time() logger.info( diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index 4a0f0c2f2..3d19f8903 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -14,12 +14,12 @@ from brainiak.fcma.classifier import Classifier from sklearn import svm -from sklearn.linear_model import LogisticRegression +#from sklearn.linear_model import LogisticRegression import sys import logging from file_io import prepareData import numpy as np -from sklearn.externals import joblib +#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, @@ -40,10 +40,11 @@ # no shrinking, set C=1 use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) #use_clf = LogisticRegression() - clf = Classifier(epochs_per_subj, use_clf) + 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)) diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py index 0d1ada27c..c6408508c 100644 --- a/tests/fcma/test_classification.py +++ b/tests/fcma/test_classification.py @@ -52,7 +52,7 @@ def test_classification(): # svm svm_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) training_data = fake_raw_data[0: 12] - clf = Classifier(epochs_per_subj, svm_clf) + 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, 1, 1, 0, 0, 0, 0, 1] @@ -61,7 +61,7 @@ def test_classification(): 'classification via SVM does not provide correct results' # logistic regression lr_clf = LogisticRegression() - clf = Classifier(epochs_per_subj, lr_clf) + 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:]) expected_output = [0, 1, 1, 0, 0, 1, 0, 1] From 8974ca696f026f9df382f54989283b6f1966034c Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 29 Sep 2016 12:45:42 -0700 Subject: [PATCH 14/17] typo --- brainiak/fcma/voxelselector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From a2dd999c89a5682a52e7e1e80ea8e8ee0d875817 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 29 Sep 2016 14:57:11 -0700 Subject: [PATCH 15/17] impose a pattern to test_classification; address comments of PR reviews --- brainiak/fcma/classifier.py | 55 +++++++++++++++--------------- tests/fcma/test_classification.py | 19 ++++------- tests/fcma/test_voxel_selection.py | 5 +-- 3 files changed, 34 insertions(+), 45 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index f4b953149..8b53a8032 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -74,24 +74,21 @@ def __init__(self, epochs_per_subj=0): self.clf = clf self.epochs_per_subj = epochs_per_subj - self.training_data = None - self.num_voxels = -1 - self.num_samples = -1 return def fit(self, X, y): """ use correlation data to train a model - the input data X is activity data filtered by top voxels - and prepared for correlation computation. - X needs to be first converted to correlation, - and then normalized within subject + 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 [nun_TRs, num_voxels] + 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) @@ -104,8 +101,8 @@ def fit(self, X, 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 + 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 @@ -159,13 +156,15 @@ def fit(self, X, y): 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) + 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( @@ -177,32 +176,32 @@ def fit(self, X, y): def predict(self, X): """ use a trained model to predict correlation data - the input data X is activity data filtered by top voxels - and prepared for correlation computation. - X needs to be first converted to correlation, - and then normalized across all samples in the list + first compute the correlation of the input data, + and then normalize across all samples in the list if len(X) > 1, - and then predicted via self.clf. + and then predict via self.clf. Parameters ---------- - X: a list of numpy array in shape [nun_TRs, num_voxels] - len(X) equals num_test_samples - if num_test_samples > 0: normalization is done + 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. + len(X) is the number of test samples + if len(X) > 0: normalization is done on all test samples - num_voxels equals the one used in the model + num_voxels must be consistent with the one used in training Returns ------- - y_pred: the predicted label of X, in shape [num_test_samples,] + 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), + self.num_voxels_, + self.num_voxels_), np.float32, order='C') # compute correlation @@ -210,7 +209,7 @@ def predict(self, X): for data in X: num_TRs = data.shape[0] num_voxels = data.shape[1] - assert self.num_voxels == num_voxels, \ + 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', @@ -241,10 +240,10 @@ def predict(self, X): # predict if isinstance(self.clf, sklearn.svm.SVC) \ and self.clf.kernel == 'precomputed': - assert self.training_data is not None, \ + 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] + num_training_samples = self.training_data_.shape[0] data = np.zeros((num_test_samples, num_training_samples), np.float32, order='C') @@ -256,7 +255,7 @@ def predict(self, X): num_test_samples, num_voxels * num_voxels, 1.0, - self.training_data, + self.training_data_, num_voxels * num_voxels, corr_data, num_voxels * num_voxels, diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py index c6408508c..4b3d85070 100644 --- a/tests/fcma/test_classification.py +++ b/tests/fcma/test_classification.py @@ -24,10 +24,13 @@ # specify the random state to fix the random numbers prng = RandomState(1234567890) -def create_epoch(): +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 @@ -36,16 +39,7 @@ def create_epoch(): return mat def test_classification(): - fake_raw_data = [create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch(), - create_epoch(), create_epoch()] + 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 @@ -55,7 +49,7 @@ def test_classification(): 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, 1, 1, 0, 0, 0, 0, 1] + 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' @@ -64,7 +58,6 @@ def test_classification(): 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:]) - expected_output = [0, 1, 1, 0, 0, 1, 0, 1] hamming_distance = hamming(y_pred, expected_output) * len(y_pred) assert hamming_distance <= 1, \ 'classification via logistic regression ' \ diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index f93cb798d..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) From c80b884ef44fd4d41cc43326fb1d50d5a524bc71 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 29 Sep 2016 15:53:09 -0700 Subject: [PATCH 16/17] add Attributes to the class docstring --- brainiak/fcma/classifier.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 8b53a8032..1959b8014 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -54,19 +54,23 @@ class Classifier(BaseEstimator): classifier training if epochs_per_subj is specified default 0 means no within-subject normalization - training_data: 2D numpy array in shape [num_samples, num_features] + + 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 + 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 + num_samples_: int The number of samples of the training set """ def __init__(self, From d1581df0370c612cd7b081453d91690e311c0db4 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 29 Sep 2016 16:20:17 -0700 Subject: [PATCH 17/17] use PEP8 name convention; escape for in docstring --- brainiak/fcma/classifier.py | 3 +-- examples/fcma/classification.py | 4 ++-- examples/fcma/file_io.py | 10 +++++----- examples/fcma/voxel_selection.py | 4 ++-- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 1959b8014..56f79a62e 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -187,13 +187,12 @@ def predict(self, X): Parameters ---------- - X: a list of numpy array in shape [num_TRs, num_voxels] + 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 - num_voxels must be consistent with the one used in training Returns ------- diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index 3d19f8903..7e6896d53 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -17,7 +17,7 @@ #from sklearn.linear_model import LogisticRegression import sys import logging -from file_io import prepareData +from file_io import prepare_data import numpy as np #from sklearn.externals import joblib @@ -35,7 +35,7 @@ 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]) # no shrinking, set C=1 use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) diff --git a/examples/fcma/file_io.py b/examples/fcma/file_io.py index 600486468..7dfe25f00 100644 --- a/examples/fcma/file_io.py +++ b/examples/fcma/file_io.py @@ -23,7 +23,7 @@ logger = logging.getLogger(__name__) -def readActivityData(dir, file_extension, mask_file): +def read_activity_data(dir, file_extension, mask_file): """ read data in NIfTI format and apply the spatial mask to them Parameters @@ -76,7 +76,7 @@ def readActivityData(dir, file_extension, mask_file): return activity_data -def separateEpochs(activity_data, epoch_list): +def separate_epochs(activity_data, epoch_list): """ separate data into epochs of interest specified in epoch_list and z-score them for computing correlation @@ -127,7 +127,7 @@ def separateEpochs(activity_data, epoch_list): return raw_data, labels -def prepareData(data_dir, extension, mask_file, epoch_file): +def prepare_data(data_dir, extension, mask_file, epoch_file): """ read the data in and generate epochs of interests, then broadcast to all workers @@ -157,10 +157,10 @@ def prepareData(data_dir, extension, mask_file, epoch_file): labels = [] raw_data = [] if rank == 0: - activity_data = readActivityData(data_dir, extension, mask_file) + 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 = separateEpochs(activity_data, epoch_list) + 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) diff --git a/examples/fcma/voxel_selection.py b/examples/fcma/voxel_selection.py index c2fa0fc8f..88a1e9dfc 100644 --- a/examples/fcma/voxel_selection.py +++ b/examples/fcma/voxel_selection.py @@ -17,7 +17,7 @@ import sys from mpi4py import MPI import logging -from file_io import prepareData +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, @@ -40,7 +40,7 @@ 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)