From 2eced9ff1fa48862e735821cba8b6062e7e50370 Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Mon, 23 Jan 2017 13:19:45 -0500
Subject: [PATCH 1/8] start to compute SVM kernel matrices portion by portion

---
 brainiak/fcma/classifier.py          |  8 ++++----
 examples/fcma/classification.py      | 12 +++++++++++-
 examples/fcma/mvpa_classification.py | 10 +++++++++-
 3 files changed, 24 insertions(+), 6 deletions(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index eb0ac707b..8f5fa34ec 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -87,7 +87,7 @@ def __init__(self,
         self.epochs_per_subj = epochs_per_subj
         return
 
-    def _prepare_corerelation_data(self, X):
+    def _prepare_auto_corerelation_data(self, X):
         """ compute auto-correlation for the input data X
 
         Parameters
@@ -258,7 +258,7 @@ def fit(self, X, y):
         self.num_voxels_ = num_voxels
         self.num_samples_ = num_samples
         # 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,
@@ -326,7 +326,7 @@ def predict(self, X):
         time1 = time.time()
         self.test_raw_data_ = X
         # 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,
@@ -389,7 +389,7 @@ def decision_function(self, 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/examples/fcma/classification.py b/examples/fcma/classification.py
index 56c714df8..561f532de 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -19,6 +19,8 @@
 import sys
 import logging
 import numpy as np
+from scipy.spatial.distance import hamming
+import gc
 #from sklearn.externals import joblib
 
 format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
@@ -58,6 +60,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)
+        )
+        gc.collect()
diff --git a/examples/fcma/mvpa_classification.py b/examples/fcma/mvpa_classification.py
index cc2d71791..7de98fd26 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')

From 31eb8a3e5f6ff7d30b8d8463c543e51ff00dc920 Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Tue, 24 Jan 2017 11:03:20 -0500
Subject: [PATCH 2/8] partial similarity matrix computation in training

---
 brainiak/fcma/classifier.py     | 141 +++++++++++++++++++-------------
 brainiak/fcma/cython_blas.pyx   | 138 +++++++++++++++++++++++++------
 brainiak/fcma/voxelselector.py  |  15 ++--
 examples/fcma/classification.py |   2 +-
 4 files changed, 207 insertions(+), 89 deletions(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index 8f5fa34ec..6bf0f08ee 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -43,6 +43,12 @@ class Classifier(BaseEstimator):
     clf: class
         The classifier used, normally a classifier class of sklearn
 
+    num_processed_voxels: int, default 1000
+        Used for SVM with precomputed kernel,
+        every time only compute correlation between num_process_voxels and
+        the whole brain to aggregate the kernel matrices.
+        This is to better use the memory
+
     epochs_per_subj: int, default 0
         The number of epochs of each subject
         within-subject normalization will be performed during
@@ -57,7 +63,10 @@ class Classifier(BaseEstimator):
         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 +78,10 @@ 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 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
@@ -82,24 +94,34 @@ class Classifier(BaseEstimator):
     """
     def __init__(self,
                  clf,
+                 num_processed_voxels=1000,
                  epochs_per_subj=0):
         self.clf = clf
+        self.num_processed_voxels = num_processed_voxels
         self.epochs_per_subj = epochs_per_subj
         return
 
-    def _prepare_auto_corerelation_data(self, X):
+    def _prepare_auto_corerelation_data(self, X,
+                                        start_voxel=0,
+                                        num_voxels_1=0):
         """ 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_voxels_1: int, default self.num_voxels_
+            the number of voxels it computes for correlation computation
 
         Returns
         -------
         corr_data: the correlation data
-                    in shape [len(X), num_voxels, num_voxels]
+                    in shape [len(X), num_voxels_1, num_voxels]
         """
         num_samples = len(X)
         assert num_samples > 0, \
@@ -108,29 +130,19 @@ def _prepare_auto_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_voxels_1 == 0:
+            num_voxels_1 = num_voxels
+        corr_data = np.zeros((num_samples, num_voxels_1, 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_voxels_1, 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 +157,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_voxels_1, num_voxels]
         norm_unit: int
                     the number of samples on which the normalization
                     is performed
@@ -158,14 +170,14 @@ 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
             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'
             )
@@ -257,41 +269,58 @@ 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_auto_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) \
+        if not isinstance(self.clf, sklearn.svm.SVC) \
                 and self.clf.kernel == 'precomputed':
+            # 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
             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)
