diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index eb0ac707b..f6a04a54c 100644 --- a/brainiak/fcma/classifier.py +++ b/brainiak/fcma/classifier.py @@ -37,12 +37,27 @@ class Classifier(BaseEstimator): """Correlation-based classification component of FCMA + The classifier first computes correlation of the input data, + and normalizes them if needed, then uses the given classifier + to train and/or predict the correlation data. + NOTE: if the classifier is sklearn.svm.SVC with precomputed kernel, + the test data may be provided in the fit method to compute + the kernel matrix together with the training data to save the memory usage, + but the test data will NEVER be seen in the model training. + Parameters ---------- clf: class The classifier used, normally a classifier class of sklearn + num_processed_voxels: int, default 2000 + Used for SVM with precomputed kernel, + every time it only computes correlation between num_process_voxels and + the whole mask to aggregate the kernel matrices. + This is to save the memory + so as to handle correlations at a larger scale. + epochs_per_subj: int, default 0 The number of epochs of each subject within-subject normalization will be performed during @@ -54,10 +69,12 @@ class Classifier(BaseEstimator): ---------- 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 + the similarity vector for each sample to be classified. + However, if the test samples are also provided during the fit, + the similarity vectors can be precomputed too + and then training_data\ is None test_raw_data_: a list of 2D array in shape [num_TRs, num_voxels] default None @@ -69,7 +86,11 @@ class Classifier(BaseEstimator): default None test_data\_ is set after a prediction is called, so that the test data does not need to be regenerated in the - subsequent operations, e.g. getting decision values of the prediction + subsequent operations, e.g. getting decision values of the prediction. + test_data\_ may also be set in the fit method + if sklearn.svm.SVC with precomputed kernel + and the test samples are known. + NOTE: the test samples will never be used to fit the model. num_voxels_: int The number of voxels per brain used in this classifier @@ -79,27 +100,43 @@ class Classifier(BaseEstimator): num_samples_: int The number of samples of the training set + + num_digits_: int + The number of digits of the first value of the kernel matrix, + for normalizing the similarity values accordingly """ def __init__(self, clf, + num_processed_voxels=2000, epochs_per_subj=0): self.clf = clf + self.num_processed_voxels = num_processed_voxels self.epochs_per_subj = epochs_per_subj + self.num_digits_ = 0 return - def _prepare_corerelation_data(self, X): + def _prepare_auto_corerelation_data(self, X, + start_voxel=0, + num_processed_voxels=None): """ compute auto-correlation for the input data X + it will generate the correlation between some voxels and all voxels + Parameters ---------- X: a list of numpy array in shape [num_TRs, num_voxels] len(X) is the number of samples assuming all elements of X has the same num_voxels value + start_voxel: int, default 0 + the starting voxel id for correlation computation + num_processed_voxels: int, default None + the number of voxels it computes for correlation computation + if it is None, it is set to self.num_voxels Returns ------- corr_data: the correlation data - in shape [len(X), num_voxels, num_voxels] + in shape [len(X), num_processed_voxels, num_voxels] """ num_samples = len(X) assert num_samples > 0, \ @@ -108,29 +145,20 @@ def _prepare_corerelation_data(self, X): assert num_voxels == self.num_voxels_, \ 'the number of voxels provided by X does not match ' \ 'the number of voxels defined in the model' - corr_data = np.zeros((num_samples, num_voxels, num_voxels), + if num_processed_voxels is None: + num_processed_voxels = num_voxels + corr_data = np.zeros((num_samples, num_processed_voxels, num_voxels), np.float32, order='C') # compute correlation - count = 0 - for data in X: + for idx, data in enumerate(X): num_TRs = data.shape[0] - # syrk performs slower in this case - # blas.compute_single_self_correlation_syrk('L', 'N', - # num_voxels, - # num_TRs, - # 1.0, data, - # num_voxels, 0.0, - # corr_data, - # num_voxels, count) - blas.compute_single_self_correlation_gemm('N', 'T', - num_voxels, - num_voxels, - num_TRs, - 1.0, data, - num_voxels, num_voxels, - 0.0, corr_data, - num_voxels, count) - count += 1 + blas.compute_corr_vectors('N', 'T', + num_voxels, num_processed_voxels, + num_TRs, + 1.0, data, num_voxels, + data, num_voxels, + 0.0, corr_data, num_voxels, + start_voxel, idx) logger.debug( 'correlation computation done' ) @@ -145,7 +173,7 @@ def _normalize_correlation_data(self, corr_data, norm_unit): Parameters ---------- corr_data: the correlation data - in shape [num_samples, num_voxels, num_voxels] + in shape [num_samples, num_processed_voxels, num_voxels] norm_unit: int the number of samples on which the normalization is performed @@ -158,14 +186,15 @@ def _normalize_correlation_data(self, corr_data, norm_unit): # normalize if necessary if norm_unit > 1: num_samples = len(corr_data) - num_voxels = self.num_voxels_ + [_, d2, d3] = corr_data.shape + second_dimension = d2 * d3 + # this is a shallow copy normalized_corr_data = corr_data.reshape(1, num_samples, - num_voxels * num_voxels) + second_dimension) fcma_extension.normalization(normalized_corr_data, norm_unit) normalized_corr_data = normalized_corr_data.reshape(num_samples, - num_voxels, - num_voxels) + d2, d3) logger.debug( 'normalization done' ) @@ -218,7 +247,7 @@ def _prepare_test_data(self, corr_data): num_training_samples) # shrink the values for getting more stable alpha values # in SVM training iteration - num_digits = len(str(int(data[0, 0]))) + num_digits = self.num_digits_ if num_digits > 2: proportion = 10**(2-num_digits) data *= proportion @@ -230,7 +259,46 @@ def _prepare_test_data(self, corr_data): num_voxels * num_voxels) return data - def fit(self, X, y): + def _compute_similarity_matrix_in_portion(self, X): + kernel_matrix = np.zeros((self.num_samples_, self.num_samples_), + np.float32, + order='C') + sr = 0 + row_length = self.num_processed_voxels + normalized_corr_data = None + while sr < self.num_voxels_: + if row_length >= self.num_voxels_ - sr: + row_length = self.num_voxels_ - sr + # compute sub-correlation + corr_data = self._prepare_auto_corerelation_data(X, sr, + row_length) + # normalization + normalized_corr_data = self._normalize_correlation_data( + corr_data, + self.epochs_per_subj) + # compute partial similarity matrices + # for using kernel matrix computation from voxel selection + normalized_corr_data = normalized_corr_data.reshape( + 1, + self.num_samples_, + row_length * self.num_voxels_) + blas.compute_kernel_matrix('L', 'T', + self.num_samples_, + row_length * self.num_voxels_, + 1.0, normalized_corr_data, + 0, row_length * self.num_voxels_, + 1.0, kernel_matrix, self.num_samples_) + sr += row_length + # shrink the values for getting more stable alpha values + # in SVM training iteration + num_digits = len(str(int(kernel_matrix[0, 0]))) + self.num_digits_ = num_digits + if num_digits > 2: + proportion = 10**(2-num_digits) + kernel_matrix *= proportion + return kernel_matrix, normalized_corr_data + + def fit(self, X, y, num_training_samples=None): """ use correlation data to train a model first compute the correlation of the input data, @@ -245,6 +313,14 @@ def fit(self, X, y): and prepared for correlation computation. assuming all elements of X has the same num_voxels value y: labels, len(X) equals len(Y) + num_training_samples: int, default None + The number of samples used in the training, + which is set when the similarity matrix is constructed + portion by portion so the similarity vectors of the + test data have to be computed here. + This is ONLY set when sklearn.svm.SVC with + precomputed kernel is used. + If it is set, only those samples will be used to fit the model. Returns ------- @@ -257,44 +333,64 @@ def fit(self, X, y): num_voxels = X[0].shape[1] # see assumption above self.num_voxels_ = num_voxels self.num_samples_ = num_samples - # correlation computation - corr_data = self._prepare_corerelation_data(X) - # normalization - normalized_corr_data = self._normalize_correlation_data( - corr_data, - self.epochs_per_subj) - # 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 - normalized_corr_data = normalized_corr_data.reshape( - 1, - num_samples, - num_voxels * num_voxels) - blas.compute_kernel_matrix('L', 'T', - num_samples, num_voxels * num_voxels, - 1.0, corr_data, - 0, num_voxels * num_voxels, - 0.0, kernel_matrix, num_samples) - data = kernel_matrix - # training data is in shape [num_samples, num_voxels * num_voxels] - self.training_data_ = normalized_corr_data.reshape( - num_samples, - num_voxels * num_voxels) - logger.debug( - 'kernel computation done' - ) - else: + if not (isinstance(self.clf, sklearn.svm.SVC) + and self.clf.kernel == 'precomputed'): + if num_training_samples is not None: + num_training_samples = None + logger.warn( + 'num_training_samples should not be set for classifiers ' + 'other than SVM with precomputed kernels' + ) + # correlation computation + corr_data = self._prepare_auto_corerelation_data(X) + # normalization + normalized_corr_data = self._normalize_correlation_data( + corr_data, + self.epochs_per_subj) + # training data prepare data = normalized_corr_data.reshape(num_samples, num_voxels * num_voxels) self.training_data_ = None + else: # SVM with precomputed kernel + if self.num_processed_voxels < num_voxels: + if num_training_samples is None: + raise RuntimeError('the similarity matrix will be ' + 'computed portion by portion, ' + 'the test samples must be predefined ' + 'by specifying ' + 'num_training_samples') + if num_training_samples >= num_samples: + raise ValueError('the number of training samples ' + 'must be smaller than ' + 'the number of total samples') + data, normalized_corr_data = \ + self._compute_similarity_matrix_in_portion(X) + if self.num_processed_voxels >= num_voxels: + # training data is in shape + # [num_samples, num_voxels * num_voxels] + self.training_data_ = normalized_corr_data.reshape( + num_samples, + num_voxels * num_voxels) + else: + # do not store training data because it was partially computed + self.training_data_ = None + logger.debug( + 'kernel computation done' + ) - self.clf = self.clf.fit(data, y) - self.test_raw_data_ = None - self.test_data_ = None + if num_training_samples is not None: + self.test_raw_data_ = data[num_training_samples:, + 0:num_training_samples] + self.test_data_ = data[num_training_samples:, + 0:num_training_samples] + # limit training to the data specified by num_training_samples + data = data[0:num_training_samples, 0:num_training_samples] + # training + self.clf = self.clf.fit(data, y[0:num_training_samples]) + # set the test data + if num_training_samples is None: + self.test_raw_data_ = None + self.test_data_ = None time2 = time.time() logger.info( 'training done, takes %.2f s' % @@ -302,18 +398,22 @@ def fit(self, X, y): ) return self - def predict(self, X): + def predict(self, X=None): """ use a trained model to predict correlation data first compute the correlation of the input data, and then normalize across all samples in the list if len(X) > 1, and then predict via self.clf. + If X is None, use the similarity vectors produced in fit + to predict Parameters ---------- X: a list of numpy array in shape [num_TRs, self.num_voxels\_] - X contains the activity data filtered by top voxels + default None, meaning that the data to be predicted + have been processed in the fit method. + Otherwise, 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 @@ -324,15 +424,16 @@ def predict(self, X): y_pred: the predicted label of X, in shape [len(X),] """ time1 = time.time() - self.test_raw_data_ = X - # correlation computation - corr_data = self._prepare_corerelation_data(X) - # normalization - normalized_corr_data = self._normalize_correlation_data( - corr_data, - len(X)) - # test data generation - self.test_data_ = self._prepare_test_data(normalized_corr_data) + if X is not None: + self.test_raw_data_ = X + # correlation computation + corr_data = self._prepare_auto_corerelation_data(X) + # normalization + normalized_corr_data = self._normalize_correlation_data( + corr_data, + len(X)) + # test data generation + self.test_data_ = self._prepare_test_data(normalized_corr_data) # prediction y_pred = self.clf.predict(self.test_data_) time2 = time.time() @@ -365,17 +466,20 @@ def _is_equal_to_test_raw_data(self, X): return False return True - def decision_function(self, X): + def decision_function(self, X=None): """ output the decision value of the prediction if X is not equal to self.test_raw_data\_, i.e. predict is not called, first generate the test_data after getting the test_data, get the decision value via self.clf. + if X is None, test_data\_ is ready to be used Parameters ---------- X: a list of numpy array in shape [num_TRs, self.num_voxels\_] - X contains the activity data filtered by top voxels + default None, meaning that the data to be predicted + have been processed in the fit method. + Otherwise, 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) > 1: normalization is done @@ -385,11 +489,11 @@ def decision_function(self, X): ------- confidence: the predictions confidence values of X, in shape [len(X),] """ - if not self._is_equal_to_test_raw_data(X): + if X is not None and not self._is_equal_to_test_raw_data(X): self.test_raw_data_ = X # generate the test_data first # correlation computation - corr_data = self._prepare_corerelation_data(X) + corr_data = self._prepare_auto_corerelation_data(X) # normalization normalized_corr_data = self._normalize_correlation_data( corr_data, diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index f246928fa..d673f892f 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -17,11 +17,13 @@ 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): +def compute_self_corr_for_voxel_sel(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 + This method is limited to process self-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. @@ -187,12 +189,6 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, 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]): @@ -203,10 +199,10 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, for k in range(j): py_c[k, j] = py_c[j, k] -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): +def compute_single_self_corr_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 @@ -292,10 +288,10 @@ def compute_single_self_correlation_syrk(py_uplo, py_trans, py_n, py_k, for k in range(j): py_c[py_start_sample, k, j] = py_c[py_start_sample, j, k] -def compute_single_self_correlation_gemm(py_trans_a, py_trans_b, py_m, py_n, - py_k, py_alpha, py_a, py_lda, - py_ldb, py_beta, py_c, py_ldc, - int py_start_sample): +def compute_single_self_corr_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 @@ -382,13 +378,105 @@ def compute_single_self_correlation_gemm(py_trans_a, py_trans_b, py_m, py_n, 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_corr_vectors(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, + int py_start_voxel, + int py_start_sample): + """ use blas API gemm wrapped by scipy to construct a correlation vector + + The correlation vector is essentially correlation matrices computed + from two activity matrices. It will be placed in the corresponding place + of the resulting correlation data set. + 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 + + 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 + + py_n: int + the column of the resulting matrix C + + 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 input matrix A + + py_a: 2D array + + py_lda: int + the stride of the input matrix A + + py_b: 2D array + + 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 [py_m, py_n] of column-major + in fact it is + in shape [py_n, py_m] of row-major + + py_ldc: int + the stride of the resulting matrix + + py_start_voxel: int + the starting voxel of assigned voxels + used to locate the second matrix B + + py_start_sample: int + the processed sample + used to locate the resulting matrix C + + Returns + ------- + py_c: 2D array + in shape [py_m, py_n] of column-major + write the resulting matrix to the place indicated by py_start_sample + """ + cdef bytes by_trans_a=py_trans_a.encode() + cdef bytes by_trans_b=py_trans_b.encode() + cdef char* trans_a = by_trans_a + cdef char* trans_b = by_trans_b + cdef int M, N, K, lda, ldb, ldc + M = py_m + N = py_n + K = py_k + lda = py_lda + ldb = py_ldb + ldc = py_ldc + cdef float alpha, beta + alpha = py_alpha + beta = py_beta + cdef float[:, ::1] A + A = py_a + cdef float[:, ::1] B + B = py_b + cdef float[:, :, ::1] C + C = py_c + blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda, + &B[0, py_start_voxel], &ldb, &beta, &C[py_start_sample, 0, 0], &ldc) + def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a, py_lda, py_b, py_ldb, py_beta, py_c, py_ldc): - """ use blas API gemm wrapped by scipy to compute similarity matrix of SVM + """ use blas API gemm wrapped by scipy to do matrix multiplication - This is to compute the similarity matrix between test sample(s) and - training samples. The blas APIs process matrices in column-major, + This is to compute the matrix multiplication. + 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 @@ -402,29 +490,24 @@ def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, 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 @@ -433,19 +516,18 @@ def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n, the weight applied to the resulting matrix C py_c: 2D array - in shape [num_training_samples, num_test_samples] of column-major + in shape [py_m, py_n] 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 + in shape [py_n, py_m] of row-major 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 + py_c: 2D array + in shape [py_m, py_n] of column-major + write the resulting matrix """ cdef bytes by_trans_a=py_trans_a.encode() cdef bytes by_trans_b=py_trans_b.encode() diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index e3891e78e..6234e01c2 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -248,13 +248,14 @@ def _correlationComputation(self, task): 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) + blas.compute_self_corr_for_voxel_sel(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( @@ -332,6 +333,12 @@ def _crossValidation(self, task, corr, clf): 1.0, corr, i, self.num_voxels, 0.0, kernel_matrix, e) + # shrink the values for getting more stable alpha values + # in SVM training iteration + num_digits = len(str(int(kernel_matrix[0, 0]))) + if num_digits > 2: + proportion = 10**(2-num_digits) + kernel_matrix *= proportion data = kernel_matrix else: data = corr[i, :, :] diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py index 56c714df8..7fb127140 100644 --- a/examples/fcma/classification.py +++ b/examples/fcma/classification.py @@ -19,6 +19,7 @@ import sys import logging import numpy as np +from scipy.spatial.distance import hamming #from sklearn.externals import joblib format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s' @@ -27,25 +28,32 @@ logging.basicConfig(level=logging.INFO, format=format, stream=sys.stdout) logger = logging.getLogger(__name__) +def example_of_aggregating_sim_matrix(raw_data, labels): + # aggregate the similarity matrix to save memory + svm_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + clf = Classifier(svm_clf, num_processed_voxels=1000, epochs_per_subj=num_epochs_per_subj) + rearranged_data = raw_data[num_epochs_per_subj:] + raw_data[0:num_epochs_per_subj] + rearranged_labels = labels[num_epochs_per_subj:] + labels[0:num_epochs_per_subj] + clf.fit(rearranged_data, rearranged_labels, + num_training_samples=num_epochs_per_subj*(num_subjects-1)) + predict = clf.predict() + print(predict) + print(clf.decision_function()) + test_labels = labels[0:num_epochs_per_subj] + print(np.asanyarray(test_labels)) + incorrect_predict = hamming(predict, np.asanyarray(test_labels)) * num_epochs_per_subj + logger.info( + 'when aggregating the similarity matrix to save memory, ' + 'the accuracy is %d / %d = %.2f' % + (num_epochs_per_subj-incorrect_predict, num_epochs_per_subj, + (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj) + ) -# python classification.py face_scene bet.nii.gz face_scene/prefrontal_top_mask.nii.gz face_scene/fs_epoch_labels.npy -if __name__ == '__main__': - data_dir = sys.argv[1] - extension = sys.argv[2] - mask_file = sys.argv[3] - epoch_file = sys.argv[4] - - epoch_list = np.load(epoch_file) - num_subjects = len(epoch_list) - num_epochs_per_subj = epoch_list[0].shape[1] - - raw_data, labels = prepare_fcma_data(data_dir, extension, mask_file, epoch_file) - +def example_of_cross_validation_on_regular_classifier(raw_data, labels): # no shrinking, set C=1 - use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) - #use_clf = LogisticRegression() - clf = Classifier(use_clf, epochs_per_subj=num_epochs_per_subj) - + svm_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + #logit_clf = LogisticRegression() + clf = Classifier(svm_clf, epochs_per_subj=num_epochs_per_subj) # doing leave-one-subject-out cross validation for i in range(num_subjects): leave_start = i * num_epochs_per_subj @@ -58,6 +66,30 @@ # 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)) + predict = clf.predict(test_data) + print(predict) print(clf.decision_function(test_data)) print(np.asanyarray(test_labels)) + incorrect_predict = hamming(predict, np.asanyarray(test_labels)) * num_epochs_per_subj + logger.info( + 'when leaving subject %d out for testing, the accuracy is %d / %d = %.2f' % + (i, num_epochs_per_subj-incorrect_predict, num_epochs_per_subj, + (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj) + ) + +# python classification.py face_scene bet.nii.gz face_scene/prefrontal_top_mask.nii.gz face_scene/fs_epoch_labels.npy +if __name__ == '__main__': + data_dir = sys.argv[1] + extension = sys.argv[2] + mask_file = sys.argv[3] + epoch_file = sys.argv[4] + + epoch_list = np.load(epoch_file) + num_subjects = len(epoch_list) + num_epochs_per_subj = epoch_list[0].shape[1] + + raw_data, labels = prepare_fcma_data(data_dir, extension, mask_file, epoch_file) + + example_of_aggregating_sim_matrix(raw_data, labels) + + example_of_cross_validation_on_regular_classifier(raw_data, labels) diff --git a/examples/fcma/mvpa_classification.py b/examples/fcma/mvpa_classification.py index cc2d71791..4717440c1 100644 --- a/examples/fcma/mvpa_classification.py +++ b/examples/fcma/mvpa_classification.py @@ -18,6 +18,7 @@ import logging import brainiak.fcma.io as io import numpy as np +from scipy.spatial.distance import hamming #from sklearn.externals import joblib format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s' @@ -60,7 +61,14 @@ # 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)) + predict = clf.predict(test_data) + print(predict) print(clf.decision_function(test_data)) print(np.asanyarray(test_labels)) + incorrect_predict = hamming(predict, np.asanyarray(test_labels)) * num_epochs_per_subj + logger.info( + 'when leaving subject %d out for testing, the accuracy is %d / %d = %.2f' % + (i, num_epochs_per_subj-incorrect_predict, num_epochs_per_subj, + (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj) + ) logger.info('MVPA training and classification done') diff --git a/tests/fcma/test_classification.py b/tests/fcma/test_classification.py index 6cde5a6e8..60b56626f 100644 --- a/tests/fcma/test_classification.py +++ b/tests/fcma/test_classification.py @@ -40,20 +40,21 @@ def create_epoch(idx): def test_classification(): fake_raw_data = [create_epoch(i) for i in range(20)] - labels = [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1] - # 4 subjects, 4 epochs per subject + labels = [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1] + # 5 subjects, 4 epochs per subject epochs_per_subj = 4 # svm svm_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) training_data = fake_raw_data[0: 12] clf = Classifier(svm_clf, epochs_per_subj=epochs_per_subj) - clf.fit(training_data, labels) + clf.fit(training_data, labels[0:12]) expected_confidence = np.array([-1.18234421, 0.97403604, -1.04005679, 0.92403019, -0.95567738, 1.11746593, -0.83275891, 0.9486868]) recomputed_confidence = clf.decision_function(fake_raw_data[12:]) - hamming_distance = hamming(np.sign(expected_confidence), - np.sign(recomputed_confidence)) + hamming_distance = hamming(np.sign(expected_confidence), + np.sign(recomputed_confidence)) \ + * expected_confidence.size assert hamming_distance <= 1, \ 'decision function of SVM with recomputation ' \ 'does not provide correct results' @@ -64,10 +65,26 @@ def test_classification(): 'classification via SVM does not provide correct results' confidence = clf.decision_function(fake_raw_data[12:]) hamming_distance = hamming(np.sign(expected_confidence), - np.sign(confidence)) + np.sign(confidence)) * confidence.size assert hamming_distance <= 1, \ 'decision function of SVM without recomputation ' \ 'does not provide correct results' + # svm with partial similarity matrix computation + clf = Classifier(svm_clf, num_processed_voxels=2, + epochs_per_subj=epochs_per_subj) + clf.fit(fake_raw_data, labels, num_training_samples=12) + y_pred = clf.predict() + 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 (partial sim) does not ' \ + 'provide correct results' + confidence = clf.decision_function() + hamming_distance = hamming(np.sign(expected_confidence), + np.sign(confidence)) * confidence.size + assert hamming_distance <= 1, \ + 'decision function of SVM (partial sim) without recomputation ' \ + 'does not provide correct results' # logistic regression lr_clf = LogisticRegression() clf = Classifier(lr_clf, epochs_per_subj=epochs_per_subj) @@ -77,7 +94,7 @@ def test_classification(): -3.35616616, 3.77716609]) recomputed_confidence = clf.decision_function(fake_raw_data[12:]) hamming_distance = hamming(np.sign(expected_confidence), - np.sign(recomputed_confidence)) + np.sign(recomputed_confidence)) * expected_confidence.size assert hamming_distance <= 1, \ 'decision function of logistic regression with recomputation ' \ 'does not provide correct results' @@ -88,8 +105,8 @@ def test_classification(): 'classification via logistic regression ' \ 'does not provide correct results' confidence = clf.decision_function(fake_raw_data[12:]) - hamming_distance = hamming(np.sign(expected_confidence), - np.sign(confidence)) + hamming_distance = hamming(np.sign(expected_confidence), + np.sign(confidence)) * confidence.size assert hamming_distance <= 1, \ 'decision function of logistic regression without precomputation ' \ 'does not provide correct results'