mihaic · Jan 13, 2017 · Jan 12, 2017 · Jan 12, 2017 · Jan 13, 2017 · Jan 13, 2017
diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py
@@ -216,6 +216,12 @@ def _prepare_test_data(self, corr_data):
                                                       0.0,
                                                       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])))
+            if num_digits > 2:
+                proportion = 10**(2-num_digits)
+                data *= proportion
             logger.debug(
                 'similarity matrix computation done'
             )

diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx
@@ -469,9 +469,3 @@ def compute_single_matrix_multiplication(py_trans_a, py_trans_b, py_m, py_n,
     C = py_c
     blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda,
                &B[0, 0], &ldb, &beta, &C[0, 0], &ldc)
-    # shrink the values for getting more stable alpha values
-    # in SVM training iteration
-    num_digits = len(str(int(py_c[0, 0])))
-    if num_digits > 2:
-        proportion = 10**(2-num_digits)
-        py_c *= proportion
diff --git a/brainiak/fcma/mvpa_voxelselector.py b/brainiak/fcma/mvpa_voxelselector.py
@@ -58,6 +58,9 @@ class MVPAVoxelSelector:
     sl: Searchlight
         the distributed Searchlight object
 
+    Attributes
+    ----------
+
     processed\_data\_: 4D array in shape [brain 3D + epoch]
         contains the averaged and normalized brain data epoch by epoch
         it is generated from raw\_data and epoch\_info

diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py
@@ -0,0 +1,118 @@
+#  Copyright 2016 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+"""Full Correlation Matrix Analysis (FCMA)
+
+Correlation related high performance routines
+"""
+
+# Authors: Yida Wang
+# (Intel Labs), 2017
+
+import numpy as np
+from . import cython_blas as blas
+from scipy.stats.mstats import zscore
+import math
+
+
+def _normalize_for_correlation(data, axis):
+    """normalize the data before computing correlation
+
+    The data will be z-scored and divided by sqrt(n)
+    along the assigned axis
+
+    Parameters
+    ----------
+    data: 2D array
+
+    axis: int
+        specify which dimension of the data should be normalized
+
+    Returns
+    -------
+    data: 2D array
+        the normalized data
+    """
+    shape = data.shape
+    data = zscore(data, axis=axis, ddof=0)
+    # if zscore fails (standard deviation is zero),
+    # set all values to be zero
+    data = np.nan_to_num(data)
+    data = data / math.sqrt(shape[axis])
+    return data
+
+
+def compute_correlation(matrix1, matrix2):
+    """compute correlation between two sets of variables
+
+    Correlate the rows of matrix1 with the rows of matrix2.
+    If matrix1 == matrix2, it is auto-correlation computation
+    resulting in a symmetric correlation matrix.
+    The number of columns MUST agree between set1 and set2.
+    The correlation being computed here is
+    the Pearson's correlation coefficient, which can be expressed as
+
+    .. math:: corr(X, Y) = \\frac{cov(X, Y)}{\\sigma_X\\sigma_Y}
+
+    where cov(X, Y) is the covariance of variable X and Y, and
+
+    .. math:: \\sigma_X
+
+    is the standard deviation of variable X
+
+    Reducing the correlation computation to matrix multiplication
+    and using BLAS GEMM API wrapped by Scipy can speedup the numpy built-in
+    correlation computation (numpy.corrcoef) by one order of magnitude
+
+    .. math::
+        corr(X, Y)
+        &= \\frac{\\sum\\limits_{i=1}^n (x_i-\\bar{x})(y_i-\\bar{y})}{(n-1)
+        \\sqrt{\\frac{\\sum\\limits_{j=1}^n x_j^2-n\\bar{x}}{n-1}}
+        \\sqrt{\\frac{\\sum\\limits_{j=1}^{n} y_j^2-n\\bar{y}}{n-1}}}\\\\
+        &= \\sum\\limits_{i=1}^n(\\frac{(x_i-\\bar{x})}
+        {\\sqrt{\\sum\\limits_{j=1}^n x_j^2-n\\bar{x}}}
+        \\frac{(y_i-\\bar{y})}{\\sqrt{\\sum\\limits_{j=1}^n y_j^2-n\\bar{y}}})
+
+    Parameters
+    ----------
+    matrix1: 2D array in shape [r1, c]
+        MUST be continuous and row-major
+
+    matrix2: 2D array in shape [r2, c]
+        MUST be continuous and row-major
+
+    Returns
+    -------
+    corr_data: 2D array in shape [r1, r2]
+        continuous and row-major in np.float32
+    """
+    matrix1 = matrix1.astype(np.float32)
+    matrix2 = matrix2.astype(np.float32)
+    [r1, d1] = matrix1.shape
+    [r2, d2] = matrix2.shape
+    if d1 != d2:
+        raise ValueError('Dimension discrepancy')
+    # preprocess two components
+    matrix1 = _normalize_for_correlation(matrix1, 1)
+    matrix2 = _normalize_for_correlation(matrix2, 1)
+    corr_data = np.empty((r1, r2), dtype=np.float32, order='C')
+    # blas routine is column-major
+    blas.compute_single_matrix_multiplication('T', 'N',
+                                              r2, r1, d1,
+                                              1.0,
+                                              matrix2, d2,
+                                              matrix1, d1,
+                                              0.0,
+                                              corr_data,
+                                              r2)
+    return corr_data
diff --git a/examples/fcma/corr_comp.py b/examples/fcma/corr_comp.py
@@ -0,0 +1,79 @@
+#  Copyright 2016 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import sys
+import os
+import logging
+from file_io import generate_epochs_info
+import numpy as np
+from brainiak.fcma.util import compute_correlation
+import nibabel as nib
+import scipy.io
+
+format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
+# if want to output log to a file instead of outputting log to the console,
+# replace "stream=sys.stdout" with "filename='fcma.log'"
+logging.basicConfig(level=logging.INFO, format=format, stream=sys.stdout)
+logger = logging.getLogger(__name__)
+
+# python corr_comp.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]
+
+    raw_data = []
+    mask_img = nib.load(mask_file)
+    mask = mask_img.get_data()
+    mask = mask.astype(np.bool)
+    logger.info(
+        'mask size: %d' %
+        np.sum(mask)
+    )
+    files = [f for f in sorted(os.listdir(data_dir))
+             if os.path.isfile(os.path.join(data_dir, f))
+             and f.endswith(extension)]
+    for f in files:
+        img = nib.load(os.path.join(data_dir, f))
+        data = img.get_data()
+        # apply mask
+        data = data[mask, :]
+        raw_data.append(data)
+        logger.info(
+            'file %s is loaded, with data shape %s' % (f, data.shape)
+        )
+    epoch_list = np.load(epoch_file)
+    epoch_info = generate_epochs_info(epoch_list)
+
+    for idx, epoch in enumerate(epoch_info):
+        label = epoch[0]
+        sid = epoch[1]
+        start = epoch[2]
+        end = epoch[3]
+        mat = raw_data[sid][:, start:end]
+        mat = np.ascontiguousarray(mat, dtype=np.float32)
+        logger.info(
+            'start to compute correlation for subject %d epoch %d with label %d' %
+            (sid, idx, label)
+        )
+        corr = compute_correlation(mat, mat)
+        mdict = {}
+        mdict['corr'] = corr
+        filename = str(label) + '_' + str(sid) + '_' + str(idx)
+        logger.info(
+            'start to write the correlation matrix to disk as %s' %
+            filename
+        )
+        scipy.io.savemat(filename, mdict)
diff --git a/examples/fcma/mvpa_classification.py b/examples/fcma/mvpa_classification.py
@@ -88,7 +88,7 @@ def prepare_data(raw_data, epoch_info):
 
     return processed_data, labels
 
-# python mvpa_classification.py face_scene bet.nii.gz face_scene/visual_top_mask.nii.gz face_scene/fs_epoch_labels.npy 18
+# python mvpa_classification.py face_scene bet.nii.gz face_scene/visual_top_mask.nii.gz face_scene/fs_epoch_labels.npy
 if __name__ == '__main__':
     data_dir = sys.argv[1]
     extension = sys.argv[2]
@@ -118,8 +118,6 @@ def prepare_data(raw_data, epoch_info):
     epoch_list = np.load(epoch_file)
     epoch_info = generate_epochs_info(epoch_list)
 
-    num_subjs = int(sys.argv[5])
-
     # preparing data
     processed_data, labels = prepare_data(raw_data, epoch_info)
 

diff --git a/tests/fcma/test_util.py b/tests/fcma/test_util.py
@@ -0,0 +1,39 @@
+#  Copyright 2016 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import numpy as np
+from numpy.random import RandomState
+from brainiak.fcma.util import compute_correlation
+
+# specify the random state to fix the random numbers
+prng = RandomState(1234567890)
+
+def test_correlation_computation():
+    row1 = 5
+    col = 10
+    row2 = 6
+    mat1 = prng.rand(row1, col).astype(np.float32)
+    mat2 = prng.rand(row2, col).astype(np.float32)
+    corr = compute_correlation(mat1, mat1)
+    expected_corr = np.corrcoef(mat1)
+    assert np.allclose(corr, expected_corr, atol=1e-5), \
+        'high performance correlation computation does not provide correct correlation results within the same set'
+    corr = compute_correlation(mat1, mat2)
+    mat = np.concatenate((mat1, mat2), axis=0)
+    expected_corr = np.corrcoef(mat)[0:row1, row1:]
+    assert np.allclose(corr, expected_corr, atol=1e-5), \
+        'high performance correlation computation does not provide correct correlation results between two sets'
+
+if __name__ == '__main__':
+    test_correlation_computation()