+            sr = 0
+            row_length = self.num_processed_voxels
+            while sr < num_voxels:
+                if row_length >= num_voxels - sr:
+                    row_length = 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,
+                    num_samples,
+                    row_length * num_voxels)
+                blas.compute_kernel_matrix('L', 'T',
+                                           num_samples, row_length * num_voxels,
+                                           1.0, corr_data,
+                                           0, row_length * num_voxels,
+                                           1.0, kernel_matrix, num_samples)
+                sr += row_length
             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)
+            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'
             )
-        else:
-            data = normalized_corr_data.reshape(num_samples,
-                                                num_voxels * num_voxels)
-            self.training_data_ = None
-
+        # training
         self.clf = self.clf.fit(data, y)
         self.test_raw_data_ = None
         self.test_data_ = None
diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx
index f246928fa..a41a52d94 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.
@@ -203,10 +205,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 +294,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 +384,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 +496,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 +522,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..eff90f5f7 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(
diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py
index 561f532de..376273dc7 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -46,7 +46,7 @@
     # 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)
+    clf = Classifier(use_clf, num_processed_voxels=2000, epochs_per_subj=num_epochs_per_subj)
 
     # doing leave-one-subject-out cross validation
     for i in range(num_subjects):

From e2000417bfe0d0a8d049f6ee451392180f669cf2 Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Tue, 24 Jan 2017 12:47:13 -0500
Subject: [PATCH 3/8] partial similarity matrix computation in prediction

---
 brainiak/fcma/classifier.py     | 99 +++++++++++++++++++++++++--------
 brainiak/fcma/cython_blas.pyx   |  6 --
 brainiak/fcma/voxelselector.py  |  6 ++
 examples/fcma/classification.py | 24 +++++++-
 4 files changed, 102 insertions(+), 33 deletions(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index 6bf0f08ee..d19e909e5 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -43,7 +43,7 @@ class Classifier(BaseEstimator):
     clf: class
         The classifier used, normally a classifier class of sklearn
 
-    num_processed_voxels: int, default 1000
+    num_processed_voxels: int, default 2000
         Used for SVM with precomputed kernel,
         every time only compute correlation between num_process_voxels and
         the whole brain to aggregate the kernel matrices.
@@ -91,19 +91,24 @@ class Classifier(BaseEstimator):
 
     num_samples_: int
         The number of samples of the training set
+
+    num_digits_: int
+        The number of digit of the first value of the kernel matrix,
+        for normalizing the similarity values accordingly
     """
     def __init__(self,
                  clf,
-                 num_processed_voxels=1000,
+                 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_auto_corerelation_data(self, X,
                                         start_voxel=0,
-                                        num_voxels_1=0):
+                                        num_voxels_1=None):
         """ compute auto-correlation for the input data X
 
         it will generate the correlation between some voxels and all voxels
@@ -115,8 +120,9 @@ def _prepare_auto_corerelation_data(self, X,
             assuming all elements of X has the same num_voxels value
         start_voxel: int, default 0
             the starting voxel id for correlation computation
-        num_voxels_1: int, default self.num_voxels_
+        num_voxels_1: int, default None
             the number of voxels it computes for correlation computation
+            if it is None, it is set to self.num_voxels
 
         Returns
         -------
@@ -130,7 +136,7 @@ def _prepare_auto_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'
-        if num_voxels_1 == 0:
+        if num_voxels_1 is None:
             num_voxels_1 = num_voxels
         corr_data = np.zeros((num_samples, num_voxels_1, num_voxels),
                              np.float32, order='C')
@@ -230,7 +236,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
@@ -242,7 +248,7 @@ def _prepare_test_data(self, corr_data):
                                      num_voxels * num_voxels)
         return data
 
-    def fit(self, X, y):
+    def fit(self, X, y, num_training_samples=None):
         """ use correlation data to train a model
 
         first compute the correlation of the input data,
@@ -257,6 +263,12 @@ 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 that 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.
+            If it is set, only those samples will be used to fit the model
 
         Returns
         -------
@@ -271,6 +283,12 @@ def fit(self, X, y):
         self.num_samples_ = num_samples
         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
@@ -282,6 +300,14 @@ def fit(self, X, y):
                                                 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')
             kernel_matrix = np.zeros((num_samples, num_samples),
                                      np.float32,
                                      order='C')
@@ -308,6 +334,13 @@ def fit(self, X, y):
                                            0, row_length * num_voxels,
                                            1.0, kernel_matrix, 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
             data = kernel_matrix
             if self.num_processed_voxels >= num_voxels:
                 # training data is in shape [num_samples, num_voxels * num_voxels]
@@ -320,10 +353,20 @@ def fit(self, X, y):
             logger.debug(
                 'kernel computation done'
             )
+
+        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)
-        self.test_raw_data_ = None
-        self.test_data_ = None
+        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' %
@@ -331,18 +374,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
@@ -353,15 +400,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_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)
+        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()
@@ -394,17 +442,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
@@ -414,7 +465,7 @@ 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
diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx
index a41a52d94..d673f892f 100644
--- a/brainiak/fcma/cython_blas.pyx
+++ b/brainiak/fcma/cython_blas.pyx
@@ -189,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]):
diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py
index eff90f5f7..6234e01c2 100644
--- a/brainiak/fcma/voxelselector.py
+++ b/brainiak/fcma/voxelselector.py
@@ -333,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 376273dc7..4b62c95ba 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -20,7 +20,6 @@
 import logging
 import numpy as np
 from scipy.spatial.distance import hamming
-import gc
 #from sklearn.externals import joblib
 
 format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
@@ -43,10 +42,30 @@
 
     raw_data, labels = prepare_fcma_data(data_dir, extension, mask_file, epoch_file)
 
+    # aggregate the similarity matrix to save memory
+    use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1)
+    clf = Classifier(use_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)
+    )
+
     # no shrinking, set C=1
     use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1)
     #use_clf = LogisticRegression()
-    clf = Classifier(use_clf, num_processed_voxels=2000, epochs_per_subj=num_epochs_per_subj)
+    clf = Classifier(use_clf, epochs_per_subj=num_epochs_per_subj)
 
     # doing leave-one-subject-out cross validation
     for i in range(num_subjects):
@@ -70,4 +89,3 @@
             (i, num_epochs_per_subj-incorrect_predict, num_epochs_per_subj,
              (num_epochs_per_subj-incorrect_predict)*1.0/num_epochs_per_subj)
         )
-        gc.collect()

From b6cc9c1967c7bc163514ebeed81e802e2d7a84be Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Tue, 24 Jan 2017 21:36:43 -0500
Subject: [PATCH 4/8] simplify the fit method and address format issues and add
 test code

---
 brainiak/fcma/classifier.py       | 102 +++++++++++++++++-------------
 tests/fcma/test_classification.py |  35 +++++++---
 2 files changed, 83 insertions(+), 54 deletions(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index d19e909e5..aa0adb931 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -103,7 +103,7 @@ def __init__(self,
         self.clf = clf
         self.num_processed_voxels = num_processed_voxels
         self.epochs_per_subj = epochs_per_subj
-        self.num_digits_= 0
+        self.num_digits_ = 0
         return
 
     def _prepare_auto_corerelation_data(self, X,
@@ -178,6 +178,7 @@ def _normalize_correlation_data(self, corr_data, norm_unit):
             num_samples = len(corr_data)
             [_, d2, d3] = corr_data.shape
             second_dimension = d2 * d3
+            # this is a shallow copy
             normalized_corr_data = corr_data.reshape(1,
                                                      num_samples,
                                                      second_dimension)
@@ -248,6 +249,45 @@ def _prepare_test_data(self, corr_data):
                                      num_voxels * num_voxels)
         return data
 
+    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
 
@@ -281,8 +321,8 @@ def fit(self, X, y, num_training_samples=None):
         num_voxels = X[0].shape[1]  # see assumption above
         self.num_voxels_ = num_voxels
         self.num_samples_ = num_samples
-        if not isinstance(self.clf, sklearn.svm.SVC) \
-                and self.clf.kernel == 'precomputed':
+        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(
@@ -299,51 +339,23 @@ def fit(self, X, y, num_training_samples=None):
             data = normalized_corr_data.reshape(num_samples,
                                                 num_voxels * num_voxels)
             self.training_data_ = None
-        else: # SVM with precomputed kernel
+        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')
+                    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 '
+                    raise ValueError('the number of training samples '
+                                     'must be smaller than '
                                      'the number of total samples')
-            kernel_matrix = np.zeros((num_samples, num_samples),
-                                     np.float32,
-                                     order='C')
-            sr = 0
-            row_length = self.num_processed_voxels
-            while sr < num_voxels:
-                if row_length >= num_voxels - sr:
-                    row_length = 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,
-                    num_samples,
-                    row_length * num_voxels)
-                blas.compute_kernel_matrix('L', 'T',
-                                           num_samples, row_length * num_voxels,
-                                           1.0, corr_data,
-                                           0, row_length * num_voxels,
-                                           1.0, kernel_matrix, 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
-            data = kernel_matrix
+            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]
+                # training data is in shape
+                # [num_samples, num_voxels * num_voxels]
                 self.training_data_ = normalized_corr_data.reshape(
                     num_samples,
                     num_voxels * num_voxels)
@@ -356,9 +368,9 @@ def fit(self, X, y, num_training_samples=None):
 
         if num_training_samples is not None:
             self.test_raw_data_ = data[num_training_samples:,
-                                  0:num_training_samples]
+                                       0:num_training_samples]
             self.test_data_ = data[num_training_samples:,
-                              0: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
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'

From c3b8b2a1c1041ea05139916e17d9c6f73dbc55fb Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Tue, 24 Jan 2017 21:40:51 -0500
Subject: [PATCH 5/8] format tweak

---
 examples/fcma/classification.py      | 2 +-
 examples/fcma/mvpa_classification.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py
index 4b62c95ba..8686dc4e0 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -87,5 +87,5 @@
         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)
+             (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj)
         )
diff --git a/examples/fcma/mvpa_classification.py b/examples/fcma/mvpa_classification.py
index 7de98fd26..4717440c1 100644
--- a/examples/fcma/mvpa_classification.py
+++ b/examples/fcma/mvpa_classification.py
@@ -69,6 +69,6 @@
         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)
+             (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj)
         )
     logger.info('MVPA training and classification done')

From cb6ed0090eca40e126708a5f5f755af5e0dc3327 Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Wed, 25 Jan 2017 08:53:22 -0500
Subject: [PATCH 6/8] improve correlation-based classification example

---
 examples/fcma/classification.py | 36 ++++++++++++++++++---------------
 1 file changed, 20 insertions(+), 16 deletions(-)

diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py
index 8686dc4e0..462664ba2 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -28,20 +28,7 @@
 logging.basicConfig(level=logging.INFO, format=format, stream=sys.stdout)
 logger = logging.getLogger(__name__)
 
-
-# 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_aggregating_sim_matrix(raw_data, labels):
     # aggregate the similarity matrix to save memory
     use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1)
     clf = Classifier(use_clf, num_processed_voxels=1000, epochs_per_subj=num_epochs_per_subj)
@@ -59,14 +46,14 @@
         '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)
+         (num_epochs_per_subj-incorrect_predict) * 1.0 / num_epochs_per_subj)
     )
 
+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)
-
     # doing leave-one-subject-out cross validation
     for i in range(num_subjects):
         leave_start = i * num_epochs_per_subj
@@ -89,3 +76,20 @@
             (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)

From d4b0feec9ac7ca438c79c6ac723c8cd87bab5c0c Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Thu, 26 Jan 2017 10:33:43 -0800
Subject: [PATCH 7/8] address review comments

---
 brainiak/fcma/classifier.py     | 43 +++++++++++++++++++++------------
 examples/fcma/classification.py | 10 ++++----
 2 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index aa0adb931..fcfe87ea3 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -37,6 +37,14 @@
 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
     ----------
 
@@ -45,9 +53,10 @@ class Classifier(BaseEstimator):
 
     num_processed_voxels: int, default 2000
         Used for SVM with precomputed kernel,
-        every time only compute correlation between num_process_voxels and
-        the whole brain to aggregate the kernel matrices.
-        This is to better use the memory
+        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
@@ -60,7 +69,6 @@ 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.
@@ -80,7 +88,8 @@ class Classifier(BaseEstimator):
         so that the test data does not need to be regenerated in the
         subsequent operations, e.g. getting decision values of the prediction.
         test_data\_ may also be set in the fit method
-        if SVM.SVC with precomputed kernel and the test samples are known.
+        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
@@ -93,7 +102,7 @@ class Classifier(BaseEstimator):
         The number of samples of the training set
 
     num_digits_: int
-        The number of digit of the first value of the kernel matrix,
+        The number of digits of the first value of the kernel matrix,
         for normalizing the similarity values accordingly
     """
     def __init__(self,
@@ -108,7 +117,7 @@ def __init__(self,
 
     def _prepare_auto_corerelation_data(self, X,
                                         start_voxel=0,
-                                        num_voxels_1=None):
+                                        num_processed_voxels=None):
         """ compute auto-correlation for the input data X
 
         it will generate the correlation between some voxels and all voxels
@@ -120,14 +129,14 @@ def _prepare_auto_corerelation_data(self, X,
             assuming all elements of X has the same num_voxels value
         start_voxel: int, default 0
             the starting voxel id for correlation computation
-        num_voxels_1: int, default None
+        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_1, num_voxels]
+                    in shape [len(X), num_processed_voxels, num_voxels]
         """
         num_samples = len(X)
         assert num_samples > 0, \
@@ -136,15 +145,15 @@ def _prepare_auto_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'
-        if num_voxels_1 is None:
-            num_voxels_1 = num_voxels
-        corr_data = np.zeros((num_samples, num_voxels_1, 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
         for idx, data in enumerate(X):
             num_TRs = data.shape[0]
             blas.compute_corr_vectors('N', 'T',
-                                      num_voxels, num_voxels_1, num_TRs,
+                                      num_voxels, num_processed_voxels, num_TRs,
                                       1.0, data, num_voxels,
                                       data, num_voxels,
                                       0.0, corr_data, num_voxels,
@@ -163,7 +172,7 @@ def _normalize_correlation_data(self, corr_data, norm_unit):
         Parameters
         ----------
         corr_data: the correlation data
-                    in shape [num_samples, num_voxels_1, num_voxels]
+                    in shape [num_samples, num_processed_voxels, num_voxels]
         norm_unit: int
                     the number of samples on which the normalization
                     is performed
@@ -304,11 +313,13 @@ def fit(self, X, y, num_training_samples=None):
             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 that used in the training,
+            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.
-            If it is set, only those samples will be used to fit the model
+            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
         -------
diff --git a/examples/fcma/classification.py b/examples/fcma/classification.py
index 462664ba2..7fb127140 100644
--- a/examples/fcma/classification.py
+++ b/examples/fcma/classification.py
@@ -30,8 +30,8 @@
 
 def example_of_aggregating_sim_matrix(raw_data, labels):
     # aggregate the similarity matrix to save memory
-    use_clf = svm.SVC(kernel='precomputed', shrinking=False, C=1)
-    clf = Classifier(use_clf, num_processed_voxels=1000, epochs_per_subj=num_epochs_per_subj)
+    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,
@@ -51,9 +51,9 @@ def example_of_aggregating_sim_matrix(raw_data, labels):
 
 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

From 06274498296019f6f77151bcc817a5f1fc34169f Mon Sep 17 00:00:00 2001
From: Yida Wang <yidawa@gmail.com>
Date: Thu, 26 Jan 2017 11:13:03 -0800
Subject: [PATCH 8/8] fix a format issue

---
 brainiak/fcma/classifier.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
index fcfe87ea3..f6a04a54c 100644
--- a/brainiak/fcma/classifier.py
+++ b/brainiak/fcma/classifier.py
@@ -153,7 +153,8 @@ def _prepare_auto_corerelation_data(self, X,
         for idx, data in enumerate(X):
             num_TRs = data.shape[0]
             blas.compute_corr_vectors('N', 'T',
-                                      num_voxels, num_processed_voxels, num_TRs,
+                                      num_voxels, num_processed_voxels,
+                                      num_TRs,
                                       1.0, data, num_voxels,
                                       data, num_voxels,
                                       0.0, corr_data, num_voxels,