From c62f9214160ab04f0fe4cb52255d131c85957943 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 15:28:56 -0800 Subject: [PATCH 01/11] add pure high performance correlation computation function --- brainiak/fcma/classifier.py | 6 ++++ brainiak/fcma/cython_blas.pyx | 6 ---- brainiak/fcma/util.py | 53 +++++++++++++++++++++++++++++++++++ tests/fcma/test_util.py | 32 +++++++++++++++++++++ 4 files changed, 91 insertions(+), 6 deletions(-) create mode 100644 brainiak/fcma/util.py create mode 100644 tests/fcma/test_util.py 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/util.py b/brainiak/fcma/util.py new file mode 100644 index 000000000..6020e4069 --- /dev/null +++ b/brainiak/fcma/util.py @@ -0,0 +1,53 @@ +# 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): + 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(component1, component2): + [r1, d1] = component1.shape + [r2, d2] = component2.shape + if d1 != d2: + raise ValueError('Dimension discrepancy') + # preprocess two components + component1 = normalize_for_correlation(component1, 1) + component2 = normalize_for_correlation(component2, 1) + corr_data = np.empty((r1, r2), dtype=np.float32, order='C') + blas.compute_single_matrix_multiplication('T', 'N', + r1, r2, d1, + 1.0, + component1, d1, + component2, d2, + 0.0, + corr_data, + r2) + return corr_data diff --git a/tests/fcma/test_util.py b/tests/fcma/test_util.py new file mode 100644 index 000000000..0b1c15c5a --- /dev/null +++ b/tests/fcma/test_util.py @@ -0,0 +1,32 @@ +# 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(): + row = 5 + col = 5 + mat = prng.rand(row, col).astype(np.float32) + corr = compute_correlation(mat, mat) + expected_corr = np.corrcoef(mat) + assert np.allclose(corr, expected_corr, atol=1e-5), \ + 'high performance correlation computation does not provide correct results' + +if __name__ == '__main__': + test_correlation_computation() From da16ad0be4a7eeca6279295c0dc589cf3f95d00c Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 16:14:57 -0800 Subject: [PATCH 02/11] add test case for correlation between two sets --- brainiak/fcma/util.py | 3 ++- tests/fcma/test_util.py | 15 +++++++++++---- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 6020e4069..79ae30c2f 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -42,11 +42,12 @@ def compute_correlation(component1, component2): component1 = normalize_for_correlation(component1, 1) component2 = normalize_for_correlation(component2, 1) corr_data = np.empty((r1, r2), dtype=np.float32, order='C') + # blas routine is column-major blas.compute_single_matrix_multiplication('T', 'N', r1, r2, d1, 1.0, - component1, d1, component2, d2, + component1, d1, 0.0, corr_data, r2) diff --git a/tests/fcma/test_util.py b/tests/fcma/test_util.py index 0b1c15c5a..85562c7e5 100644 --- a/tests/fcma/test_util.py +++ b/tests/fcma/test_util.py @@ -15,6 +15,7 @@ import numpy as np from numpy.random import RandomState from brainiak.fcma.util import compute_correlation +from scipy.signal import correlate # specify the random state to fix the random numbers prng = RandomState(1234567890) @@ -22,11 +23,17 @@ def test_correlation_computation(): row = 5 col = 5 - mat = prng.rand(row, col).astype(np.float32) - corr = compute_correlation(mat, mat) - expected_corr = np.corrcoef(mat) + mat1 = prng.rand(row, col).astype(np.float32) + mat2 = prng.rand(row, 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 results' + '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:row, col:] + 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() From 43a8222ae5cf54d1d9b928adf5ae6d246bcd42b8 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 16:23:24 -0800 Subject: [PATCH 03/11] fix a bug when two components have different number of items --- brainiak/fcma/util.py | 2 +- tests/fcma/test_util.py | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 79ae30c2f..7f59bef31 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -44,7 +44,7 @@ def compute_correlation(component1, component2): corr_data = np.empty((r1, r2), dtype=np.float32, order='C') # blas routine is column-major blas.compute_single_matrix_multiplication('T', 'N', - r1, r2, d1, + r2, r1, d1, 1.0, component2, d2, component1, d1, diff --git a/tests/fcma/test_util.py b/tests/fcma/test_util.py index 85562c7e5..843acf06d 100644 --- a/tests/fcma/test_util.py +++ b/tests/fcma/test_util.py @@ -21,17 +21,18 @@ prng = RandomState(1234567890) def test_correlation_computation(): - row = 5 - col = 5 - mat1 = prng.rand(row, col).astype(np.float32) - mat2 = prng.rand(row, col).astype(np.float32) + 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:row, col:] + 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' From f93fd28fd39fe5771c884063dadc81321364640a Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 16:43:33 -0800 Subject: [PATCH 04/11] add an example of using correlation computation and writing out mat files --- examples/fcma/corr_comp.py | 79 ++++++++++++++++++++++++++++ examples/fcma/mvpa_classification.py | 4 +- 2 files changed, 80 insertions(+), 3 deletions(-) create mode 100644 examples/fcma/corr_comp.py 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) From 1513e2c0993969efccd3ba073dee9c4cb6acb2ac Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 16:56:56 -0800 Subject: [PATCH 05/11] add docstrings --- brainiak/fcma/util.py | 52 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 7f59bef31..ebf3ae660 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -25,6 +25,23 @@ 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), @@ -33,21 +50,42 @@ def normalize_for_correlation(data, axis): data = data / math.sqrt(shape[axis]) return data -def compute_correlation(component1, component2): - [r1, d1] = component1.shape - [r2, d2] = component2.shape +def compute_correlation(set1, set2): + """compute correlation between two sets + + Correlate the rows of set1 with the rows of set2. + If set1 == set2, it is auto-correlation computation + resulting in a symmetric correlation matrix. + The number of columns MUST agree between set1 and set2 + + Parameters + ---------- + set1: 2D array in shape [r1, c] + MUST be continuous and row-major + + set2: 2D array in shape [r2, c] + MUST be continuous and row-major + + Returns + ------- + corr_data: 2D array in shape [r1, r2] + """ + set1 = set1.astype(np.float32) + set2 = set2.astype(np.float32) + [r1, d1] = set1.shape + [r2, d2] = set2.shape if d1 != d2: raise ValueError('Dimension discrepancy') # preprocess two components - component1 = normalize_for_correlation(component1, 1) - component2 = normalize_for_correlation(component2, 1) + set1 = normalize_for_correlation(set1, 1) + set2 = normalize_for_correlation(set2, 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, - component2, d2, - component1, d1, + set2, d2, + set1, d1, 0.0, corr_data, r2) From 21e8551c896aa6ce15b1f15dea20acb460814a00 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 17:00:36 -0800 Subject: [PATCH 06/11] formatting issues --- brainiak/fcma/util.py | 2 ++ tests/fcma/test_util.py | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index ebf3ae660..5e2d88935 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -24,6 +24,7 @@ from scipy.stats.mstats import zscore import math + def normalize_for_correlation(data, axis): """normalize the data before computing correlation @@ -50,6 +51,7 @@ def normalize_for_correlation(data, axis): data = data / math.sqrt(shape[axis]) return data + def compute_correlation(set1, set2): """compute correlation between two sets diff --git a/tests/fcma/test_util.py b/tests/fcma/test_util.py index 843acf06d..19b4d05f9 100644 --- a/tests/fcma/test_util.py +++ b/tests/fcma/test_util.py @@ -15,7 +15,6 @@ import numpy as np from numpy.random import RandomState from brainiak.fcma.util import compute_correlation -from scipy.signal import correlate # specify the random state to fix the random numbers prng = RandomState(1234567890) From b40d001e401e2ba1c0e951423fd6faaed690382c Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 17:07:41 -0800 Subject: [PATCH 07/11] some more docstrings --- brainiak/fcma/util.py | 1 + 1 file changed, 1 insertion(+) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 5e2d88935..6cbfed506 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -71,6 +71,7 @@ def compute_correlation(set1, set2): Returns ------- corr_data: 2D array in shape [r1, r2] + continuous and row-major in np.float32 """ set1 = set1.astype(np.float32) set2 = set2.astype(np.float32) From 58279563fc1add44c35211e86bb3ac5ab1e87a33 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 19:45:34 -0800 Subject: [PATCH 08/11] addressing review comments --- brainiak/fcma/util.py | 40 +++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 6cbfed506..12a40f7b9 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -52,20 +52,30 @@ def normalize_for_correlation(data, axis): return data -def compute_correlation(set1, set2): - """compute correlation between two sets +def compute_correlation(matrix1, matrix2): + """compute correlation between two sets of variables - Correlate the rows of set1 with the rows of set2. - If set1 == set2, it is auto-correlation computation + 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 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 Parameters ---------- - set1: 2D array in shape [r1, c] + matrix1: 2D array in shape [r1, c] MUST be continuous and row-major - set2: 2D array in shape [r2, c] + matrix2: 2D array in shape [r2, c] MUST be continuous and row-major Returns @@ -73,22 +83,22 @@ def compute_correlation(set1, set2): corr_data: 2D array in shape [r1, r2] continuous and row-major in np.float32 """ - set1 = set1.astype(np.float32) - set2 = set2.astype(np.float32) - [r1, d1] = set1.shape - [r2, d2] = set2.shape + 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 - set1 = normalize_for_correlation(set1, 1) - set2 = normalize_for_correlation(set2, 1) + 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, - set2, d2, - set1, d1, + matrix2, d2, + matrix1, d1, 0.0, corr_data, r2) From 530b185a2cc3986f5e6b85590e22a9c46ecef98b Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 20:11:48 -0800 Subject: [PATCH 09/11] add more math derivation --- brainiak/fcma/util.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 12a40f7b9..566aaf89b 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -70,6 +70,19 @@ def compute_correlation(matrix1, matrix2): is the standard deviation of variable X + Reducing the correlation computation to matrix multiplication + and use BLAS GEMM API wrapped by Scipy can speedup the numpy build-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] From 285b1ea70926ac957be515be82a30c3b5180779c Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 22:06:08 -0800 Subject: [PATCH 10/11] more docstring issues --- brainiak/fcma/mvpa_voxelselector.py | 3 +++ 1 file changed, 3 insertions(+) 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 From d0c998452fb66e6e5d050f190ec1687104426f68 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 12 Jan 2017 22:43:52 -0800 Subject: [PATCH 11/11] more review comments --- brainiak/fcma/util.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/brainiak/fcma/util.py b/brainiak/fcma/util.py index 566aaf89b..76f5f885e 100644 --- a/brainiak/fcma/util.py +++ b/brainiak/fcma/util.py @@ -25,7 +25,7 @@ import math -def normalize_for_correlation(data, axis): +def _normalize_for_correlation(data, axis): """normalize the data before computing correlation The data will be z-scored and divided by sqrt(n) @@ -71,7 +71,7 @@ def compute_correlation(matrix1, matrix2): is the standard deviation of variable X Reducing the correlation computation to matrix multiplication - and use BLAS GEMM API wrapped by Scipy can speedup the numpy build-in + and using BLAS GEMM API wrapped by Scipy can speedup the numpy built-in correlation computation (numpy.corrcoef) by one order of magnitude .. math:: @@ -103,8 +103,8 @@ def compute_correlation(matrix1, matrix2): if d1 != d2: raise ValueError('Dimension discrepancy') # preprocess two components - matrix1 = normalize_for_correlation(matrix1, 1) - matrix2 = normalize_for_correlation(matrix2, 1) + 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',