diff --git a/brainiak/fcma/classifier.py b/brainiak/fcma/classifier.py index 672e5a294..eb0ac707b 100644 --- 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 index ffc1617cb..f246928fa 100644 --- 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 index d1f2ecdb2..1429a5a06 100644 --- 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 new file mode 100644 index 000000000..76f5f885e --- /dev/null +++ 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 new file mode 100644 index 000000000..1922515ba --- /dev/null +++ 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 index e74816697..578d87520 100644 --- 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 new file mode 100644 index 000000000..19b4d05f9 --- /dev/null +++ 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()