From 08764d3223fabfb3cf51e06e725a30a9c4cd180e Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 16 Aug 2016 19:08:02 -0700 Subject: [PATCH 01/59] add CMake support to the isntallation --- setup.py | 60 ++++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 54 insertions(+), 6 deletions(-) diff --git a/setup.py b/setup.py index b2757fad0..f85d88220 100644 --- a/setup.py +++ b/setup.py @@ -3,6 +3,7 @@ import os import sys import setuptools +import subprocess __version__ = '0.2' @@ -32,6 +33,10 @@ def __str__(self): import pybind11 return pybind11.get_include(self.user) +class CMakeExtension(Extension): + def __init__(self, name, sourcedir=''): + Extension.__init__(self, name, sources=[]) + self.sourcedir = os.path.abspath(sourcedir) ext_modules = [ Extension( @@ -43,6 +48,7 @@ def __str__(self): get_pybind_include(user=True) ], ), + CMakeExtension('brainiak.fcma', 'brainiak/fcma/'), ] @@ -77,16 +83,40 @@ def cpp_flag(compiler): class BuildExt(build_ext): - """A custom build extension for adding compiler-specific options.""" + """A custom build extension for adding compiler-specific options. + Support both regular compiling and CMake compiling. + Support Unix/Linux and MacOS. + """ + if sys.platform == 'Windows': + raise RuntimeError("BrainIAK cannot be built on Windows") + + # configuration for regular compiling c_opts = { 'unix': ['-g0', '-fopenmp'], } - if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', '-ftemplate-depth-1024'] + # configuration for CMake + try: + out = subprocess.check_output(['cmake', '--version']) + except OSError: + raise RuntimeError("CMake must be installed to build the following extensions: " + + ", ".join(e.name for e in self.extensions)) def build_extensions(self): + """the system will execute the run functions of the base class (distutils.command.build_ext) + and get to this function which is override in the sub class + """ + # First, sanity-check the 'extensions' list, i.e. ext_module + self.check_extensions_list(self.extensions) + for ext in self.extensions: + if isinstance(ext, CMakeExtension): + self.cmake_compiling(ext) + elif isinstance(ext, Extension): + self.regular_compiling(ext) + + def regular_compiling(self, ext): ct = self.compiler.compiler_type opts = self.c_opts.get(ct, []) if ct == 'unix': @@ -95,10 +125,28 @@ def build_extensions(self): opts.append(cpp_flag(self.compiler)) if has_flag(self.compiler, '-fvisibility=hidden'): opts.append('-fvisibility=hidden') - for ext in self.extensions: - ext.extra_compile_args = opts - ext.extra_link_args = opts - build_ext.build_extensions(self) + ext.extra_compile_args = opts + ext.extra_link_args = opts + self.build_extension(ext) + + def cmake_compiling(self, ext): + extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) + cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, + '-DPYTHON_EXECUTABLE=' + sys.executable] + + cfg = 'Debug' if self.debug else 'Release' + build_args = ['--config', cfg] + + cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] + build_args += ['--', '-j'] + + env = os.environ.copy() + env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format(env.get('CXXFLAGS', ''), + self.distribution.get_version()) + if not os.path.exists(self.build_temp): + os.makedirs(self.build_temp) + subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, cwd=self.build_temp, env=env) + subprocess.check_call(['cmake', '--build', '.'] + build_args, cwd=self.build_temp) setup( name='brainiak', From 17ab76a1c7cb74280a1bc5e5436c5728b7b8de94 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 17 Aug 2016 14:03:20 -0700 Subject: [PATCH 02/59] initialize fcma --- brainiak/fcma/CMakeLists.txt | 59 ++++++++++++++++++++++++++++++++++++ brainiak/fcma/__init__.py | 14 +++++++++ brainiak/fcma/fcma.py | 40 ++++++++++++++++++++++++ 3 files changed, 113 insertions(+) create mode 100644 brainiak/fcma/CMakeLists.txt create mode 100644 brainiak/fcma/__init__.py create mode 100644 brainiak/fcma/fcma.py diff --git a/brainiak/fcma/CMakeLists.txt b/brainiak/fcma/CMakeLists.txt new file mode 100644 index 000000000..e443676a5 --- /dev/null +++ b/brainiak/fcma/CMakeLists.txt @@ -0,0 +1,59 @@ +cmake_minimum_required (VERSION 2.6) +project (FCMA) + +if ($ENV{VERBOSE}) + set (CMAKE_VERBOSE_MAKEFILE ON) +endif ($ENV{VERBOSE}) + +# add the binary tree to the search path for include files +# so that we will find TutorialConfig.h +include_directories("${PROJECT_BINARY_DIR}") + +link_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") + +set (LD_LIBRARY_PATH "${LD_LIBRARY_PATH}:${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") + +# configure a header file to pass some of the CMake settings +# to the source code +#configure_file ( +# "${PROJECT_SOURCE_DIR}/FCMAConfig.h.in" +# "${PROJECT_BINARY_DIR}/FCMAConfig.h" +# ) + +# should we use MKL? otherwise we use OpenBLAS +#option (USE_MKL "Use Intel MKL" ON) + +set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lz -lhdf5 -lhdf5_hl -lboost_serialization -lboost_iostreams -lmpi") +set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lz -lhdf5 -lhdf5_hl -lboost_serialization -lboost_iostreams -lmpi") +# set (CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") +set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -DCMM_INT_IS_LONG -std=c++14 -fPIC -flto") + +if ($ENV{USE_MKL}) + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_MKL") + set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -mkl=parallel") +else ($ENV{USE_MKL}) + set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lopenblas -lomp") + set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lopenblas -lomp") +endif ($ENV{USE_MKL}) + +set(CMAKE_C_COMPILER ${CC}) +set(CMAKE_CXX_COMPILER ${CXX}) + + +include_directories ("${PROJECT_SOURCE_DIR}/deps/SVM-phi") +add_subdirectory (deps/SVM-phi) +include_directories ("${PROJECT_SOURCE_DIR}/deps/SVM-Simple") +add_subdirectory (deps/SVM-Simple) +include_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/include") +#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -I${CMAKE_INCLUDE_PATH}") # either this or last line +link_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") +set (EXTRA_LIBS ${EXTRA_LIBS} svmphi svmsimple niftiio znz) + +# get all *.cc files recursively +file(GLOB_RECURSE SRC_FILES src/*.cc) + +# add the executable +add_subdirectory(pybind11) +pybind11_add_module(fcma ${SRC_FILES}) +#add_executable (fcma ${SRC_FILES}) +target_link_libraries (fcma ${EXTRA_LIBS}) diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py new file mode 100644 index 000000000..5f9a9900f --- /dev/null +++ b/brainiak/fcma/__init__.py @@ -0,0 +1,14 @@ +# 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""" diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py new file mode 100644 index 000000000..b4bdc96cb --- /dev/null +++ b/brainiak/fcma/fcma.py @@ -0,0 +1,40 @@ +# 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) + +This implementation is based on the following publications: + +.. [Wang2015] Full correlation matrix analysis (FCMA): An unbiased method for + task-related functional connectivity", + Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. + Journal of Neuroscience Methods, 2015. + +.. [Wang2015] "Full correlation matrix analysis of fMRI data on Intel® Xeon + Phi™ coprocessors", + Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, Kai Li, + Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, Theodore L. Willke. + In Proceedings of the International Conference for High Performance Computing, + Networking, Storage and Analysis. 2015. +""" + +# Authors: Yida Wang +# (Intel Labs), 2016 + +import numpy as np +import nibabel as nib + +# __all__ = ['fcma'] + +class fcma: + def run(self): From 1bfccc799dca7bfaad3e32daa1acf98c4f4596ac Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 23 Aug 2016 17:19:10 -0700 Subject: [PATCH 03/59] first path of voxel selection with hard-coded path --- brainiak/fcma/.gitignore | 1 + brainiak/fcma/fcma.py | 299 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 296 insertions(+), 4 deletions(-) create mode 100644 brainiak/fcma/.gitignore diff --git a/brainiak/fcma/.gitignore b/brainiak/fcma/.gitignore new file mode 100644 index 000000000..56fa9649b --- /dev/null +++ b/brainiak/fcma/.gitignore @@ -0,0 +1 @@ +*.npy diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index b4bdc96cb..cc60dc606 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -15,12 +15,12 @@ This implementation is based on the following publications: -.. [Wang2015] Full correlation matrix analysis (FCMA): An unbiased method for +.. [Wang2015-1] Full correlation matrix analysis (FCMA): An unbiased method for task-related functional connectivity", Yida Wang, Jonathan D Cohen, Kai Li, Nicholas B Turk-Browne. Journal of Neuroscience Methods, 2015. -.. [Wang2015] "Full correlation matrix analysis of fMRI data on Intel® Xeon +.. [Wang2015-2] "Full correlation matrix analysis of fMRI data on Intel® Xeon Phi™ coprocessors", Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, Kai Li, Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, Theodore L. Willke. @@ -33,8 +33,299 @@ import numpy as np import nibabel as nib +import os, math, sys, time +from mpi4py import MPI +from nilearn.input_data import NiftiMasker +from scipy.stats.mstats import zscore +import ctypes +from ctypes.util import find_library +from ctypes import byref, c_char, c_float, c_int +from sklearn import cross_validation +from sklearn import svm -# __all__ = ['fcma'] -class fcma: +# __all__ = [''] + +WORKTAG = 0 +TERMINATETAG = 1 + +def readActivityData(dir, file_extension, mask_file): + """ read in data in NIfTI format and apply the spacial mask to them + :param dir: the path to all input file + :param file_extension: the file extension, usually nii.gz or nii + :param mask_file: the absolute path of the mask file, we apply the mask right after reading a file for saving memory + :return: activity_data: array of matrices in (nTRs, nVoxels) shape (numpy array), + len(activity_data) equals the number of subjects + """ + mask_img = nib.load(mask_file) + nifti_masker = NiftiMasker(mask_img=mask_img) + files = [f for f in os.listdir(dir) if os.path.isfile(os.path.join(dir, f)) and f.endswith(file_extension)] + activity_data = [] + for f in files: + # NOTE: nifti_masker.fit_transform runs one order of magnitude slower than C++ version + # it returns data in (nTRs, nVoxels) shape + img = nib.load(os.path.join(dir, f)) + masked_data = nifti_masker.fit_transform(img) + activity_data.append(masked_data) + print(f, masked_data.shape) + sys.stdout.flush() + return activity_data + +def separateEpochs(activity_data, epoch_map): + """ keep data in epochs of interest specified in epoch_map and z-score them for computing correlation + :param activity_data: array of matrices in (nTRs, nVoxels) shape, + consisting of the masked activity data of all subjects + :param epoch_map: array of (subject id, starting TR, ending TR) tuple + :return: raw_data: array of matrices in (epoch length, nVoxels) shape, + len(raw_data) equals the number of epochs + """ + raw_data = [] + for (sid, s, e) in epoch_map: + mat = activity_data[sid][s:e+1,:] + (r, c) = mat.shape + mat = zscore(mat, axis=0, ddof=0) + mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero + mat = mat / math.sqrt(r) + raw_data.append(mat) + return raw_data + + +def preprocessActivityData(): + """Read, filter (applying masks and eliminating TRs not in epochs of interest) the activity data + + """ + +class VoxelSelector: + """Correlation-based voxel selection component of FCMA + Parameters + ---------- + + raw_data : numpy array of array of float + Each array of float contains the activity data of an epoch of a subject, + which essentially is a nVoxels x nTRs matrix but is serialized as one dimensional vector + Assumption: 1. all activity data contains the same number of voxels + 2. the activity data has been z-scored, ready to compute correlation as matrix multiplication + 3. the epochs belonging to the same subject are continuous in the array + 4. all subjects have the same number of epochs + 5. voxel selection is always done in the self-correlation, i.e. raw_data correlate with themselves + 6. the classifier used is SVM with linear kernel + + epochs_per_subj : int + The number of epochs of each subject + + num_voxels: int + The number of voxels participating in the voxel selection + + labels: 1D array of int + Describing the condition labels of the epochs, the length of labels equals the number of epochs, + i.e. the length of raw_data + + num_folds: int + The number of folds to be conducted in the cross validation + + voxel_unit: int, default 100 + The number of voxel assigned to a worker each time + """ + def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100): + self.raw_data = raw_data + self.epochs_per_subj = epochs_per_subj + self.num_voxels = raw_data[0].shape[1] + self.labels = labels + self.num_folds = num_folds + self.voxel_unit = voxel_unit + libblas_so_name = find_library('blas') if find_library('blas') else find_library('mkl_rt') + if libblas_so_name is None: + raise RuntimeError("No Blas library is found in the system") + self.blas_library = ctypes.cdll.LoadLibrary(libblas_so_name) + if self.num_voxels == 0: + raise RuntimeError("Zero processed voxels") + def run(self): + """run voxel selection in master-worker model + Sort the voxels based on the cross-validation accuracy of their correlation vectors + """ + rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + results = self.master() + else: + self.worker() + # TODO: sort results + results.sort(key=lambda tup: tup[1], reverse=True) + return results + + def master(self): + """ + return value + results: an array of tuple (voxel id, accuracy), the length of array equals the number of voxels + """ + results = [] + comm = MPI.COMM_WORLD + size = comm.Get_size() + sending_voxels = self.voxel_unit if self.voxel_unit Date: Tue, 23 Aug 2016 17:58:25 -0700 Subject: [PATCH 04/59] sort files before processing; only process the results in master --- brainiak/fcma/fcma.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index cc60dc606..525f678db 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -59,7 +59,7 @@ def readActivityData(dir, file_extension, mask_file): """ mask_img = nib.load(mask_file) nifti_masker = NiftiMasker(mask_img=mask_img) - files = [f for f in os.listdir(dir) if os.path.isfile(os.path.join(dir, f)) and f.endswith(file_extension)] + files = [f for f in sorted(os.listdir(dir)) if os.path.isfile(os.path.join(dir, f)) and f.endswith(file_extension)] activity_data = [] for f in files: # NOTE: nifti_masker.fit_transform runs one order of magnitude slower than C++ version @@ -147,10 +147,11 @@ def run(self): rank = MPI.COMM_WORLD.Get_rank() if rank == 0: results = self.master() + # TODO: sort results + results.sort(key=lambda tup: tup[1], reverse=True) else: self.worker() - # TODO: sort results - results.sort(key=lambda tup: tup[1], reverse=True) + results = [] return results def master(self): @@ -192,7 +193,7 @@ def master(self): results += result for i in range(1, size): - comm.send(None, dest=i, tag=DIETAG) + comm.send(None, dest=i, tag=TERMINATETAG) return results @@ -328,4 +329,5 @@ def voxelScoring(self, task): labels = comm.bcast(labels, root=0) vs = VoxelSelector(raw_data, 12, labels, 18) results = vs.run() - print(results[0:100]) \ No newline at end of file + if rank==0: + print(results[0:100]) From 570bd1525722ac9d4348f6c1ec51b2faad6e03e4 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 24 Aug 2016 00:25:18 -0700 Subject: [PATCH 05/59] use self-defined mask applying method, instead of NiftiMasker --- brainiak/fcma/fcma.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 525f678db..fff60093f 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -35,7 +35,6 @@ import nibabel as nib import os, math, sys, time from mpi4py import MPI -from nilearn.input_data import NiftiMasker from scipy.stats.mstats import zscore import ctypes from ctypes.util import find_library @@ -58,14 +57,23 @@ def readActivityData(dir, file_extension, mask_file): len(activity_data) equals the number of subjects """ mask_img = nib.load(mask_file) - nifti_masker = NiftiMasker(mask_img=mask_img) + mask = mask_img.get_data() + count = 0 + for index in np.ndindex(mask.shape): + if mask[index]!=0: + count+=1 files = [f for f in sorted(os.listdir(dir)) if os.path.isfile(os.path.join(dir, f)) and f.endswith(file_extension)] activity_data = [] for f in files: - # NOTE: nifti_masker.fit_transform runs one order of magnitude slower than C++ version - # it returns data in (nTRs, nVoxels) shape img = nib.load(os.path.join(dir, f)) - masked_data = nifti_masker.fit_transform(img) + data = img.get_data() + (d1, d2, d3, d4) = data.shape + masked_data = np.zeros([d4, count], np.float32, order='C') + count1=0 + for index in np.ndindex(mask.shape): + if mask[index]!=0: + masked_data[:,count1] = np.copy(data[index]) + count1+=1 activity_data.append(masked_data) print(f, masked_data.shape) sys.stdout.flush() @@ -298,6 +306,7 @@ def voxelScoring(self, task): results = self.crossValidation(task, corr) time2 = time.time() print('task:', task[0], time2-time1) + sys.stdout.flush() return results if __name__ == '__main__': From 0c19f91850150f1f2f83a777b18832cb3d875fde Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 24 Aug 2016 16:56:24 -0700 Subject: [PATCH 06/59] add Fisher tranformation to the within subject normalization --- brainiak/fcma/fcma.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index fff60093f..978107595 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -216,7 +216,7 @@ def worker(self): def correlationComputation(self, task): s = task[0] e = s+task[1] - #time1 = time.time() + time1 = time.time() nEpochs = len(self.raw_data) corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') count=0 @@ -235,28 +235,31 @@ def correlationComputation(self, task): mat.ctypes.data_as(ctypes.c_void_p), n2, zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) count += 1 - #time2 = time.time() - #print('correlation computation time:', (time2-time1)*1000.0) - #sys.stdout.flush() + time2 = time.time() + print('correlation computation time:', (time2-time1)*1000.0) + sys.stdout.flush() return corr def correlationNormalization(self, corr): - #time1 = time.time() + time1 = time.time() (sv, e, av) = corr.shape for i in range(sv): start = 0 while start Date: Wed, 24 Aug 2016 17:12:34 -0700 Subject: [PATCH 07/59] fix a bug in Fisher Transformation --- brainiak/fcma/fcma.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 978107595..d6a9cf77a 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -248,7 +248,7 @@ def correlationNormalization(self, corr): while start Date: Thu, 25 Aug 2016 16:04:15 -0700 Subject: [PATCH 08/59] use C++ to do within-subject normalization, add running example --- .gitmodules | 3 + brainiak/fcma/CMakeLists.txt | 51 ++++------------- brainiak/fcma/fcma.py | 44 ++------------- brainiak/fcma/pybind11 | 1 + brainiak/fcma/src/fcma_extension.cc | 88 +++++++++++++++++++++++++++++ examples/fcma/run.py | 35 ++++++++++++ setup.py | 2 +- 7 files changed, 144 insertions(+), 80 deletions(-) create mode 100644 .gitmodules create mode 160000 brainiak/fcma/pybind11 create mode 100644 brainiak/fcma/src/fcma_extension.cc create mode 100644 examples/fcma/run.py diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 000000000..d19b82e8b --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "brainiak/fcma/pybind11"] + path = brainiak/fcma/pybind11 + url = https://github.com/pybind/pybind11.git diff --git a/brainiak/fcma/CMakeLists.txt b/brainiak/fcma/CMakeLists.txt index e443676a5..1b8e6bec1 100644 --- a/brainiak/fcma/CMakeLists.txt +++ b/brainiak/fcma/CMakeLists.txt @@ -1,59 +1,32 @@ cmake_minimum_required (VERSION 2.6) -project (FCMA) +project (FCMA_EXT) if ($ENV{VERBOSE}) set (CMAKE_VERBOSE_MAKEFILE ON) endif ($ENV{VERBOSE}) # add the binary tree to the search path for include files -# so that we will find TutorialConfig.h include_directories("${PROJECT_BINARY_DIR}") -link_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") - -set (LD_LIBRARY_PATH "${LD_LIBRARY_PATH}:${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") - -# configure a header file to pass some of the CMake settings -# to the source code -#configure_file ( -# "${PROJECT_SOURCE_DIR}/FCMAConfig.h.in" -# "${PROJECT_BINARY_DIR}/FCMAConfig.h" -# ) - -# should we use MKL? otherwise we use OpenBLAS -#option (USE_MKL "Use Intel MKL" ON) - -set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lz -lhdf5 -lhdf5_hl -lboost_serialization -lboost_iostreams -lmpi") -set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lz -lhdf5 -lhdf5_hl -lboost_serialization -lboost_iostreams -lmpi") +#set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") +set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") # set (CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -DCMM_INT_IS_LONG -std=c++14 -fPIC -flto") -if ($ENV{USE_MKL}) - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_MKL") - set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -mkl=parallel") -else ($ENV{USE_MKL}) - set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lopenblas -lomp") - set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lopenblas -lomp") -endif ($ENV{USE_MKL}) +#if ($ENV{USE_MKL}) +# set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_MKL") +# set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -mkl=parallel") +#else ($ENV{USE_MKL}) +# set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lopenblas -lomp") +# set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lopenblas -lomp") +#endif ($ENV{USE_MKL}) set(CMAKE_C_COMPILER ${CC}) set(CMAKE_CXX_COMPILER ${CXX}) - -include_directories ("${PROJECT_SOURCE_DIR}/deps/SVM-phi") -add_subdirectory (deps/SVM-phi) -include_directories ("${PROJECT_SOURCE_DIR}/deps/SVM-Simple") -add_subdirectory (deps/SVM-Simple) -include_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/include") -#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -I${CMAKE_INCLUDE_PATH}") # either this or last line -link_directories ("${PROJECT_SOURCE_DIR}/deps/nifticlib-2.0.0/lib") -set (EXTRA_LIBS ${EXTRA_LIBS} svmphi svmsimple niftiio znz) - # get all *.cc files recursively file(GLOB_RECURSE SRC_FILES src/*.cc) -# add the executable +# add the library add_subdirectory(pybind11) -pybind11_add_module(fcma ${SRC_FILES}) -#add_executable (fcma ${SRC_FILES}) -target_link_libraries (fcma ${EXTRA_LIBS}) +pybind11_add_module(fcma_extension ${SRC_FILES}) \ No newline at end of file diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index d6a9cf77a..1c5b61ab5 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -41,9 +41,10 @@ from ctypes import byref, c_char, c_float, c_int from sklearn import cross_validation from sklearn import svm +from . import fcma_extension -# __all__ = [''] +#__all__ = ['readActivityData', 'separateEpochs', 'VoxelSelector'] WORKTAG = 0 TERMINATETAG = 1 @@ -97,12 +98,6 @@ def separateEpochs(activity_data, epoch_map): raw_data.append(mat) return raw_data - -def preprocessActivityData(): - """Read, filter (applying masks and eliminating TRs not in epochs of interest) the activity data - - """ - class VoxelSelector: """Correlation-based voxel selection component of FCMA Parameters @@ -305,41 +300,10 @@ def voxelScoring(self, task): corr = self.correlationComputation(task) # corr is a 3D array in row major, # in (selected_voxels, epochs, all_voxels) shape # corr[i,e,s+j] = corr[j,e,s+i] - corr = self.correlationNormalization(corr) # in-place z-score, the result is still in corr + #corr = self.correlationNormalization(corr) # in-place z-score, the result is still in corr + fcma_extension.normalization(corr, self.epochs_per_subj) results = self.crossValidation(task, corr) time2 = time.time() print('task:', task[0], time2-time1) sys.stdout.flush() return results - -if __name__ == '__main__': - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - sys.stdout.flush() - raw_data = [] - labels = [] - if rank==0: - activity_data = readActivityData('/Users/yidawang/data/face_scene/raw', 'nii.gz', - '/Users/yidawang/data/face_scene/mask.nii.gz') - epoch_map = [] - for i in range(18): - epoch_map.append((i, 4, 15)) - epoch_map.append((i, 24, 35)) - epoch_map.append((i, 44, 55)) - epoch_map.append((i, 64, 75)) - epoch_map.append((i, 84, 95)) - epoch_map.append((i, 104, 115)) - epoch_map.append((i, 124, 135)) - epoch_map.append((i, 144, 155)) - epoch_map.append((i, 164, 175)) - epoch_map.append((i, 184, 195)) - epoch_map.append((i, 204, 215)) - epoch_map.append((i, 224, 235)) - raw_data=separateEpochs(activity_data, epoch_map) - labels = [i%2 for i in range(216)] - raw_data = comm.bcast(raw_data, root=0) - labels = comm.bcast(labels, root=0) - vs = VoxelSelector(raw_data, 12, labels, 18) - results = vs.run() - if rank==0: - print(results[0:100]) diff --git a/brainiak/fcma/pybind11 b/brainiak/fcma/pybind11 new file mode 160000 index 000000000..69b624667 --- /dev/null +++ b/brainiak/fcma/pybind11 @@ -0,0 +1 @@ +Subproject commit 69b6246677e43b8e2c2ee6be6970d2b5dbb156c6 diff --git a/brainiak/fcma/src/fcma_extension.cc b/brainiak/fcma/src/fcma_extension.cc new file mode 100644 index 000000000..0522ea53e --- /dev/null +++ b/brainiak/fcma/src/fcma_extension.cc @@ -0,0 +1,88 @@ +/* +// 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. +// */ +// +#include "Python.h" +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +void within_subject_norm_native(py::array_t py_data, int epochsPerSubj) +{ + // the input data is a 3D array in row major, in (selected_voxels, epochs, all_voxels) shape + py::buffer_info buf = py_data.request(); + float* data = (float*)buf.ptr; + // sanity check + if (buf.ndim != 3) + throw std::runtime_error("The multi-subject correlation data structure must be 3D"); + size_t nSelectedVoxels = buf.shape[0]; + size_t nEpochs = buf.shape[1]; + size_t nVoxels = buf.shape[2]; + size_t nSubjs = nEpochs / epochsPerSubj; + //array_t inArray = {(float *)inBuf.ptr, (int)inBuf.strides[0] >> 2, + // (int)inBuf.strides[1] >> 2, (int)inBuf.strides[2] >> 2, + // (int)inBuf.shape[3]}; + #pragma omp parallel for + for(int v = 0 ; v < nSelectedVoxels*nSubjs ; v++) + { + int s = v % nSubjs; // subject id + int i = v / nSubjs; // voxel id + //float (*mat)[nVoxels] = (float(*)[nVoxels])&(data[i*nEpochs*nVoxels]); + float *mat = data+i*nEpochs*nVoxels; + #pragma simd + for(int j = 0 ; j < nVoxels ; j++) + { + float mean = 0.0f; + float std_dev = 0.0f; + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + float num = 1.0f + mat[b*nVoxels+j]; + float den = 1.0f - mat[b*nVoxels+j]; + num = (num <= 0.0f) ? 1e-4 : num; + den = (den <= 0.0f) ? 1e-4 : den; + mat[b*nVoxels+j] = 0.5f * logf(num/den); + mean += mat[b*nVoxels+j]; + std_dev += mat[b*nVoxels+j] * mat[b*nVoxels+j]; + } + mean = mean / (float)epochsPerSubj; + std_dev = std_dev / (float)epochsPerSubj - mean*mean; + float inv_std_dev = (std_dev <= 0.0f) ? 0.0f : 1.0f / sqrt(std_dev); + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + mat[b*nVoxels+j] = (mat[b*nVoxels+j] - mean) * inv_std_dev; + } + } + } + return; +} + +PYBIND11_PLUGIN(fcma_extension) { + py::module m("fcma_extension"); + m.def("normalization", &within_subject_norm_native, "Within-subejct correlation normalization"); + +#ifdef VERSION_INFO + m.attr("__version__") = py::str(VERSION_INFO); +#else + m.attr("__version__") = py::str("dev"); +#endif + return m.ptr(); +} \ No newline at end of file diff --git a/examples/fcma/run.py b/examples/fcma/run.py new file mode 100644 index 000000000..eff21c706 --- /dev/null +++ b/examples/fcma/run.py @@ -0,0 +1,35 @@ +from brainiak.fcma.fcma import * +from mpi4py import MPI +import sys + +# mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz +if __name__ == '__main__': + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + sys.stdout.flush() + raw_data = [] + labels = [] + if rank==0: + activity_data = readActivityData(sys.argv[1], sys.argv[2], sys.argv[3]) + epoch_map = [] + for i in range(18): + epoch_map.append((i, 4, 15)) + epoch_map.append((i, 24, 35)) + epoch_map.append((i, 44, 55)) + epoch_map.append((i, 64, 75)) + epoch_map.append((i, 84, 95)) + epoch_map.append((i, 104, 115)) + epoch_map.append((i, 124, 135)) + epoch_map.append((i, 144, 155)) + epoch_map.append((i, 164, 175)) + epoch_map.append((i, 184, 195)) + epoch_map.append((i, 204, 215)) + epoch_map.append((i, 224, 235)) + raw_data=separateEpochs(activity_data, epoch_map) + labels = [i%2 for i in range(216)] + raw_data = comm.bcast(raw_data, root=0) + labels = comm.bcast(labels, root=0) + vs = VoxelSelector(raw_data, 12, labels, 18) + results = vs.run() + if rank==0: + print(results[0:100]) diff --git a/setup.py b/setup.py index f85d88220..c233c6383 100644 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ def __init__(self, name, sourcedir=''): get_pybind_include(user=True) ], ), - CMakeExtension('brainiak.fcma', 'brainiak/fcma/'), + CMakeExtension('brainiak.fcma.fcma_extension', 'brainiak/fcma/'), ] From 1bb52823d443a11b3c9a62fa2a2cff5d89bbdc32 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 25 Aug 2016 16:29:20 -0700 Subject: [PATCH 09/59] check compiler type in CMake --- brainiak/fcma/CMakeLists.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/brainiak/fcma/CMakeLists.txt b/brainiak/fcma/CMakeLists.txt index 1b8e6bec1..5858b129c 100644 --- a/brainiak/fcma/CMakeLists.txt +++ b/brainiak/fcma/CMakeLists.txt @@ -9,7 +9,7 @@ endif ($ENV{VERBOSE}) include_directories("${PROJECT_BINARY_DIR}") #set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") -set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") +#set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") # set (CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -DCMM_INT_IS_LONG -std=c++14 -fPIC -flto") @@ -24,6 +24,10 @@ set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -DCMM_INT_IS_LONG -std=c++ set(CMAKE_C_COMPILER ${CC}) set(CMAKE_CXX_COMPILER ${CXX}) +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "clang") + set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") +endif() + # get all *.cc files recursively file(GLOB_RECURSE SRC_FILES src/*.cc) From 7cbad638cdb12a831aa358640992e4808a71ebc2 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 25 Aug 2016 18:40:46 -0700 Subject: [PATCH 10/59] enable finding blas library/MKL in different OSes --- brainiak/fcma/CMakeLists.txt | 4 ++- brainiak/fcma/fcma.py | 55 ++++++++++++++++-------------------- 2 files changed, 27 insertions(+), 32 deletions(-) diff --git a/brainiak/fcma/CMakeLists.txt b/brainiak/fcma/CMakeLists.txt index 5858b129c..e675bdf1d 100644 --- a/brainiak/fcma/CMakeLists.txt +++ b/brainiak/fcma/CMakeLists.txt @@ -11,7 +11,7 @@ include_directories("${PROJECT_BINARY_DIR}") #set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") #set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") # set (CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") -set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -DCMM_INT_IS_LONG -std=c++14 -fPIC -flto") +set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -std=c++14 -fPIC -flto") #if ($ENV{USE_MKL}) # set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_MKL") @@ -26,6 +26,8 @@ set(CMAKE_CXX_COMPILER ${CXX}) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "clang") set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lgomp") endif() # get all *.cc files recursively diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 1c5b61ab5..e62ee1f7b 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -136,10 +136,23 @@ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100) self.labels = labels self.num_folds = num_folds self.voxel_unit = voxel_unit - libblas_so_name = find_library('blas') if find_library('blas') else find_library('mkl_rt') - if libblas_so_name is None: - raise RuntimeError("No Blas library is found in the system") - self.blas_library = ctypes.cdll.LoadLibrary(libblas_so_name) + #libblas_so_name = find_library('blas') if find_library('blas') else find_library('mkl_rt') + #if libblas_so_name is None: + # raise RuntimeError("No Blas library is found in the system") + # use the following code since find_library doesn't work for libraries in customized directory in Linux + if sys.platform=='darwin': + extension = '.dylib' + elif sys.platform=='linux': + extension = '.so' + else: + raise RuntimeError("Unsupported operating system") + try: + self.blas_library = ctypes.cdll.LoadLibrary('libblas'+extension) + except: + try: + self.blas_library = ctypes.cdll.LoadLibrary('libmkl_rt'+extension) + except: + raise RuntimeError("No Blas library is found in the system") if self.num_voxels == 0: raise RuntimeError("Zero processed voxels") @@ -211,7 +224,6 @@ def worker(self): def correlationComputation(self, task): s = task[0] e = s+task[1] - time1 = time.time() nEpochs = len(self.raw_data) corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') count=0 @@ -230,13 +242,9 @@ def correlationComputation(self, task): mat.ctypes.data_as(ctypes.c_void_p), n2, zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) count += 1 - time2 = time.time() - print('correlation computation time:', (time2-time1)*1000.0) - sys.stdout.flush() return corr def correlationNormalization(self, corr): - time1 = time.time() (sv, e, av) = corr.shape for i in range(sv): start = 0 @@ -248,13 +256,9 @@ def correlationNormalization(self, corr): zscore(corr[i, start:start+self.epochs_per_subj, :], axis = 0, ddof = 0) start += self.epochs_per_subj corr = np.nan_to_num(corr) # if zscore fails (standard deviation is zero), set all values to be zero - time2 = time.time() - print('normalization time:', (time2-time1)*1000.0) - sys.stdout.flush() return corr def crossValidation(self, task, corr): - time1 = time.time() (sv, e, av) = corr.shape kernel_matrix = np.zeros((e, e), np.float32, order='C') results = [] @@ -272,38 +276,27 @@ def crossValidation(self, task, corr): for j in range(kernel_matrix.shape[0]): for k in range(j): kernel_matrix[k,j] = kernel_matrix[j,k] - clf = svm.SVC(kernel='precomputed', shrinking=False) - scores = cross_validation.cross_val_score(clf, kernel_matrix, self.labels, cv=18) # by default, no shuffling + clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) + skf = cross_validation.StratifiedKFold(self.labels, n_folds=self.num_folds, shuffle=False) # no shuffling in cv + scores = cross_validation.cross_val_score(clf, kernel_matrix, self.labels, cv=skf) results.append((i+task[0], scores.mean())) - time2 = time.time() - print('cross validation time:', (time2-time1)*1000.0) - #print(results) - sys.stdout.flush() return results def voxelScoring(self, task): """Voxel selection worker Take the task in, do analysis on voxels specified by the task (voxel id, number of voxels) """ - # TODO: use the C++ routines to do the following: - # 1. correlation computation - # 2. normalization - # 3. cross validation - # input: activity matrices, starting voxel, number of voxels, epochs_per_subj, labels, number of folds - # output: array of float, accuracy numbers - # scores = voxelScoringingCPP(self.raw_data, task[0], task[1], self.epochs_per_subj, self.labels, self.num_folds) - #result = () - #for i in range(len(scores)): - # result.append((task[0]+i, scores[i])) - #return result time1 = time.time() + # 1. correlation computation corr = self.correlationComputation(task) # corr is a 3D array in row major, # in (selected_voxels, epochs, all_voxels) shape # corr[i,e,s+j] = corr[j,e,s+i] + # 2. normalization #corr = self.correlationNormalization(corr) # in-place z-score, the result is still in corr fcma_extension.normalization(corr, self.epochs_per_subj) + # 3. cross validation results = self.crossValidation(task, corr) time2 = time.time() - print('task:', task[0], time2-time1) + print('task:', task[0]%self.num_voxels, time2-time1) sys.stdout.flush() return results From b9f4d9cc02d536e2b39d0193b040e6a8c41e4100 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 26 Aug 2016 13:00:00 -0700 Subject: [PATCH 11/59] use new data structure to depict epoch infomation, add an example epoch cube list --- brainiak/fcma/fcma.py | 30 ++++++++++++++++--------- examples/fcma/data/fs_epoch_labels.npy | Bin 0 -> 105504 bytes examples/fcma/run.py | 18 ++------------- 3 files changed, 22 insertions(+), 26 deletions(-) create mode 100644 examples/fcma/data/fs_epoch_labels.npy diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index e62ee1f7b..fa9581462 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -51,6 +51,7 @@ def readActivityData(dir, file_extension, mask_file): """ read in data in NIfTI format and apply the spacial mask to them + :param dir: the path to all input file :param file_extension: the file extension, usually nii.gz or nii :param mask_file: the absolute path of the mask file, we apply the mask right after reading a file for saving memory @@ -80,26 +81,35 @@ def readActivityData(dir, file_extension, mask_file): sys.stdout.flush() return activity_data -def separateEpochs(activity_data, epoch_map): +def separateEpochs(activity_data, epoch_list): """ keep data in epochs of interest specified in epoch_map and z-score them for computing correlation + :param activity_data: array of matrices in (nTRs, nVoxels) shape, consisting of the masked activity data of all subjects - :param epoch_map: array of (subject id, starting TR, ending TR) tuple + :param epoch_list: array of (subject id, starting TR, ending TR) tuple :return: raw_data: array of matrices in (epoch length, nVoxels) shape, len(raw_data) equals the number of epochs """ raw_data = [] - for (sid, s, e) in epoch_map: - mat = activity_data[sid][s:e+1,:] - (r, c) = mat.shape - mat = zscore(mat, axis=0, ddof=0) - mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero - mat = mat / math.sqrt(r) - raw_data.append(mat) - return raw_data + labels = [] + for sid in range(len(epoch_list)): + epoch = epoch_list[sid] # epoch is a numpy array in shape (condition, nEpochs, nTRs) + for cond in range(epoch.shape[0]): + sub_epoch = epoch[cond,:,:] + for eid in range(epoch.shape[1]): + r = np.sum(sub_epoch[eid,:]) + if r > 0: # there is an epoch in this condition + mat = activity_data[sid][sub_epoch[eid,:]==1,:] + mat = zscore(mat, axis=0, ddof=0) + mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero + mat = mat / math.sqrt(r) + raw_data.append(mat) + labels.append(cond) + return raw_data, labels class VoxelSelector: """Correlation-based voxel selection component of FCMA + Parameters ---------- diff --git a/examples/fcma/data/fs_epoch_labels.npy b/examples/fcma/data/fs_epoch_labels.npy new file mode 100644 index 0000000000000000000000000000000000000000..0602190b2beabf7a5cee2df4dcaedce89d8ee062 GIT binary patch literal 105504 zcmeIwv5KB$6ot{G^(v+tg)|8%#0#)fY!s{{n23THN$k|C@Plw-IG5bU%z?(KNjGDgDk*9)-dzHKo;O3YnWRBPFk7GJrA4j-md1MmD$|$u=(!oY93md&3*IG z_M4Sl%>x5jfQPJM=7E7Mz(dwB^T5ER04JA9JKylI^VXM3dAL;C`G$v`x4vA;!==*B zH_b!a@04+A7T_Ukn0a6z3-FLN%seoV1$f9B=2n1{R%Ub0!{)oUt9fW;HupSizI(fx zhgN2D-#oPaW+hkiz(5w@A#0d*Ar950J?R?WbwEa#Qmu3MTvWA%l2C@JTS;Nc&16hEFtYK~iIB8`z_dIOA zd%K#4R%Ub0!{)oUt9fW;Huud#+izBKH4hAA0Uol3nFj{401sKi%mV|L0-Rhb?R>+- z&RbtD<>69k=Nle&-uiMW50^?i-!uEWkt7F!R7b7T_Ukm|Fo( zTA9r~51a4auI8bY+1&H6`R?s%9$J~ree=-vo0VM60|Qxrhpb`dfq^W*L)I|!z`&&d zCznb)-|(>W)|X3pxK!HthKHTEzFf+~rP9tf%|qMolyPYm;2~?6d0-$5@Q^jkJTQ<2 zc*q*&R)CXMW^>QO=DWA6d1z%e_dIOAd%K#4R%UbGJhc60C0FynKo;O3YnXXpAPex2 zHOxFPa4EparP9tfJnX#nQO=DWA6d1z%e_sv7wZ&q?O z4-8}h99oA2JP=Ao6@-1D&c?(J$GTA9s# z^U(I2m0ZmO16hEFtYPMXfh@p7)-dzHz@-2umr6U|@UZjNmrHrLRNDE5hn=^+T*||x K(#|*k=iy(7!Xh^S literal 0 HcmV?d00001 diff --git a/examples/fcma/run.py b/examples/fcma/run.py index eff21c706..2717658be 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -11,22 +11,8 @@ labels = [] if rank==0: activity_data = readActivityData(sys.argv[1], sys.argv[2], sys.argv[3]) - epoch_map = [] - for i in range(18): - epoch_map.append((i, 4, 15)) - epoch_map.append((i, 24, 35)) - epoch_map.append((i, 44, 55)) - epoch_map.append((i, 64, 75)) - epoch_map.append((i, 84, 95)) - epoch_map.append((i, 104, 115)) - epoch_map.append((i, 124, 135)) - epoch_map.append((i, 144, 155)) - epoch_map.append((i, 164, 175)) - epoch_map.append((i, 184, 195)) - epoch_map.append((i, 204, 215)) - epoch_map.append((i, 224, 235)) - raw_data=separateEpochs(activity_data, epoch_map) - labels = [i%2 for i in range(216)] + epoch_map = np.load('data/fs_epoch_labels.npy') # a list of numpy array in shape (condition, nEpochs, nTRs) + raw_data, labels=separateEpochs(activity_data, epoch_map) raw_data = comm.bcast(raw_data, root=0) labels = comm.bcast(labels, root=0) vs = VoxelSelector(raw_data, 12, labels, 18) From 54da6f0d42fe845d87fdcb2073e7d80d7c0a8365 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 26 Aug 2016 14:47:48 -0700 Subject: [PATCH 12/59] simplify the fcma example, encapsulate MPI communicatation --- brainiak/fcma/fcma.py | 25 +++++++++++++++++++++++++ examples/fcma/run.py | 23 +++++++++-------------- 2 files changed, 34 insertions(+), 14 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index fa9581462..81d4a7e55 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -58,6 +58,7 @@ def readActivityData(dir, file_extension, mask_file): :return: activity_data: array of matrices in (nTRs, nVoxels) shape (numpy array), len(activity_data) equals the number of subjects """ + time1 = time.time() mask_img = nib.load(mask_file) mask = mask_img.get_data() count = 0 @@ -79,6 +80,9 @@ def readActivityData(dir, file_extension, mask_file): activity_data.append(masked_data) print(f, masked_data.shape) sys.stdout.flush() + time2 = time.time() + print('data reading done, takes', time2-time1, 's') + sys.stdout.flush() return activity_data def separateEpochs(activity_data, epoch_list): @@ -90,6 +94,7 @@ def separateEpochs(activity_data, epoch_list): :return: raw_data: array of matrices in (epoch length, nVoxels) shape, len(raw_data) equals the number of epochs """ + time1 = time.time() raw_data = [] labels = [] for sid in range(len(epoch_list)): @@ -105,6 +110,26 @@ def separateEpochs(activity_data, epoch_list): mat = mat / math.sqrt(r) raw_data.append(mat) labels.append(cond) + time2 = time.time() + print('epoch separation done, takes', time2-time1, 's') + sys.stdout.flush() + return raw_data, labels + +def prepareData(data_dir, extension, mask_file, epoch_file): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + raw_data = [] + labels = [] + if rank==0: + activity_data = readActivityData(data_dir, extension, mask_file) + epoch_list = np.load(epoch_file) # a list of numpy array in shape (condition, nEpochs, nTRs) + raw_data, labels=separateEpochs(activity_data, epoch_list) + time1 = time.time() + raw_data = comm.bcast(raw_data, root=0) + labels = comm.bcast(labels, root=0) + if rank == 0: + time2 = time.time() + print('data broadcasting done, takes', time2-time1, 's') return raw_data, labels class VoxelSelector: diff --git a/examples/fcma/run.py b/examples/fcma/run.py index 2717658be..0efe3b63c 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -1,21 +1,16 @@ from brainiak.fcma.fcma import * -from mpi4py import MPI import sys -# mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz +# mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz data/fs_epoch_labels.npy if __name__ == '__main__': - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - sys.stdout.flush() - raw_data = [] - labels = [] - if rank==0: - activity_data = readActivityData(sys.argv[1], sys.argv[2], sys.argv[3]) - epoch_map = np.load('data/fs_epoch_labels.npy') # a list of numpy array in shape (condition, nEpochs, nTRs) - raw_data, labels=separateEpochs(activity_data, epoch_map) - raw_data = comm.bcast(raw_data, root=0) - labels = comm.bcast(labels, root=0) - vs = VoxelSelector(raw_data, 12, labels, 18) + data_dir = sys.argv[1] + extension = sys.argv[2] + mask_file = sys.argv[3] + epoch_file = sys.argv[4] + raw_data, labels = prepareData(data_dir, extension, mask_file, epoch_file) + epochs_per_subj = 12 + num_subjs = 18 + vs = VoxelSelector(raw_data, epochs_per_subj, labels, num_subjs) results = vs.run() if rank==0: print(results[0:100]) From d9ea11270b4f35fc3b14b007ca8aa7befa6718de Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 26 Aug 2016 16:35:51 -0700 Subject: [PATCH 13/59] broadcast the data subject by subject to prevent size overflow --- brainiak/fcma/fcma.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 81d4a7e55..cb6c5d28d 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -118,14 +118,19 @@ def separateEpochs(activity_data, epoch_list): def prepareData(data_dir, extension, mask_file, epoch_file): comm = MPI.COMM_WORLD rank = comm.Get_rank() - raw_data = [] labels = [] + raw_data = [] if rank==0: activity_data = readActivityData(data_dir, extension, mask_file) epoch_list = np.load(epoch_file) # a list of numpy array in shape (condition, nEpochs, nTRs) raw_data, labels=separateEpochs(activity_data, epoch_list) time1 = time.time() - raw_data = comm.bcast(raw_data, root=0) + raw_data_length = len(raw_data) + raw_data_length = comm.bcast(raw_data_length, root=0) + for i in range(raw_data_length): # broadcast the data subject by subject to prevent size overflow + if rank!=0: + raw_data.append(None) + raw_data[i] = comm.bcast(raw_data[i], root=0) labels = comm.bcast(labels, root=0) if rank == 0: time2 = time.time() From 38db4e3bfcadcf86d3b77e5475b517035a5b05a0 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 26 Aug 2016 16:43:58 -0700 Subject: [PATCH 14/59] minor bug fix in fcma example --- examples/fcma/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fcma/run.py b/examples/fcma/run.py index 0efe3b63c..ad5b5465d 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -12,5 +12,5 @@ num_subjs = 18 vs = VoxelSelector(raw_data, epochs_per_subj, labels, num_subjs) results = vs.run() - if rank==0: + if MPI.COMM_WORLD.Get_rank()==0: print(results[0:100]) From 3935570fc7c87a6358a1b198f7e0a42ec4d19685 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Sat, 27 Aug 2016 21:35:58 -0700 Subject: [PATCH 15/59] formalize comments --- brainiak/fcma/fcma.py | 95 +++++++++++++++++++++++----------- examples/fcma/requirements.txt | 1 + examples/fcma/run.py | 26 ++++++++-- 3 files changed, 89 insertions(+), 33 deletions(-) create mode 100644 examples/fcma/requirements.txt diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index cb6c5d28d..e98f362e0 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -37,25 +37,21 @@ from mpi4py import MPI from scipy.stats.mstats import zscore import ctypes -from ctypes.util import find_library from ctypes import byref, c_char, c_float, c_int from sklearn import cross_validation from sklearn import svm from . import fcma_extension - -#__all__ = ['readActivityData', 'separateEpochs', 'VoxelSelector'] - WORKTAG = 0 TERMINATETAG = 1 def readActivityData(dir, file_extension, mask_file): - """ read in data in NIfTI format and apply the spacial mask to them + """ read data in NIfTI format and apply the spatial mask to them - :param dir: the path to all input file + :param dir: the path to all subject files :param file_extension: the file extension, usually nii.gz or nii :param mask_file: the absolute path of the mask file, we apply the mask right after reading a file for saving memory - :return: activity_data: array of matrices in (nTRs, nVoxels) shape (numpy array), + :return: activity_data: list of matrices (numpy array) in (nTRs, nVoxels) shape, len(activity_data) equals the number of subjects """ time1 = time.time() @@ -86,13 +82,15 @@ def readActivityData(dir, file_extension, mask_file): return activity_data def separateEpochs(activity_data, epoch_list): - """ keep data in epochs of interest specified in epoch_map and z-score them for computing correlation + """ separate data into epochs of interest specified in epoch_list and z-score them for computing correlation - :param activity_data: array of matrices in (nTRs, nVoxels) shape, + :param activity_data: list of matrices in (nTRs, nVoxels) shape, consisting of the masked activity data of all subjects - :param epoch_list: array of (subject id, starting TR, ending TR) tuple - :return: raw_data: array of matrices in (epoch length, nVoxels) shape, + :param epoch_list: list of cubes (numpy array) in shape (condition, nEpochs, nTRs) + assuming all subjects have the same number of epochs + :return: raw_data: list of matrices (numpy array) in (epoch length, nVoxels) shape, len(raw_data) equals the number of epochs + labels: list of the condition labels of the epochs, the length of labels equals the number of epochs """ time1 = time.time() raw_data = [] @@ -116,6 +114,16 @@ def separateEpochs(activity_data, epoch_list): return raw_data, labels def prepareData(data_dir, extension, mask_file, epoch_file): + """ read the data in and generate epochs of interests, then broadcast to all workers + + :param data_dir: the path to all subject files + :param extension: the file extension, usually nii.gz or nii + :param mask_file: the absolute path of the mask file + :param epoch_file: the absolute path of the epoch file + :return: raw_data: list of matrices (numpy array) in (epoch length, nVoxels) shape, + len(raw_data) equals the number of epochs + labels: list of the condition labels of the epochs, the length of labels equals the number of epochs + """ comm = MPI.COMM_WORLD rank = comm.Get_rank() labels = [] @@ -148,10 +156,9 @@ class VoxelSelector: which essentially is a nVoxels x nTRs matrix but is serialized as one dimensional vector Assumption: 1. all activity data contains the same number of voxels 2. the activity data has been z-scored, ready to compute correlation as matrix multiplication - 3. the epochs belonging to the same subject are continuous in the array - 4. all subjects have the same number of epochs - 5. voxel selection is always done in the self-correlation, i.e. raw_data correlate with themselves - 6. the classifier used is SVM with linear kernel + 3. all subjects have the same number of epochs + 4. voxel selection is always done in the auto-correlation, i.e. raw_data correlate with themselves + 5. the classifier used is SVM with linear kernel epochs_per_subj : int The number of epochs of each subject @@ -176,9 +183,6 @@ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100) self.labels = labels self.num_folds = num_folds self.voxel_unit = voxel_unit - #libblas_so_name = find_library('blas') if find_library('blas') else find_library('mkl_rt') - #if libblas_so_name is None: - # raise RuntimeError("No Blas library is found in the system") # use the following code since find_library doesn't work for libraries in customized directory in Linux if sys.platform=='darwin': extension = '.dylib' @@ -197,13 +201,15 @@ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100) raise RuntimeError("Zero processed voxels") def run(self): - """run voxel selection in master-worker model + """ run correlation-based voxel selection in master-worker model + Sort the voxels based on the cross-validation accuracy of their correlation vectors + :return: results: list of tuple (voxel_id, accuracy) in accuracy descending order """ rank = MPI.COMM_WORLD.Get_rank() if rank == 0: results = self.master() - # TODO: sort results + # Sort the voxels based on the cross-validation accuracy of their correlation vectors results.sort(key=lambda tup: tup[1], reverse=True) else: self.worker() @@ -211,9 +217,9 @@ def run(self): return results def master(self): - """ - return value - results: an array of tuple (voxel id, accuracy), the length of array equals the number of voxels + """ master node's operation, assigning tasks and collecting results + + :return: results: list of tuple (voxel_id, accuracy), the length of array equals the number of voxels """ results = [] comm = MPI.COMM_WORLD @@ -254,6 +260,10 @@ def master(self): return results def worker(self): + """ worker node's operation, receiving tasks from the master to process and sending the result back + + :return: none + """ comm = MPI.COMM_WORLD status = MPI.Status() while 1: @@ -262,6 +272,11 @@ def worker(self): comm.send(self.voxelScoring(task), dest=0) def correlationComputation(self, task): + """ use BLAS API to do correlation computation (matrix multiplication) + + :param task: a tuple (start_voxel_id, num_assigned_voxels) depicting the voxels assigned to compute + :return: corr: the correlation values + """ s = task[0] e = s+task[1] nEpochs = len(self.raw_data) @@ -285,6 +300,12 @@ def correlationComputation(self, task): return corr def correlationNormalization(self, corr): + """ within-subject normalization + + this method uses scipy.zscore to normalize the data, but is much slower than its C++ counterpart + :param corr: the raw correlation values + :return: corr: the normalized correlation values + """ (sv, e, av) = corr.shape for i in range(sv): start = 0 @@ -299,6 +320,12 @@ def correlationNormalization(self, corr): return corr def crossValidation(self, task, corr): + """ voxelwise cross validation based on correlation vectors + + :param task: a tuple (start_voxel_id, num_assigned_voxels) depicting the voxels assigned to compute + :param corr: the normalized correlation values + :return: results: list of tuple (voxel_id, accuracy), the length of array equals the number of assigned voxels + """ (sv, e, av) = corr.shape kernel_matrix = np.zeros((e, e), np.float32, order='C') results = [] @@ -310,31 +337,41 @@ def crossValidation(self, task, corr): n2 = c_int(self.num_voxels) one = c_float(1.0) zero = c_float(0.0) - self.blas_library.cblas_ssyrk(row_order, lower, no_trans, n1, n2, one, corr[i,:,:].ctypes.data_as(ctypes.c_void_p), + self.blas_library.cblas_ssyrk(row_order, lower, no_trans, n1, n2, one, + corr[i,:,:].ctypes.data_as(ctypes.c_void_p), n2, zero, kernel_matrix.ctypes.data_as(ctypes.c_void_p), n1) kernel_matrix *= .001 for j in range(kernel_matrix.shape[0]): for k in range(j): kernel_matrix[k,j] = kernel_matrix[j,k] + # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) - skf = cross_validation.StratifiedKFold(self.labels, n_folds=self.num_folds, shuffle=False) # no shuffling in cv + # no shuffling in cv + skf = cross_validation.StratifiedKFold(self.labels, n_folds=self.num_folds, shuffle=False) scores = cross_validation.cross_val_score(clf, kernel_matrix, self.labels, cv=skf) results.append((i+task[0], scores.mean())) return results def voxelScoring(self, task): - """Voxel selection worker + """ voxel selection processing done in the worker node + Take the task in, do analysis on voxels specified by the task (voxel id, number of voxels) + It is a three-stage pipeline consisting of: + 1. correlation computation + 2. within-subject normalization + 3. voxelwise cross validaion + :param task: a tuple (start_voxel_id, num_assigned_voxels) depicting the voxels assigned to compute + :return: results: list of tuple (voxel_id, accuracy), the length of array equals the number of assigned voxels """ time1 = time.time() - # 1. correlation computation + # correlation computation corr = self.correlationComputation(task) # corr is a 3D array in row major, # in (selected_voxels, epochs, all_voxels) shape # corr[i,e,s+j] = corr[j,e,s+i] - # 2. normalization + # normalization #corr = self.correlationNormalization(corr) # in-place z-score, the result is still in corr fcma_extension.normalization(corr, self.epochs_per_subj) - # 3. cross validation + # cross validation results = self.crossValidation(task, corr) time2 = time.time() print('task:', task[0]%self.num_voxels, time2-time1) diff --git a/examples/fcma/requirements.txt b/examples/fcma/requirements.txt new file mode 100644 index 000000000..66c5ba046 --- /dev/null +++ b/examples/fcma/requirements.txt @@ -0,0 +1 @@ +mpi4py diff --git a/examples/fcma/run.py b/examples/fcma/run.py index ad5b5465d..a916a5580 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -1,15 +1,33 @@ +# 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. + from brainiak.fcma.fcma import * import sys - -# mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz data/fs_epoch_labels.npy +from mpi4py import MPI +""" +example running command: +mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz + data/fs_epoch_labels.npy 12 18 +""" if __name__ == '__main__': data_dir = sys.argv[1] extension = sys.argv[2] mask_file = sys.argv[3] epoch_file = sys.argv[4] raw_data, labels = prepareData(data_dir, extension, mask_file, epoch_file) - epochs_per_subj = 12 - num_subjs = 18 + epochs_per_subj = int(sys.argv[5]) + num_subjs = int(sys.argv[6]) vs = VoxelSelector(raw_data, epochs_per_subj, labels, num_subjs) results = vs.run() if MPI.COMM_WORLD.Get_rank()==0: From aab70dfde2bdcc82310d29f5ad6b9c2aa934b918 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 29 Aug 2016 10:45:16 -0700 Subject: [PATCH 16/59] impose scikit-learn svm cross validation to use all CPUs --- brainiak/fcma/fcma.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index e98f362e0..3d6fd6c7c 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -96,7 +96,7 @@ def separateEpochs(activity_data, epoch_list): raw_data = [] labels = [] for sid in range(len(epoch_list)): - epoch = epoch_list[sid] # epoch is a numpy array in shape (condition, nEpochs, nTRs) + epoch = epoch_list[sid] for cond in range(epoch.shape[0]): sub_epoch = epoch[cond,:,:] for eid in range(epoch.shape[1]): @@ -348,7 +348,7 @@ def crossValidation(self, task, corr): clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) # no shuffling in cv skf = cross_validation.StratifiedKFold(self.labels, n_folds=self.num_folds, shuffle=False) - scores = cross_validation.cross_val_score(clf, kernel_matrix, self.labels, cv=skf) + scores = cross_validation.cross_val_score(clf, kernel_matrix, self.labels, cv=skf, n_jobs=-1) results.append((i+task[0], scores.mean())) return results @@ -374,6 +374,6 @@ def voxelScoring(self, task): # cross validation results = self.crossValidation(task, corr) time2 = time.time() - print('task:', task[0]%self.num_voxels, time2-time1) + print('task:', int(task[0]%self.num_voxels), time2-time1) sys.stdout.flush() return results From 67d97d120d1e4bdd21c3ed5f98cca83b0a596d92 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 29 Aug 2016 11:30:27 -0700 Subject: [PATCH 17/59] tiny bug fixed in the worker output --- brainiak/fcma/fcma.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 3d6fd6c7c..57f60872f 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -374,6 +374,6 @@ def voxelScoring(self, task): # cross validation results = self.crossValidation(task, corr) time2 = time.time() - print('task:', int(task[0]%self.num_voxels), time2-time1) + print('task:', int(task[0]/self.voxel_unit), time2-time1) sys.stdout.flush() return results From d0af529400b203412815ff563b1209cef3ac6fd6 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 29 Aug 2016 19:46:47 -0700 Subject: [PATCH 18/59] change to column major in preparation for using BLAS APIs under SciPy Cython wrapper --- brainiak/fcma/fcma.py | 57 ++++++++++++++---------- brainiak/fcma/src/fcma_extension.cc | 67 +++++++++++++++-------------- 2 files changed, 69 insertions(+), 55 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 57f60872f..bd8fcd1b6 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -67,7 +67,7 @@ def readActivityData(dir, file_extension, mask_file): img = nib.load(os.path.join(dir, f)) data = img.get_data() (d1, d2, d3, d4) = data.shape - masked_data = np.zeros([d4, count], np.float32, order='C') + masked_data = np.zeros([d4, count], np.float32, order='F') count1=0 for index in np.ndindex(mask.shape): if mask[index]!=0: @@ -102,7 +102,7 @@ def separateEpochs(activity_data, epoch_list): for eid in range(epoch.shape[1]): r = np.sum(sub_epoch[eid,:]) if r > 0: # there is an epoch in this condition - mat = activity_data[sid][sub_epoch[eid,:]==1,:] + mat = np.asfortranarray(activity_data[sid][sub_epoch[eid,:]==1,:]) mat = zscore(mat, axis=0, ddof=0) mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero mat = mat / math.sqrt(r) @@ -275,27 +275,33 @@ def correlationComputation(self, task): """ use BLAS API to do correlation computation (matrix multiplication) :param task: a tuple (start_voxel_id, num_assigned_voxels) depicting the voxels assigned to compute - :return: corr: the correlation values + :return: corr: the correlation values, column major, in shape (num_voxels, num_epochs, num_selected_voxels) """ s = task[0] e = s+task[1] nEpochs = len(self.raw_data) - corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') + #corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='F') + corr = np.zeros((self.num_voxels, nEpochs, task[1]), np.float32, order='F') count=0 for mat in self.raw_data: - row_order = c_int(101) + col_order = c_int(102) no_trans = c_int(111) trans = c_int(112) n1 = c_int(task[1]) n2 = c_int(self.num_voxels) n3 = c_int(mat.shape[0]) n4 = c_int(self.num_voxels*nEpochs) + #print(mat.flags['F_CONTIGUOUS'], mat.shape, n1, n2, n3, n4, s, e, mat[:,s:e].shape) one = c_float(1.0) zero = c_float(0.0) - self.blas_library.cblas_sgemm(row_order, trans, no_trans, n1, n2, n3, one, - mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, - mat.ctypes.data_as(ctypes.c_void_p), n2, - zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) + #self.blas_library.cblas_sgemm(row_order, trans, no_trans, n1, n2, n3, one, + # mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, + # mat.ctypes.data_as(ctypes.c_void_p), n2, + # zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) + self.blas_library.cblas_sgemm(col_order, trans, no_trans, n2, n1, n3, one, + mat.ctypes.data_as(ctypes.c_void_p), n3, + mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n3, + zero, corr[:,count,0].ctypes.data_as(ctypes.c_void_p), n4) count += 1 return corr @@ -306,15 +312,15 @@ def correlationNormalization(self, corr): :param corr: the raw correlation values :return: corr: the normalized correlation values """ - (sv, e, av) = corr.shape + (av, e, sv) = corr.shape for i in range(sv): start = 0 while start py_data, int epochsPerSubj) { // the input data is a 3D array in row major, in (selected_voxels, epochs, all_voxels) shape @@ -34,45 +34,46 @@ void within_subject_norm_native(py::array_t> 2, // (int)inBuf.strides[1] >> 2, (int)inBuf.strides[2] >> 2, // (int)inBuf.shape[3]}; - #pragma omp parallel for - for(int v = 0 ; v < nSelectedVoxels*nSubjs ; v++) - { - int s = v % nSubjs; // subject id - int i = v / nSubjs; // voxel id - //float (*mat)[nVoxels] = (float(*)[nVoxels])&(data[i*nEpochs*nVoxels]); - float *mat = data+i*nEpochs*nVoxels; - #pragma simd - for(int j = 0 ; j < nVoxels ; j++) + #pragma omp parallel for + for(int v = 0 ; v < nSelectedVoxels*nSubjs ; v++) { - float mean = 0.0f; - float std_dev = 0.0f; - for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) - { - float num = 1.0f + mat[b*nVoxels+j]; - float den = 1.0f - mat[b*nVoxels+j]; - num = (num <= 0.0f) ? 1e-4 : num; - den = (den <= 0.0f) ? 1e-4 : den; - mat[b*nVoxels+j] = 0.5f * logf(num/den); - mean += mat[b*nVoxels+j]; - std_dev += mat[b*nVoxels+j] * mat[b*nVoxels+j]; - } - mean = mean / (float)epochsPerSubj; - std_dev = std_dev / (float)epochsPerSubj - mean*mean; - float inv_std_dev = (std_dev <= 0.0f) ? 0.0f : 1.0f / sqrt(std_dev); - for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) - { - mat[b*nVoxels+j] = (mat[b*nVoxels+j] - mean) * inv_std_dev; - } + int s = v % nSubjs; // subject id + int i = v / nSubjs; // voxel id + //float (*mat)[nVoxels] = (float(*)[nVoxels])&(data[i*nEpochs*nVoxels]); + float *mat = data+i*nEpochs*nVoxels; + #pragma simd + for(int j = 0 ; j < nVoxels ; j++) + { + float mean = 0.0f; + float std_dev = 0.0f; + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + float num = 1.0f + mat[b*nVoxels+j]; + float den = 1.0f - mat[b*nVoxels+j]; + num = (num <= 0.0f) ? 1e-4 : num; + den = (den <= 0.0f) ? 1e-4 : den; + mat[b*nVoxels+j] = 0.5f * logf(num/den); + mean += mat[b*nVoxels+j]; + std_dev += mat[b*nVoxels+j] * mat[b*nVoxels+j]; + } + mean = mean / (float)epochsPerSubj; + std_dev = std_dev / (float)epochsPerSubj - mean*mean; + float inv_std_dev = (std_dev <= 0.0f) ? 0.0f : 1.0f / sqrt(std_dev); + for(int b = s*epochsPerSubj; b < (s+1)*epochsPerSubj; b++) + { + mat[b*nVoxels+j] = (mat[b*nVoxels+j] - mean) * inv_std_dev; + } + } } - } - return; + return; } PYBIND11_PLUGIN(fcma_extension) { From 55565af654ce4bc7d978312a278745f4f2cfb085 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 30 Aug 2016 08:24:54 -0700 Subject: [PATCH 19/59] change the input data to row-major --- brainiak/fcma/fcma.py | 40 +++++++++++++++++------------ brainiak/fcma/src/fcma_extension.cc | 15 +++++------ 2 files changed, 30 insertions(+), 25 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index bd8fcd1b6..73ebdf04b 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -41,6 +41,9 @@ from sklearn import cross_validation from sklearn import svm from . import fcma_extension +#import pyximport +#pyximport.install() +#from . import cython_blas WORKTAG = 0 TERMINATETAG = 1 @@ -67,7 +70,7 @@ def readActivityData(dir, file_extension, mask_file): img = nib.load(os.path.join(dir, f)) data = img.get_data() (d1, d2, d3, d4) = data.shape - masked_data = np.zeros([d4, count], np.float32, order='F') + masked_data = np.zeros([d4, count], np.float32, order='C') count1=0 for index in np.ndindex(mask.shape): if mask[index]!=0: @@ -102,7 +105,8 @@ def separateEpochs(activity_data, epoch_list): for eid in range(epoch.shape[1]): r = np.sum(sub_epoch[eid,:]) if r > 0: # there is an epoch in this condition - mat = np.asfortranarray(activity_data[sid][sub_epoch[eid,:]==1,:]) + # mat is row-major no matter what order acitvity_data[sid] is + mat = activity_data[sid][sub_epoch[eid,:]==1,:] mat = zscore(mat, axis=0, ddof=0) mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero mat = mat / math.sqrt(r) @@ -281,7 +285,7 @@ def correlationComputation(self, task): e = s+task[1] nEpochs = len(self.raw_data) #corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='F') - corr = np.zeros((self.num_voxels, nEpochs, task[1]), np.float32, order='F') + corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') count=0 for mat in self.raw_data: col_order = c_int(102) @@ -298,10 +302,14 @@ def correlationComputation(self, task): # mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, # mat.ctypes.data_as(ctypes.c_void_p), n2, # zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) - self.blas_library.cblas_sgemm(col_order, trans, no_trans, n2, n1, n3, one, - mat.ctypes.data_as(ctypes.c_void_p), n3, - mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n3, - zero, corr[:,count,0].ctypes.data_as(ctypes.c_void_p), n4) + self.blas_library.cblas_sgemm(col_order, no_trans, trans, n2, n1, n3, one, + mat.ctypes.data_as(ctypes.c_void_p), n2, + mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n1, + zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) + #blas.sgemm('T', 'N', n2, n1, n3, one, + # mat.ctypes.data_as(ctypes.c_void_p), n3, + # mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n3, + # zero, corr[:,count,0].ctypes.data_as(ctypes.c_void_p), n4) count += 1 return corr @@ -312,15 +320,15 @@ def correlationNormalization(self, corr): :param corr: the raw correlation values :return: corr: the normalized correlation values """ - (av, e, sv) = corr.shape + (sv, e, av) = corr.shape for i in range(sv): start = 0 while start py_data, int epochsPerSubj) { // the input data is a 3D array in row major, in (selected_voxels, epochs, all_voxels) shape @@ -34,20 +34,17 @@ void within_subject_norm_native(py::array_t> 2, - // (int)inBuf.strides[1] >> 2, (int)inBuf.strides[2] >> 2, - // (int)inBuf.shape[3]}; + #pragma omp parallel for for(int v = 0 ; v < nSelectedVoxels*nSubjs ; v++) { int s = v % nSubjs; // subject id int i = v / nSubjs; // voxel id - //float (*mat)[nVoxels] = (float(*)[nVoxels])&(data[i*nEpochs*nVoxels]); + float *mat = data+i*nEpochs*nVoxels; #pragma simd for(int j = 0 ; j < nVoxels ; j++) @@ -86,4 +83,4 @@ PYBIND11_PLUGIN(fcma_extension) { m.attr("__version__") = py::str("dev"); #endif return m.ptr(); -} \ No newline at end of file +} From 29f6d476b6dbd35e0e8dfc9c4f4ac4a8d29586d6 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 30 Aug 2016 09:32:02 -0700 Subject: [PATCH 20/59] minor bugs fixed --- brainiak/fcma/fcma.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 73ebdf04b..9c439e51c 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -288,6 +288,7 @@ def correlationComputation(self, task): corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') count=0 for mat in self.raw_data: + row_order = c_int(101) col_order = c_int(102) no_trans = c_int(111) trans = c_int(112) @@ -304,7 +305,7 @@ def correlationComputation(self, task): # zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) self.blas_library.cblas_sgemm(col_order, no_trans, trans, n2, n1, n3, one, mat.ctypes.data_as(ctypes.c_void_p), n2, - mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n1, + mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) #blas.sgemm('T', 'N', n2, n1, n3, one, # mat.ctypes.data_as(ctypes.c_void_p), n3, @@ -357,7 +358,7 @@ def crossValidation(self, task, corr): kernel_matrix *= .001 for j in range(kernel_matrix.shape[0]): for k in range(j): - kernel_matrix[k,j] = kernel_matrix[j,k] + kernel_matrix[j,k] = kernel_matrix[k,j] # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) # no shuffling in cv From 044d34c783a8f7f185608c1cff1a5d01d27bd3d3 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 30 Aug 2016 14:56:24 -0700 Subject: [PATCH 21/59] use blas wrapper by scipy cython --- brainiak/fcma/__init__.py | 3 ++ brainiak/fcma/cython_blas.pyx | 67 +++++++++++++++++++++++++++++++++++ brainiak/fcma/fcma.py | 64 +++++---------------------------- setup.py | 3 ++ 4 files changed, 81 insertions(+), 56 deletions(-) create mode 100644 brainiak/fcma/cython_blas.pyx diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py index 5f9a9900f..3880ca640 100644 --- a/brainiak/fcma/__init__.py +++ b/brainiak/fcma/__init__.py @@ -12,3 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. """Full correlation matrix analysis""" + +import pyximport +pyximport.install() diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx new file mode 100644 index 000000000..1c4d8f7f1 --- /dev/null +++ b/brainiak/fcma/cython_blas.pyx @@ -0,0 +1,67 @@ +# 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. + +cimport scipy.linalg.cython_blas as blas + +import numpy as np + +def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a, py_lda, + py_start_voxel, py_ldb, py_beta, py_c, py_ldc, py_start_epoch): + 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] C + C = py_c + blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0,0], &lda, + &A[0,py_start_voxel], &ldb, &beta, &C[0,py_start_epoch,0], &ldc) + +def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_start_voxel, py_lda, + py_beta, py_c, py_ldc): + cdef bytes by_uplo=py_uplo.encode() + cdef bytes by_trans=py_trans.encode() + cdef char* uplo = by_uplo + cdef char* trans = by_trans + cdef int N, K, lda, ldc + N = py_n + K = py_k + lda = py_lda + ldc = py_ldc + cdef float alpha, beta + alpha = py_alpha + beta = py_beta + cdef float[:, :, ::1] A + A = py_a + cdef float[:, ::1] C + C = py_c + blas.ssyrk(uplo, trans, &N, &K, &alpha, &A[py_start_voxel,0,0], &lda, + &beta, &C[0,0], &ldc) + # complete the upper triangle of the kernel matrix + # shrink the values for getting more stable alpha values in SVM training iteration + py_c *= .001 + for j in range(py_c.shape[0]): + for k in range(j): + py_c[j,k] = py_c[k,j] diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 9c439e51c..8392a2c37 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -36,14 +36,10 @@ import os, math, sys, time from mpi4py import MPI from scipy.stats.mstats import zscore -import ctypes -from ctypes import byref, c_char, c_float, c_int from sklearn import cross_validation from sklearn import svm from . import fcma_extension -#import pyximport -#pyximport.install() -#from . import cython_blas +from . import cython_blas as blas WORKTAG = 0 TERMINATETAG = 1 @@ -187,20 +183,6 @@ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100) self.labels = labels self.num_folds = num_folds self.voxel_unit = voxel_unit - # use the following code since find_library doesn't work for libraries in customized directory in Linux - if sys.platform=='darwin': - extension = '.dylib' - elif sys.platform=='linux': - extension = '.so' - else: - raise RuntimeError("Unsupported operating system") - try: - self.blas_library = ctypes.cdll.LoadLibrary('libblas'+extension) - except: - try: - self.blas_library = ctypes.cdll.LoadLibrary('libmkl_rt'+extension) - except: - raise RuntimeError("No Blas library is found in the system") if self.num_voxels == 0: raise RuntimeError("Zero processed voxels") @@ -288,29 +270,11 @@ def correlationComputation(self, task): corr = np.zeros((task[1], nEpochs, self.num_voxels), np.float32, order='C') count=0 for mat in self.raw_data: - row_order = c_int(101) - col_order = c_int(102) - no_trans = c_int(111) - trans = c_int(112) - n1 = c_int(task[1]) - n2 = c_int(self.num_voxels) - n3 = c_int(mat.shape[0]) - n4 = c_int(self.num_voxels*nEpochs) - #print(mat.flags['F_CONTIGUOUS'], mat.shape, n1, n2, n3, n4, s, e, mat[:,s:e].shape) - one = c_float(1.0) - zero = c_float(0.0) - #self.blas_library.cblas_sgemm(row_order, trans, no_trans, n1, n2, n3, one, - # mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, - # mat.ctypes.data_as(ctypes.c_void_p), n2, - # zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) - self.blas_library.cblas_sgemm(col_order, no_trans, trans, n2, n1, n3, one, - mat.ctypes.data_as(ctypes.c_void_p), n2, - mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n2, - zero, corr[0,count,:].ctypes.data_as(ctypes.c_void_p), n4) - #blas.sgemm('T', 'N', n2, n1, n3, one, - # mat.ctypes.data_as(ctypes.c_void_p), n3, - # mat[:,s:e].ctypes.data_as(ctypes.c_void_p), n3, - # zero, corr[:,count,0].ctypes.data_as(ctypes.c_void_p), n4) + 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) count += 1 return corr @@ -345,20 +309,8 @@ def crossValidation(self, task, corr): kernel_matrix = np.zeros((e, e), np.float32, order='C') results = [] for i in range(sv): - col_order = c_int(102) - trans = c_int(112) - lower = c_int(122) - n1 = c_int(e) - n2 = c_int(self.num_voxels) - one = c_float(1.0) - zero = c_float(0.0) - self.blas_library.cblas_ssyrk(col_order, lower, trans, n1, n2, one, - corr[i,:,:].ctypes.data_as(ctypes.c_void_p), - n2, zero, kernel_matrix.ctypes.data_as(ctypes.c_void_p), n1) - kernel_matrix *= .001 - for j in range(kernel_matrix.shape[0]): - for k in range(j): - kernel_matrix[j,k] = kernel_matrix[k,j] + blas.compute_kernel_matrix('L', 'T', e, self.num_voxels, 1.0, corr, i, self.num_voxels, + 0.0, kernel_matrix, e) # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) # no shuffling in cv diff --git a/setup.py b/setup.py index c233c6383..a2cef3222 100644 --- a/setup.py +++ b/setup.py @@ -169,5 +169,8 @@ def cmake_compiling(self, ext): ext_modules=ext_modules, cmdclass={'build_ext': BuildExt}, packages=find_packages(), + package_data = { + 'brainiak.fcma': ['*.pyx'], + }, zip_safe=False, ) From 849a0616c87aa619f6d24ce149843c617af422d9 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 30 Aug 2016 15:28:36 -0700 Subject: [PATCH 22/59] format output --- brainiak/fcma/fcma.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 8392a2c37..6f5c37053 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -76,7 +76,7 @@ def readActivityData(dir, file_extension, mask_file): print(f, masked_data.shape) sys.stdout.flush() time2 = time.time() - print('data reading done, takes', time2-time1, 's') + print('data reading done, takes', round(time2-time1, 2), 's') sys.stdout.flush() return activity_data @@ -109,7 +109,7 @@ def separateEpochs(activity_data, epoch_list): raw_data.append(mat) labels.append(cond) time2 = time.time() - print('epoch separation done, takes', time2-time1, 's') + print('epoch separation done, takes', round(time2-time1, 2), 's') sys.stdout.flush() return raw_data, labels @@ -142,7 +142,7 @@ def prepareData(data_dir, extension, mask_file, epoch_file): labels = comm.bcast(labels, root=0) if rank == 0: time2 = time.time() - print('data broadcasting done, takes', time2-time1, 's') + print('data broadcasting done, takes', round(time2-time1, 2), 's') return raw_data, labels class VoxelSelector: @@ -335,19 +335,13 @@ def voxelScoring(self, task): corr = self.correlationComputation(task) # corr is a 3D array in row major, # in (selected_voxels, epochs, all_voxels) shape # corr[i,e,s+j] = corr[j,e,s+i] - time3 = time.time() - print('corr comp', time3-time1) - sys.stdout.flush() # normalization #corr = self.correlationNormalization(corr) # in-place z-score, the result is still in corr fcma_extension.normalization(corr, self.epochs_per_subj) - time4 = time.time() - print('norm', time4-time3) - sys.stdout.flush() + # cross validation results = self.crossValidation(task, corr) time2 = time.time() - print('cv', time2-time4) - print('task:', int(task[0]/self.voxel_unit), time2-time1) + print('task:', int(task[0]/self.voxel_unit), round(time2-time1, 2), 's') sys.stdout.flush() return results From 2bdd56969f1765d68a15bf031bbac6fe336b3220 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 31 Aug 2016 15:27:46 -0700 Subject: [PATCH 23/59] add in classifier flexibility, use regular compiling for fcma_extension, format the code --- .gitmodules | 3 - brainiak/fcma/cython_blas.pyx | 117 ++++++++- brainiak/fcma/fcma.py | 393 +++++++++++++++++++--------- brainiak/fcma/pybind11 | 1 - brainiak/fcma/src/fcma_extension.cc | 13 +- examples/fcma/requirements.txt | 1 + examples/fcma/run.py | 6 +- setup.py | 11 +- 8 files changed, 409 insertions(+), 136 deletions(-) delete mode 100644 .gitmodules delete mode 160000 brainiak/fcma/pybind11 diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index d19b82e8b..000000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "brainiak/fcma/pybind11"] - path = brainiak/fcma/pybind11 - url = https://github.com/pybind/pybind11.git diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 1c4d8f7f1..649061642 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -14,10 +14,60 @@ cimport scipy.linalg.cython_blas as blas -import numpy as np - def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a, py_lda, py_start_voxel, py_ldb, py_beta, py_c, py_ldc, py_start_epoch): + """ use blas API wrapped by scipy.linalg.cython_blas to compute 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 + + The resulting matrix in shape [num_assigned_voxels, num_voxels] + is stored in an alternate way to make sure that + the correlation vectors of the same voxel stored continuously + + Parameters + ---------- + 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 + in our case, is num_voxels + py_n: int + the column of the resulting matrix C + in our case, is num_assigned_voxels + 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 first matrix A + py_a: 2D array in shape [epoch_length, num_voxels] in our case + the activity data of an epoch + py_lda: int + the stride of the first matrix A + py_start_voxel: int + the starting voxel of assigned voxels + used to locate the second matrix B + py_ldb: int + the stride of the second matrix B + py_beta: float + the weight applied to the resulting matrix C + py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + place to store the resulting correlation values + py_ldc: int + the stride of the resulting matrix + in our case, num_voxels*num_epochs + py_start_epoch: int + the epoch over which the correlation is computed + + Returns + ------- + py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] + write the resulting correlation values in an alternate way + for the processing epoch + """ 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 @@ -41,6 +91,56 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_start_voxel, py_lda, py_beta, py_c, py_ldc): + """ use blas API wrapped by scipy.linalg.cython_blas to compute kernel matrix of SVM + + 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 + + In SVM with linear kernel, the distance of two samples + is essentially the dot product of them. + Therefore, the kernel matrix can be obtained by matrix multiplication. + Since the kernel matrix is symmetric, ssyrk is used, + the other half of the matrix is assigned later. + In our case, the dimension of samples is much larger than + the number samples, so we proportionally shrink the values of the kernel matrix + for getting more robust alpha values in SVM iteration. + + Parameters + ---------- + py_uplo: str + getting the upper or lower triangle of the matrix + py_trans: str + do transpose or not for the input matrix A + py_n: int + the row and column of the resulting matrix C + in our case, is num_epochs + 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_voxels + py_alpha: float + the weight applied to the input matrix A + py_a: 3D array in shape [num_assigned_voxels, num_epochs, num_voxels] in our case + the normalized correlation values of a voxel + py_start_voxel: int + the processed voxel + used to locate the input matrix A + py_lda: int + the stride of the input matrix A + py_beta: float + the weight applied to the resulting matrix C + py_c: 2D array in shape [num_epochs, num_epochs] + place to store the resulting kernel matrix + py_ldc: int + the stride of the resulting matrix + + Returns + ------- + py_c: 2D array in shape [num_epochs, num_epochs] + write the resulting kernel_matrix + for the processing voxel + """ cdef bytes by_uplo=py_uplo.encode() cdef bytes by_trans=py_trans.encode() cdef char* uplo = by_uplo @@ -59,9 +159,14 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_star C = py_c blas.ssyrk(uplo, trans, &N, &K, &alpha, &A[py_start_voxel,0,0], &lda, &beta, &C[0,0], &ldc) - # complete the upper triangle of the kernel matrix # shrink the values for getting more stable alpha values in SVM training iteration py_c *= .001 - for j in range(py_c.shape[0]): - for k in range(j): - py_c[j,k] = py_c[k,j] + # complete the other half of the kernel matrix + if (py_uplo=='L'): + for j in range(py_c.shape[0]): + for k in range(j): + py_c[j,k] = py_c[k,j] + else: + for j in range(py_c.shape[0]): + for k in range(j): + py_c[k,j] = py_c[j,k] diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 6f5c37053..5975e35aa 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -22,9 +22,11 @@ .. [Wang2015-2] "Full correlation matrix analysis of fMRI data on Intel® Xeon Phi™ coprocessors", - Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, Kai Li, - Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, Theodore L. Willke. - In Proceedings of the International Conference for High Performance Computing, + Yida Wang, Michael J. Anderson, Jonathan D. Cohen, Alexander Heinecke, + Kai Li, Nadathur Satish, Narayanan Sundaram, Nicholas B. Turk-Browne, + Theodore L. Willke. + In Proceedings of the International Conference for + High Performance Computing, Networking, Storage and Analysis. 2015. """ @@ -33,63 +35,92 @@ import numpy as np import nibabel as nib -import os, math, sys, time +import os +import math +import sys +import time from mpi4py import MPI from scipy.stats.mstats import zscore from sklearn import cross_validation -from sklearn import svm +import sklearn from . import fcma_extension from . import cython_blas as blas WORKTAG = 0 TERMINATETAG = 1 + def readActivityData(dir, file_extension, mask_file): """ read data in NIfTI format and apply the spatial mask to them - :param dir: the path to all subject files - :param file_extension: the file extension, usually nii.gz or nii - :param mask_file: the absolute path of the mask file, we apply the mask right after reading a file for saving memory - :return: activity_data: list of matrices (numpy array) in (nTRs, nVoxels) shape, - len(activity_data) equals the number of subjects + Parameters + ---------- + dir: str + the path to all subject files + file_extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, we apply the mask right after + reading a file for saving memory + + Returns + ------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats + len(activity_data) equals the number of subjects """ time1 = time.time() mask_img = nib.load(mask_file) mask = mask_img.get_data() count = 0 for index in np.ndindex(mask.shape): - if mask[index]!=0: - count+=1 - files = [f for f in sorted(os.listdir(dir)) if os.path.isfile(os.path.join(dir, f)) and f.endswith(file_extension)] + if mask[index] != 0: + count += 1 + files = [f for f in sorted(os.listdir(dir)) + if os.path.isfile(os.path.join(dir, f)) + and f.endswith(file_extension)] activity_data = [] for f in files: img = nib.load(os.path.join(dir, f)) data = img.get_data() (d1, d2, d3, d4) = data.shape masked_data = np.zeros([d4, count], np.float32, order='C') - count1=0 + count1 = 0 for index in np.ndindex(mask.shape): - if mask[index]!=0: - masked_data[:,count1] = np.copy(data[index]) - count1+=1 + if mask[index] != 0: + masked_data[:, count1] = np.copy(data[index]) + count1 += 1 activity_data.append(masked_data) print(f, masked_data.shape) sys.stdout.flush() time2 = time.time() - print('data reading done, takes', round(time2-time1, 2), 's') + print('data reading done, takes', round(time2 - time1, 2), 's') sys.stdout.flush() return activity_data + def separateEpochs(activity_data, epoch_list): - """ separate data into epochs of interest specified in epoch_list and z-score them for computing correlation - - :param activity_data: list of matrices in (nTRs, nVoxels) shape, - consisting of the masked activity data of all subjects - :param epoch_list: list of cubes (numpy array) in shape (condition, nEpochs, nTRs) - assuming all subjects have the same number of epochs - :return: raw_data: list of matrices (numpy array) in (epoch length, nVoxels) shape, - len(raw_data) equals the number of epochs - labels: list of the condition labels of the epochs, the length of labels equals the number of epochs + """ separate data into epochs of interest specified in epoch_list + and z-score them for computing correlation + + Parameters + ---------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats of all subjects + epoch_list: list of 3D array in shape [condition, nEpochs, nTRs] + specification of epochs and conditions + assuming all subjects have the same number of epochs + len(epoch_list) equals the number of subjects + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + and z-scored in preparation of correlation computation + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs """ time1 = time.time() raw_data = [] @@ -97,78 +128,98 @@ def separateEpochs(activity_data, epoch_list): for sid in range(len(epoch_list)): epoch = epoch_list[sid] for cond in range(epoch.shape[0]): - sub_epoch = epoch[cond,:,:] + sub_epoch = epoch[cond, :, :] for eid in range(epoch.shape[1]): - r = np.sum(sub_epoch[eid,:]) + r = np.sum(sub_epoch[eid, :]) if r > 0: # there is an epoch in this condition - # mat is row-major no matter what order acitvity_data[sid] is - mat = activity_data[sid][sub_epoch[eid,:]==1,:] + # mat is row-major + # regardless of the order of acitvity_data[sid] + mat = activity_data[sid][sub_epoch[eid, :] == 1, :] mat = zscore(mat, axis=0, ddof=0) - mat = np.nan_to_num(mat) # if zscore fails (standard deviation is zero), set all values to be zero + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) mat = mat / math.sqrt(r) raw_data.append(mat) labels.append(cond) time2 = time.time() - print('epoch separation done, takes', round(time2-time1, 2), 's') + print('epoch separation done, takes', round(time2 - time1, 2), 's') sys.stdout.flush() return raw_data, labels + def prepareData(data_dir, extension, mask_file, epoch_file): - """ read the data in and generate epochs of interests, then broadcast to all workers - - :param data_dir: the path to all subject files - :param extension: the file extension, usually nii.gz or nii - :param mask_file: the absolute path of the mask file - :param epoch_file: the absolute path of the epoch file - :return: raw_data: list of matrices (numpy array) in (epoch length, nVoxels) shape, - len(raw_data) equals the number of epochs - labels: list of the condition labels of the epochs, the length of labels equals the number of epochs + """ read the data in and generate epochs of interests, + then broadcast to all workers + + Parameters + ---------- + data_dir: str + the path to all subject files + extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, + we apply the mask right after reading a file for saving memory + epoch_file: str + the absolute path of the epoch file + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs """ comm = MPI.COMM_WORLD rank = comm.Get_rank() labels = [] raw_data = [] - if rank==0: + if rank == 0: activity_data = readActivityData(data_dir, extension, mask_file) - epoch_list = np.load(epoch_file) # a list of numpy array in shape (condition, nEpochs, nTRs) - raw_data, labels=separateEpochs(activity_data, epoch_list) + # a list of numpy array in shape [condition, nEpochs, nTRs] + epoch_list = np.load(epoch_file) + raw_data, labels = separateEpochs(activity_data, epoch_list) time1 = time.time() raw_data_length = len(raw_data) raw_data_length = comm.bcast(raw_data_length, root=0) - for i in range(raw_data_length): # broadcast the data subject by subject to prevent size overflow - if rank!=0: + # broadcast the data subject by subject to prevent size overflow + for i in range(raw_data_length): + if rank != 0: raw_data.append(None) raw_data[i] = comm.bcast(raw_data[i], root=0) labels = comm.bcast(labels, root=0) if rank == 0: time2 = time.time() - print('data broadcasting done, takes', round(time2-time1, 2), 's') + print('data broadcasting done, takes', round(time2 - time1, 2), 's') return raw_data, labels + class VoxelSelector: """Correlation-based voxel selection component of FCMA Parameters ---------- - raw_data : numpy array of array of float - Each array of float contains the activity data of an epoch of a subject, - which essentially is a nVoxels x nTRs matrix but is serialized as one dimensional vector + raw_data: list of 2D array in shape [epoch length, nVoxels] Assumption: 1. all activity data contains the same number of voxels - 2. the activity data has been z-scored, ready to compute correlation as matrix multiplication + 2. the activity data has been z-scored, + ready to compute correlation as matrix multiplication 3. all subjects have the same number of epochs - 4. voxel selection is always done in the auto-correlation, i.e. raw_data correlate with themselves - 5. the classifier used is SVM with linear kernel + 4. voxel selection is always done in the auto-correlation, + i.e. raw_data correlate with themselves - epochs_per_subj : int + epochs_per_subj: int The number of epochs of each subject num_voxels: int The number of voxels participating in the voxel selection - labels: 1D array of int - Describing the condition labels of the epochs, the length of labels equals the number of epochs, - i.e. the length of raw_data + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs num_folds: int The number of folds to be conducted in the cross validation @@ -176,7 +227,12 @@ class VoxelSelector: voxel_unit: int, default 100 The number of voxel assigned to a worker each time """ - def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100): + def __init__(self, + raw_data, + epochs_per_subj, + labels, + num_folds, + voxel_unit=100): self.raw_data = raw_data self.epochs_per_subj = epochs_per_subj self.num_voxels = raw_data[0].shape[1] @@ -186,54 +242,79 @@ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, voxel_unit=100) if self.num_voxels == 0: raise RuntimeError("Zero processed voxels") - def run(self): + def run(self, clf): """ run correlation-based voxel selection in master-worker model - Sort the voxels based on the cross-validation accuracy of their correlation vectors - :return: results: list of tuple (voxel_id, accuracy) in accuracy descending order + Sort the voxels based on the cross-validation accuracy + of their correlation vectors + + Parameters + ---------- + clf: classification function + the classifier to be used in cross validation + + Returns + ------- + results: list of tuple (voxel_id, accuracy) + the accuracy numbers of all voxels, in accuracy descending order + the length of array equals the number of voxels """ rank = MPI.COMM_WORLD.Get_rank() if rank == 0: results = self.master() - # Sort the voxels based on the cross-validation accuracy of their correlation vectors + # Sort the voxels results.sort(key=lambda tup: tup[1], reverse=True) else: - self.worker() + self.worker(clf) results = [] return results def master(self): """ master node's operation, assigning tasks and collecting results - :return: results: list of tuple (voxel_id, accuracy), the length of array equals the number of voxels + Parameters + ---------- + None + + Returns + ------- + results: list of tuple (voxel_id, accuracy) + the accuracy numbers of all voxels, in accuracy descending order + the length of array equals the number of voxels """ results = [] comm = MPI.COMM_WORLD size = comm.Get_size() - sending_voxels = self.voxel_unit if self.voxel_unit py_data, int epochsPerSubj) { - // the input data is a 3D array in row major, in (selected_voxels, epochs, all_voxels) shape +//This function calculates the within-subject normalization of the input correlation values +// +//Input parameters are: +//py_data: the correlation value array, in shape [num_voxels, num_epochs, num_selected_voxels] +// containing the raw correlation values +//epochsPerSubj: the number pf epochs per subject +// assuming all subjects have the same number of epochs +// +//Ouput parameters are: +//py_data: the correlation value array, in shape [num_voxels, num_epochs, num_selected_voxels] +// containing the normalized correlation values +// py::buffer_info buf = py_data.request(); float* data = (float*)buf.ptr; // sanity check diff --git a/examples/fcma/requirements.txt b/examples/fcma/requirements.txt index 66c5ba046..8808a066e 100644 --- a/examples/fcma/requirements.txt +++ b/examples/fcma/requirements.txt @@ -1 +1,2 @@ mpi4py +sklearn diff --git a/examples/fcma/run.py b/examples/fcma/run.py index a916a5580..1911b4051 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -13,6 +13,7 @@ # limitations under the License. from brainiak.fcma.fcma import * +from sklearn import svm import sys from mpi4py import MPI """ @@ -29,6 +30,9 @@ epochs_per_subj = int(sys.argv[5]) num_subjs = int(sys.argv[6]) vs = VoxelSelector(raw_data, epochs_per_subj, labels, num_subjs) - results = vs.run() + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=10 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) + results = vs.run(clf) if MPI.COMM_WORLD.Get_rank()==0: print(results[0:100]) diff --git a/setup.py b/setup.py index a2cef3222..c2ccf6bf9 100644 --- a/setup.py +++ b/setup.py @@ -48,7 +48,16 @@ def __init__(self, name, sourcedir=''): get_pybind_include(user=True) ], ), - CMakeExtension('brainiak.fcma.fcma_extension', 'brainiak/fcma/'), + #CMakeExtension('brainiak.fcma.fcma_extension', 'brainiak/fcma/'), + Extension( + 'brainiak.fcma.fcma_extension', + ['brainiak/fcma/src/fcma_extension.cc'], + include_dirs=[ + # Path to pybind11 headers + get_pybind_include(), + get_pybind_include(user=True) + ], + ), ] From 5d7ea8be4e2834bd529add018935c7a35acff219 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 31 Aug 2016 18:25:40 -0700 Subject: [PATCH 24/59] add fcma test, modify run-tests.sh for testing out MPI processes and collecting the right coverage --- brainiak/fcma/cython_blas.pyx | 109 ++++++++++++-------- brainiak/fcma/fcma.py | 154 +---------------------------- examples/fcma/requirements.txt | 2 - examples/fcma/run.py | 153 ++++++++++++++++++++++++++++ run-tests.sh | 3 +- tests/fcma/test_voxel_selection.py | 45 +++++++++ 6 files changed, 269 insertions(+), 197 deletions(-) create mode 100644 tests/fcma/test_voxel_selection.py diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 649061642..b1ff53906 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -18,9 +18,9 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a py_start_voxel, py_ldb, py_beta, py_c, py_ldc, py_start_epoch): """ use blas API wrapped by scipy.linalg.cython_blas to compute 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 - + 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. The resulting matrix in shape [num_assigned_voxels, num_voxels] is stored in an alternate way to make sure that the correlation vectors of the same voxel stored continuously @@ -28,45 +28,58 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a Parameters ---------- py_trans_a: str - do transpose or not for the first matrix A + do transpose or not for the first matrix A + py_trans_b: str - do transpose or not for the first matrix B + do transpose or not for the first matrix B + py_m: int - the row of the resulting matrix C - in our case, is num_voxels + the row of the resulting matrix C + in our case, is num_voxels + py_n: int - the column of the resulting matrix C - in our case, is num_assigned_voxels + the column of the resulting matrix C + in our case, is num_assigned_voxels + 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 + 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 first matrix A + the weight applied to the first matrix A + py_a: 2D array in shape [epoch_length, num_voxels] in our case - the activity data of an epoch + the activity data of an epoch + py_lda: int - the stride of the first matrix A + the stride of the first matrix A + py_start_voxel: int - the starting voxel of assigned voxels - used to locate the second matrix B + the starting voxel of assigned voxels + used to locate the second matrix B + py_ldb: int - the stride of the second matrix B + the stride of the second matrix B + py_beta: float - the weight applied to the resulting matrix C + the weight applied to the resulting matrix C + py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] - place to store the resulting correlation values + place to store the resulting correlation values + py_ldc: int - the stride of the resulting matrix - in our case, num_voxels*num_epochs + the stride of the resulting matrix + in our case, num_voxels*num_epochs + py_start_epoch: int - the epoch over which the correlation is computed + the epoch over which the correlation is computed Returns ------- py_c: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] - write the resulting correlation values in an alternate way - for the processing epoch + write the resulting correlation values in an alternate way + for the processing epoch """ cdef bytes by_trans_a=py_trans_a.encode() cdef bytes by_trans_b=py_trans_b.encode() @@ -108,38 +121,48 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_star Parameters ---------- py_uplo: str - getting the upper or lower triangle of the matrix + getting the upper or lower triangle of the matrix + py_trans: str - do transpose or not for the input matrix A + do transpose or not for the input matrix A + py_n: int - the row and column of the resulting matrix C - in our case, is num_epochs + the row and column of the resulting matrix C + in our case, is num_epochs + 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_voxels + 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_voxels + py_alpha: float - the weight applied to the input matrix A + the weight applied to the input matrix A + py_a: 3D array in shape [num_assigned_voxels, num_epochs, num_voxels] in our case - the normalized correlation values of a voxel + the normalized correlation values of a voxel + py_start_voxel: int - the processed voxel - used to locate the input matrix A + the processed voxel + used to locate the input matrix A + py_lda: int - the stride of the input matrix A + the stride of the input matrix A + py_beta: float - the weight applied to the resulting matrix C + the weight applied to the resulting matrix C + py_c: 2D array in shape [num_epochs, num_epochs] - place to store the resulting kernel matrix + place to store the resulting kernel matrix + py_ldc: int - the stride of the resulting matrix + the stride of the resulting matrix Returns ------- py_c: 2D array in shape [num_epochs, num_epochs] - write the resulting kernel_matrix - for the processing voxel + write the resulting kernel_matrix + for the processing voxel """ cdef bytes by_uplo=py_uplo.encode() cdef bytes by_trans=py_trans.encode() diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/fcma.py index 5975e35aa..ef8bb541b 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/fcma.py @@ -34,9 +34,6 @@ # (Intel Labs), 2016 import numpy as np -import nibabel as nib -import os -import math import sys import time from mpi4py import MPI @@ -50,153 +47,6 @@ TERMINATETAG = 1 -def readActivityData(dir, file_extension, mask_file): - """ read data in NIfTI format and apply the spatial mask to them - - Parameters - ---------- - dir: str - the path to all subject files - file_extension: str - the file extension, usually nii.gz or nii - mask_file: str - the absolute path of the mask file, we apply the mask right after - reading a file for saving memory - - Returns - ------- - activity_data: list of 2D array in shape [nTRs, nVoxels] - the masked activity data organized in TR*voxel formats - len(activity_data) equals the number of subjects - """ - time1 = time.time() - mask_img = nib.load(mask_file) - mask = mask_img.get_data() - count = 0 - for index in np.ndindex(mask.shape): - if mask[index] != 0: - count += 1 - files = [f for f in sorted(os.listdir(dir)) - if os.path.isfile(os.path.join(dir, f)) - and f.endswith(file_extension)] - activity_data = [] - for f in files: - img = nib.load(os.path.join(dir, f)) - data = img.get_data() - (d1, d2, d3, d4) = data.shape - masked_data = np.zeros([d4, count], np.float32, order='C') - count1 = 0 - for index in np.ndindex(mask.shape): - if mask[index] != 0: - masked_data[:, count1] = np.copy(data[index]) - count1 += 1 - activity_data.append(masked_data) - print(f, masked_data.shape) - sys.stdout.flush() - time2 = time.time() - print('data reading done, takes', round(time2 - time1, 2), 's') - sys.stdout.flush() - return activity_data - - -def separateEpochs(activity_data, epoch_list): - """ separate data into epochs of interest specified in epoch_list - and z-score them for computing correlation - - Parameters - ---------- - activity_data: list of 2D array in shape [nTRs, nVoxels] - the masked activity data organized in TR*voxel formats of all subjects - epoch_list: list of 3D array in shape [condition, nEpochs, nTRs] - specification of epochs and conditions - assuming all subjects have the same number of epochs - len(epoch_list) equals the number of subjects - - Returns - ------- - raw_data: list of 2D array in shape [epoch length, nVoxels] - the data organized in epochs - and z-scored in preparation of correlation computation - len(raw_data) equals the number of epochs - labels: list of 1D array - the condition labels of the epochs - len(labels) labels equals the number of epochs - """ - time1 = time.time() - raw_data = [] - labels = [] - for sid in range(len(epoch_list)): - epoch = epoch_list[sid] - for cond in range(epoch.shape[0]): - sub_epoch = epoch[cond, :, :] - for eid in range(epoch.shape[1]): - r = np.sum(sub_epoch[eid, :]) - if r > 0: # there is an epoch in this condition - # mat is row-major - # regardless of the order of acitvity_data[sid] - mat = activity_data[sid][sub_epoch[eid, :] == 1, :] - mat = zscore(mat, axis=0, ddof=0) - # if zscore fails (standard deviation is zero), - # set all values to be zero - mat = np.nan_to_num(mat) - mat = mat / math.sqrt(r) - raw_data.append(mat) - labels.append(cond) - time2 = time.time() - print('epoch separation done, takes', round(time2 - time1, 2), 's') - sys.stdout.flush() - return raw_data, labels - - -def prepareData(data_dir, extension, mask_file, epoch_file): - """ read the data in and generate epochs of interests, - then broadcast to all workers - - Parameters - ---------- - data_dir: str - the path to all subject files - extension: str - the file extension, usually nii.gz or nii - mask_file: str - the absolute path of the mask file, - we apply the mask right after reading a file for saving memory - epoch_file: str - the absolute path of the epoch file - - Returns - ------- - raw_data: list of 2D array in shape [epoch length, nVoxels] - the data organized in epochs - len(raw_data) equals the number of epochs - labels: list of 1D array - the condition labels of the epochs - len(labels) labels equals the number of epochs - """ - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - labels = [] - raw_data = [] - if rank == 0: - activity_data = readActivityData(data_dir, extension, mask_file) - # a list of numpy array in shape [condition, nEpochs, nTRs] - epoch_list = np.load(epoch_file) - raw_data, labels = separateEpochs(activity_data, epoch_list) - time1 = time.time() - raw_data_length = len(raw_data) - raw_data_length = comm.bcast(raw_data_length, root=0) - # broadcast the data subject by subject to prevent size overflow - for i in range(raw_data_length): - if rank != 0: - raw_data.append(None) - raw_data[i] = comm.bcast(raw_data[i], root=0) - labels = comm.bcast(labels, root=0) - if rank == 0: - time2 = time.time() - print('data broadcasting done, takes', round(time2 - time1, 2), 's') - return raw_data, labels - - class VoxelSelector: """Correlation-based voxel selection component of FCMA @@ -208,7 +58,9 @@ class VoxelSelector: 2. the activity data has been z-scored, ready to compute correlation as matrix multiplication 3. all subjects have the same number of epochs - 4. voxel selection is always done in the auto-correlation, + 4. epochs belonging to the same subject are adjacent + in the list + 5. voxel selection is always done in the auto-correlation, i.e. raw_data correlate with themselves epochs_per_subj: int diff --git a/examples/fcma/requirements.txt b/examples/fcma/requirements.txt index 8808a066e..e69de29bb 100644 --- a/examples/fcma/requirements.txt +++ b/examples/fcma/requirements.txt @@ -1,2 +0,0 @@ -mpi4py -sklearn diff --git a/examples/fcma/run.py b/examples/fcma/run.py index 1911b4051..7b6cd1e89 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -13,9 +13,162 @@ # limitations under the License. from brainiak.fcma.fcma import * +from scipy.stats.mstats import zscore from sklearn import svm import sys from mpi4py import MPI +import nibabel as nib +import os +import math + + +def readActivityData(dir, file_extension, mask_file): + """ read data in NIfTI format and apply the spatial mask to them + + Parameters + ---------- + dir: str + the path to all subject files + file_extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, we apply the mask right after + reading a file for saving memory + + Returns + ------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats + len(activity_data) equals the number of subjects + """ + time1 = time.time() + mask_img = nib.load(mask_file) + mask = mask_img.get_data() + count = 0 + for index in np.ndindex(mask.shape): + if mask[index] != 0: + count += 1 + files = [f for f in sorted(os.listdir(dir)) + if os.path.isfile(os.path.join(dir, f)) + and f.endswith(file_extension)] + activity_data = [] + for f in files: + img = nib.load(os.path.join(dir, f)) + data = img.get_data() + (d1, d2, d3, d4) = data.shape + masked_data = np.zeros([d4, count], np.float32, order='C') + count1 = 0 + for index in np.ndindex(mask.shape): + if mask[index] != 0: + masked_data[:, count1] = np.copy(data[index]) + count1 += 1 + activity_data.append(masked_data) + print(f, masked_data.shape) + sys.stdout.flush() + time2 = time.time() + print('data reading done, takes', round(time2 - time1, 2), 's') + sys.stdout.flush() + return activity_data + + +def separateEpochs(activity_data, epoch_list): + """ separate data into epochs of interest specified in epoch_list + and z-score them for computing correlation + + Parameters + ---------- + activity_data: list of 2D array in shape [nTRs, nVoxels] + the masked activity data organized in TR*voxel formats of all subjects + epoch_list: list of 3D array in shape [condition, nEpochs, nTRs] + specification of epochs and conditions + assuming all subjects have the same number of epochs + len(epoch_list) equals the number of subjects + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + and z-scored in preparation of correlation computation + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs + """ + time1 = time.time() + raw_data = [] + labels = [] + for sid in range(len(epoch_list)): + epoch = epoch_list[sid] + for cond in range(epoch.shape[0]): + sub_epoch = epoch[cond, :, :] + for eid in range(epoch.shape[1]): + r = np.sum(sub_epoch[eid, :]) + if r > 0: # there is an epoch in this condition + # mat is row-major + # regardless of the order of acitvity_data[sid] + mat = activity_data[sid][sub_epoch[eid, :] == 1, :] + mat = zscore(mat, axis=0, ddof=0) + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) + mat = mat / math.sqrt(r) + raw_data.append(mat) + labels.append(cond) + time2 = time.time() + print('epoch separation done, takes', round(time2 - time1, 2), 's') + sys.stdout.flush() + return raw_data, labels + + +def prepareData(data_dir, extension, mask_file, epoch_file): + """ read the data in and generate epochs of interests, + then broadcast to all workers + + Parameters + ---------- + data_dir: str + the path to all subject files + extension: str + the file extension, usually nii.gz or nii + mask_file: str + the absolute path of the mask file, + we apply the mask right after reading a file for saving memory + epoch_file: str + the absolute path of the epoch file + + Returns + ------- + raw_data: list of 2D array in shape [epoch length, nVoxels] + the data organized in epochs + len(raw_data) equals the number of epochs + labels: list of 1D array + the condition labels of the epochs + len(labels) labels equals the number of epochs + """ + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + labels = [] + raw_data = [] + if rank == 0: + activity_data = readActivityData(data_dir, extension, mask_file) + # a list of numpy array in shape [condition, nEpochs, nTRs] + epoch_list = np.load(epoch_file) + raw_data, labels = separateEpochs(activity_data, epoch_list) + time1 = time.time() + raw_data_length = len(raw_data) + raw_data_length = comm.bcast(raw_data_length, root=0) + # broadcast the data subject by subject to prevent size overflow + for i in range(raw_data_length): + if rank != 0: + raw_data.append(None) + raw_data[i] = comm.bcast(raw_data[i], root=0) + labels = comm.bcast(labels, root=0) + if rank == 0: + time2 = time.time() + print('data broadcasting done, takes', round(time2 - time1, 2), 's') + return raw_data, labels + + """ example running command: mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz diff --git a/run-tests.sh b/run-tests.sh index a88550c2d..95194ba6c 100755 --- a/run-tests.sh +++ b/run-tests.sh @@ -20,4 +20,5 @@ pip freeze | grep -qi /brainiak || { exit 1 } -py.test --cov=brainiak +mpirun -np 2 py.test --cov=brainiak --cov-append || \ + mpirun -np 2 py.test --cov=brainiak --cov-append diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py new file mode 100644 index 000000000..18d4099e5 --- /dev/null +++ b/tests/fcma/test_voxel_selection.py @@ -0,0 +1,45 @@ +# 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. + +from brainiak.fcma.fcma import * +from scipy.stats.mstats import zscore +from sklearn import svm +import numpy as np +import math + + +def create_epoch(): + row = 12 + col = 1000 + mat = np.random.rand(row, col).astype(np.float32) + mat = zscore(mat, axis=0, ddof=0) + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) + mat = mat / math.sqrt(mat.shape[0]) + return mat + +def test_voxel_selection(): + fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch()] + labels = [0, 1, 0, 1] + vs = VoxelSelector(fake_raw_data, 2, labels, 2) + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=10 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) + results = vs.run(clf) + # test scipy normalization + fake_corr = np.random.rand(1, 12, 100).astype(np.float32) + fake_corr = vs.correlationNormalization(fake_corr) + return results, fake_corr From 193ce84d6e05d22620775bdfd6dad9d2a1a36606 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 31 Aug 2016 19:56:41 -0700 Subject: [PATCH 25/59] make two processes of fcma test exit in different time --- setup.py | 12 ++++++------ tests/fcma/test_voxel_selection.py | 5 +++++ 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/setup.py b/setup.py index c2ccf6bf9..d35262e6d 100644 --- a/setup.py +++ b/setup.py @@ -106,12 +106,6 @@ class BuildExt(build_ext): if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', '-ftemplate-depth-1024'] - # configuration for CMake - try: - out = subprocess.check_output(['cmake', '--version']) - except OSError: - raise RuntimeError("CMake must be installed to build the following extensions: " + - ", ".join(e.name for e in self.extensions)) def build_extensions(self): """the system will execute the run functions of the base class (distutils.command.build_ext) @@ -139,6 +133,12 @@ def regular_compiling(self, ext): self.build_extension(ext) def cmake_compiling(self, ext): + # configuration for CMake + try: + out = subprocess.check_output(['cmake', '--version']) + except OSError: + raise RuntimeError("CMake must be installed to build the following extensions: " + + ", ".join(e.name for e in self.extensions)) extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, '-DPYTHON_EXECUTABLE=' + sys.executable] diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index 18d4099e5..02fbadf0f 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -17,6 +17,8 @@ from sklearn import svm import numpy as np import math +from mpi4py import MPI +import time def create_epoch(): @@ -42,4 +44,7 @@ def test_voxel_selection(): # test scipy normalization fake_corr = np.random.rand(1, 12, 100).astype(np.float32) fake_corr = vs.correlationNormalization(fake_corr) + # make one process sleep a while to resolve file writing competition + if MPI.COMM_WORLD.Get_rank() == 0: + time.sleep(0.5) return results, fake_corr From 7dc574db6304558b695ca1ca6dd7205366b26a86 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 31 Aug 2016 21:18:06 -0700 Subject: [PATCH 26/59] add data type to index for better performance on Cython code; allow run-time Cython code build in run-test.sh --- brainiak/fcma/cython_blas.pyx | 16 ++++++++-------- examples/fcma/run.py | 1 + run-tests.sh | 2 ++ 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index b1ff53906..4aa2c5eb6 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -15,7 +15,7 @@ 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, - py_start_voxel, py_ldb, py_beta, py_c, py_ldc, py_start_epoch): + int py_start_voxel, py_ldb, py_beta, py_c, py_ldc, int py_start_epoch): """ use blas API wrapped by scipy.linalg.cython_blas to compute correlation The blas APIs process matrices in column-major, @@ -99,10 +99,10 @@ def compute_correlation(py_trans_a, py_trans_b, py_m, py_n, py_k, py_alpha, py_a A = py_a cdef float[:, :, ::1] C C = py_c - blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0,0], &lda, - &A[0,py_start_voxel], &ldb, &beta, &C[0,py_start_epoch,0], &ldc) + blas.sgemm(trans_a, trans_b, &M, &N, &K, &alpha, &A[0, 0], &lda, + &A[0, py_start_voxel], &ldb, &beta, &C[0, py_start_epoch,0], &ldc) -def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_start_voxel, py_lda, +def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_start_voxel, py_lda, py_beta, py_c, py_ldc): """ use blas API wrapped by scipy.linalg.cython_blas to compute kernel matrix of SVM @@ -180,16 +180,16 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, py_star A = py_a cdef float[:, ::1] C C = py_c - blas.ssyrk(uplo, trans, &N, &K, &alpha, &A[py_start_voxel,0,0], &lda, - &beta, &C[0,0], &ldc) + 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 py_c *= .001 # complete the other half of the kernel matrix if (py_uplo=='L'): for j in range(py_c.shape[0]): for k in range(j): - py_c[j,k] = py_c[k,j] + py_c[j, k] = py_c[k, j] else: for j in range(py_c.shape[0]): for k in range(j): - py_c[k,j] = py_c[j,k] + py_c[k, j] = py_c[j, k] diff --git a/examples/fcma/run.py b/examples/fcma/run.py index 7b6cd1e89..c869e2eae 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -166,6 +166,7 @@ def prepareData(data_dir, extension, mask_file, epoch_file): if rank == 0: time2 = time.time() print('data broadcasting done, takes', round(time2 - time1, 2), 's') + sys.stdout.flush() return raw_data, labels diff --git a/run-tests.sh b/run-tests.sh index 95194ba6c..677d0deb4 100755 --- a/run-tests.sh +++ b/run-tests.sh @@ -20,5 +20,7 @@ pip freeze | grep -qi /brainiak || { exit 1 } +sleep 5 # allow the run-time build of cython to complete + mpirun -np 2 py.test --cov=brainiak --cov-append || \ mpirun -np 2 py.test --cov=brainiak --cov-append From 7d91308644385289fc2ae7022322f4d141cacc25 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 31 Aug 2016 21:46:16 -0700 Subject: [PATCH 27/59] increase time for cython compiling --- run-tests.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run-tests.sh b/run-tests.sh index 677d0deb4..f990ca73c 100755 --- a/run-tests.sh +++ b/run-tests.sh @@ -20,7 +20,7 @@ pip freeze | grep -qi /brainiak || { exit 1 } -sleep 5 # allow the run-time build of cython to complete +sleep 10 # allow the run-time build of cython to complete mpirun -np 2 py.test --cov=brainiak --cov-append || \ mpirun -np 2 py.test --cov=brainiak --cov-append From ffe9c44ed58df2f430d4764577d64cb2df618de4 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 1 Sep 2016 11:17:43 -0700 Subject: [PATCH 28/59] change module name to fcma.voxelselector and address other comments from code review --- brainiak/fcma/.gitignore | 2 + brainiak/fcma/{fcma.py => voxelselector.py} | 52 +++++++++++---------- examples/fcma/run.py | 4 +- pr-check.sh | 1 + run-tests.sh | 2 - tests/fcma/test_voxel_selection.py | 4 +- 6 files changed, 35 insertions(+), 30 deletions(-) rename brainiak/fcma/{fcma.py => voxelselector.py} (91%) diff --git a/brainiak/fcma/.gitignore b/brainiak/fcma/.gitignore index 56fa9649b..169bf4eb7 100644 --- a/brainiak/fcma/.gitignore +++ b/brainiak/fcma/.gitignore @@ -1 +1,3 @@ *.npy +*.so +__pycache__ diff --git a/brainiak/fcma/fcma.py b/brainiak/fcma/voxelselector.py similarity index 91% rename from brainiak/fcma/fcma.py rename to brainiak/fcma/voxelselector.py index ef8bb541b..d73b6818c 100644 --- a/brainiak/fcma/fcma.py +++ b/brainiak/fcma/voxelselector.py @@ -43,6 +43,10 @@ from . import fcma_extension from . import cython_blas as blas +__all__ = [ + "VoxelSelector", +] + WORKTAG = 0 TERMINATETAG = 1 @@ -66,9 +70,6 @@ class VoxelSelector: epochs_per_subj: int The number of epochs of each subject - num_voxels: int - The number of voxels participating in the voxel selection - labels: list of 1D array the condition labels of the epochs len(labels) labels equals the number of epochs @@ -92,7 +93,7 @@ def __init__(self, self.num_folds = num_folds self.voxel_unit = voxel_unit if self.num_voxels == 0: - raise RuntimeError("Zero processed voxels") + raise ValueError("Zero processed voxels") def run(self, clf): """ run correlation-based voxel selection in master-worker model @@ -113,16 +114,18 @@ def run(self, clf): """ rank = MPI.COMM_WORLD.Get_rank() if rank == 0: - results = self.master() + results = self._master() # Sort the voxels results.sort(key=lambda tup: tup[1], reverse=True) else: - self.worker(clf) + self._worker(clf) results = [] return results - def master(self): - """ master node's operation, assigning tasks and collecting results + def _master(self): + """ master node's operation + + Assigning tasks to workers and collecting results from them Parameters ---------- @@ -178,9 +181,10 @@ def master(self): return results - def worker(self, clf): - """ worker node's operation, - receiving tasks from the master to process and sending the result back + def _worker(self, clf): + """ worker node's operation + + Receiving tasks from the master to process and sending the result back Parameters ---------- @@ -197,9 +201,9 @@ def worker(self, clf): task = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) if status.Get_tag(): break - comm.send(self.voxelScoring(task, clf), dest=0) + comm.send(self._voxelScoring(task, clf), dest=0) - def correlationComputation(self, task): + def _correlationComputation(self, task): """ use BLAS API to do correlation computation (matrix multiplication) Parameters @@ -212,6 +216,7 @@ def correlationComputation(self, task): corr: 3D array in shape [num_selected_voxels, num_epochs, num_voxels] the correlation values of all subjects in all epochs for the assigned values, in row-major + corr[i, e, s + j] = corr[j, e, s + i] """ s = task[0] nEpochs = len(self.raw_data) @@ -231,11 +236,12 @@ def correlationComputation(self, task): count += 1 return corr - def correlationNormalization(self, corr): + def _correlationNormalization(self, corr): """ within-subject normalization - this method uses scipy.zscore to normalize the data, - but is much slower than its C++ counterpart + This method uses scipy.zscore to normalize the data, + but is much slower than its C++ counterpart。 + It is doing in-place z-score. Parameters ---------- @@ -263,7 +269,7 @@ def correlationNormalization(self, corr): corr = np.nan_to_num(corr) return corr - def crossValidation(self, task, corr, clf): + def _crossValidation(self, task, corr, clf): """ voxelwise cross validation based on correlation vectors Parameters @@ -303,7 +309,7 @@ def crossValidation(self, task, corr, clf): results.append((i + task[0], scores.mean())) return results - def voxelScoring(self, task, clf): + def _voxelScoring(self, task, clf): """ voxel selection processing done in the worker node Take the task in, @@ -328,17 +334,13 @@ def voxelScoring(self, task, clf): """ time1 = time.time() # correlation computation - # corr is a 3D array in row major, - # in (selected_voxels, epochs, all_voxels) shape - # corr[i, e, s + j] = corr[j, e, s + i] - corr = self.correlationComputation(task) + corr = self._correlationComputation(task) # normalization - # in-place z-score, the result is still in corr - # corr = self.correlationNormalization(corr) + # corr = self._correlationNormalization(corr) fcma_extension.normalization(corr, self.epochs_per_subj) # cross validation - results = self.crossValidation(task, corr, clf) + results = self._crossValidation(task, corr, clf) time2 = time.time() print('task:', int(task[0] / self.voxel_unit), round(time2 - time1, 2), 's') diff --git a/examples/fcma/run.py b/examples/fcma/run.py index c869e2eae..9a7b75738 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -from brainiak.fcma.fcma import * +from brainiak.fcma.voxelselector import VoxelSelector from scipy.stats.mstats import zscore from sklearn import svm import sys @@ -20,6 +20,8 @@ import nibabel as nib import os import math +import time +import numpy as np def readActivityData(dir, file_extension, mask_file): diff --git a/pr-check.sh b/pr-check.sh index 168d161af..0ca5826a1 100755 --- a/pr-check.sh +++ b/pr-check.sh @@ -116,6 +116,7 @@ pip3 install $ignore_installed -U -e . || \ exit_with_error_and_venv "pip3 failed to install BrainIAK" # run tests +sleep 10 # allow the run-time build of cython to complete ./run-tests.sh || \ exit_with_error_and_venv "run-tests failed" diff --git a/run-tests.sh b/run-tests.sh index f990ca73c..95194ba6c 100755 --- a/run-tests.sh +++ b/run-tests.sh @@ -20,7 +20,5 @@ pip freeze | grep -qi /brainiak || { exit 1 } -sleep 10 # allow the run-time build of cython to complete - mpirun -np 2 py.test --cov=brainiak --cov-append || \ mpirun -np 2 py.test --cov=brainiak --cov-append diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index 02fbadf0f..92f2aebe2 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -from brainiak.fcma.fcma import * +from brainiak.fcma.voxelselector import VoxelSelector from scipy.stats.mstats import zscore from sklearn import svm import numpy as np @@ -43,7 +43,7 @@ def test_voxel_selection(): results = vs.run(clf) # test scipy normalization fake_corr = np.random.rand(1, 12, 100).astype(np.float32) - fake_corr = vs.correlationNormalization(fake_corr) + fake_corr = vs._correlationNormalization(fake_corr) # make one process sleep a while to resolve file writing competition if MPI.COMM_WORLD.Get_rank() == 0: time.sleep(0.5) From e02f751b406048d0cfd8832eb5bc5df785211c54 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 1 Sep 2016 12:04:48 -0700 Subject: [PATCH 29/59] use top-level .gitignore and remove an empty file --- .gitignore | 2 ++ brainiak/fcma/.gitignore | 3 --- examples/fcma/requirements.txt | 0 3 files changed, 2 insertions(+), 3 deletions(-) delete mode 100644 brainiak/fcma/.gitignore delete mode 100644 examples/fcma/requirements.txt diff --git a/.gitignore b/.gitignore index 892caef59..dbcb04915 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,7 @@ *.swp *.lib *.o +*.npy brainiak.egg-info build @@ -18,3 +19,4 @@ htmlcov .idea .ipynb_checkpoints junit-*.xml +brainiak/fcma/__pycache__ diff --git a/brainiak/fcma/.gitignore b/brainiak/fcma/.gitignore deleted file mode 100644 index 169bf4eb7..000000000 --- a/brainiak/fcma/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.npy -*.so -__pycache__ diff --git a/examples/fcma/requirements.txt b/examples/fcma/requirements.txt deleted file mode 100644 index e69de29bb..000000000 From 3ac536833cec6a2104462d57eea6247bc051eec3 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 1 Sep 2016 12:16:36 -0700 Subject: [PATCH 30/59] make mpi tags within class VoxelSelector --- brainiak/fcma/voxelselector.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index d73b6818c..d684669b3 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -47,9 +47,6 @@ "VoxelSelector", ] -WORKTAG = 0 -TERMINATETAG = 1 - class VoxelSelector: """Correlation-based voxel selection component of FCMA @@ -95,6 +92,10 @@ def __init__(self, if self.num_voxels == 0: raise ValueError("Zero processed voxels") + # tags for MPI messages + _WORKTAG = 0 + _TERMINATETAG = 1 + def run(self, clf): """ run correlation-based voxel selection in master-worker model @@ -151,7 +152,7 @@ def _master(self): if current_task[1] == 0: using_size = i break - comm.send(current_task, dest=i, tag=WORKTAG) + comm.send(current_task, dest=i, tag=self._WORKTAG) next_start = current_task[0] + current_task[1] sending_voxels = self.voxel_unit \ if self.voxel_unit < self.num_voxels - next_start \ @@ -165,7 +166,7 @@ def _master(self): tag=MPI.ANY_TAG, status=status) results += result - comm.send(current_task, dest=status.Get_source(), tag=WORKTAG) + comm.send(current_task, dest=status.Get_source(), tag=self._WORKTAG) next_start = current_task[0] + current_task[1] sending_voxels = self.voxel_unit \ if self.voxel_unit < self.num_voxels - next_start \ @@ -177,7 +178,7 @@ def _master(self): results += result for i in range(1, size): - comm.send(None, dest=i, tag=TERMINATETAG) + comm.send(None, dest=i, tag=self._TERMINATETAG) return results From b9d9b178ea9c9b096f2c757dff9e174013651003 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 1 Sep 2016 12:22:21 -0700 Subject: [PATCH 31/59] fix code format issue --- brainiak/fcma/voxelselector.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index d684669b3..cc1534d25 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -152,7 +152,9 @@ def _master(self): if current_task[1] == 0: using_size = i break - comm.send(current_task, dest=i, tag=self._WORKTAG) + comm.send(current_task, + dest=i, + tag=self._WORKTAG) next_start = current_task[0] + current_task[1] sending_voxels = self.voxel_unit \ if self.voxel_unit < self.num_voxels - next_start \ @@ -166,7 +168,9 @@ def _master(self): tag=MPI.ANY_TAG, status=status) results += result - comm.send(current_task, dest=status.Get_source(), tag=self._WORKTAG) + comm.send(current_task, + dest=status.Get_source(), + tag=self._WORKTAG) next_start = current_task[0] + current_task[1] sending_voxels = self.voxel_unit \ if self.voxel_unit < self.num_voxels - next_start \ @@ -178,7 +182,9 @@ def _master(self): results += result for i in range(1, size): - comm.send(None, dest=i, tag=self._TERMINATETAG) + comm.send(None, + dest=i, + tag=self._TERMINATETAG) return results @@ -199,10 +205,13 @@ def _worker(self, clf): comm = MPI.COMM_WORLD status = MPI.Status() while 1: - task = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) + task = comm.recv(source=0, + tag=MPI.ANY_TAG, + status=status) if status.Get_tag(): break - comm.send(self._voxelScoring(task, clf), dest=0) + comm.send(self._voxelScoring(task, clf), + dest=0) def _correlationComputation(self, task): """ use BLAS API to do correlation computation (matrix multiplication) From 45936ee9708eb9a090a603e7869a3b99e711b281 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 7 Sep 2016 15:19:45 -0700 Subject: [PATCH 32/59] use logging instead of print --- brainiak/fcma/voxelselector.py | 40 +++++++++++++++++++++++++++++++--- examples/fcma/run.py | 36 +++++++++++++++++++++++------- 2 files changed, 65 insertions(+), 11 deletions(-) diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index cc1534d25..b52fb45e4 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -42,6 +42,9 @@ import sklearn from . import fcma_extension from . import cython_blas as blas +import logging + +logger = logging.getLogger(__name__) __all__ = [ "VoxelSelector", @@ -138,6 +141,9 @@ def _master(self): the accuracy numbers of all voxels, in accuracy descending order the length of array equals the number of voxels """ + logger.info( + 'Master starts to allocate tasks' + ) results = [] comm = MPI.COMM_WORLD size = comm.Get_size() @@ -228,6 +234,7 @@ def _correlationComputation(self, task): for the assigned values, in row-major corr[i, e, s + j] = corr[j, e, s + i] """ + time1 = time.time() s = task[0] nEpochs = len(self.raw_data) corr = np.zeros((task[1], nEpochs, self.num_voxels), @@ -244,6 +251,11 @@ def _correlationComputation(self, task): 0.0, corr, self.num_voxels * nEpochs, count) count += 1 + time2 = time.time() + logger.debug( + 'correlation computation for %d voxels, takes %.2f s' % + ((task[1] - task[0]), (time2 - time1)) + ) return corr def _correlationNormalization(self, corr): @@ -265,6 +277,7 @@ def _correlationNormalization(self, corr): the normalized correlation values of all subjects in all epochs for the assigned values, in row-major """ + time1 = time.time() (sv, e, av) = corr.shape for i in range(sv): start = 0 @@ -277,6 +290,12 @@ def _correlationNormalization(self, corr): # if zscore fails (standard deviation is zero), # set all values to be zero corr = np.nan_to_num(corr) + time2 = time.time() + logger.debug( + 'within-subject normalization for %d voxels using numpy zscore function,' + ' takes %.2f s' % + (sv, (time2 - time1)) + ) return corr def _crossValidation(self, task, corr, clf): @@ -296,6 +315,7 @@ def _crossValidation(self, task, corr, clf): the accuracy numbers of all voxels, in accuracy descending order the length of array equals the number of assigned voxels """ + time1 = time.time() (sv, e, av) = corr.shape kernel_matrix = np.zeros((e, e), np.float32, order='C') results = [] @@ -317,6 +337,11 @@ def _crossValidation(self, task, corr, clf): scores = cross_validation.cross_val_score(clf, data, self.labels, cv=skf, n_jobs=1) results.append((i + task[0], scores.mean())) + time2 = time.time() + logger.debug( + 'cross validation for %d voxels, takes %.2f s' % + ((task[1] - task[0]), (time2 - time1)) + ) return results def _voxelScoring(self, task, clf): @@ -347,12 +372,21 @@ def _voxelScoring(self, task, clf): corr = self._correlationComputation(task) # normalization # corr = self._correlationNormalization(corr) + time3 = time.time() fcma_extension.normalization(corr, self.epochs_per_subj) + time4 = time.time() + logger.debug( + 'within-subject normalization for %d voxels using C++, ' + 'takes %.2f s' % + ((task[1] - task[0]), (time4 - time3)) + ) # cross validation results = self._crossValidation(task, corr, clf) time2 = time.time() - print('task:', int(task[0] / self.voxel_unit), - round(time2 - time1, 2), 's') - sys.stdout.flush() + logger.info( + 'in rank %d, task %d takes %.2f s' % + (MPI.COMM_WORLD.Get_rank(), + (int(task[0] / self.voxel_unit)), (time2 - time1)) + ) return results diff --git a/examples/fcma/run.py b/examples/fcma/run.py index 9a7b75738..3ced3c8a0 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/run.py @@ -22,7 +22,13 @@ import math import time import numpy as np +import logging +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__) def readActivityData(dir, file_extension, mask_file): """ read data in NIfTI format and apply the spatial mask to them @@ -65,11 +71,15 @@ def readActivityData(dir, file_extension, mask_file): masked_data[:, count1] = np.copy(data[index]) count1 += 1 activity_data.append(masked_data) - print(f, masked_data.shape) - sys.stdout.flush() + logger.info( + 'file %s is loaded and masked, with data shape %s' % + (f, masked_data.shape) + ) time2 = time.time() - print('data reading done, takes', round(time2 - time1, 2), 's') - sys.stdout.flush() + logger.info( + 'data reading done, takes %.2f s' % + (time2 - time1) + ) return activity_data @@ -117,8 +127,10 @@ def separateEpochs(activity_data, epoch_list): raw_data.append(mat) labels.append(cond) time2 = time.time() - print('epoch separation done, takes', round(time2 - time1, 2), 's') - sys.stdout.flush() + logger.info( + 'epoch separation done, takes %.2f s' % + (time2 - time1) + ) return raw_data, labels @@ -167,8 +179,10 @@ def prepareData(data_dir, extension, mask_file, epoch_file): labels = comm.bcast(labels, root=0) if rank == 0: time2 = time.time() - print('data broadcasting done, takes', round(time2 - time1, 2), 's') - sys.stdout.flush() + logger.info( + 'data broadcasting done, takes %.2f s' % + (time2 - time1) + ) return raw_data, labels @@ -178,6 +192,11 @@ def prepareData(data_dir, extension, mask_file, epoch_file): data/fs_epoch_labels.npy 12 18 """ if __name__ == '__main__': + if MPI.COMM_WORLD.Get_rank()==0: + logger.info( + 'programming starts in %d process(es)' % + MPI.COMM_WORLD.Get_size() + ) data_dir = sys.argv[1] extension = sys.argv[2] mask_file = sys.argv[3] @@ -190,5 +209,6 @@ def prepareData(data_dir, extension, mask_file, epoch_file): # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) results = vs.run(clf) + # this output is just for result checking if MPI.COMM_WORLD.Get_rank()==0: print(results[0:100]) From a009225e5805ba4c3fbfca87b7234073fba63744 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 7 Sep 2016 15:23:10 -0700 Subject: [PATCH 33/59] formatting --- .gitignore | 1 + brainiak/fcma/voxelselector.py | 9 ++++----- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index dbcb04915..d80d83f96 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ *.lib *.o *.npy +*.log brainiak.egg-info build diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index b52fb45e4..b9cf2d1aa 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -34,7 +34,6 @@ # (Intel Labs), 2016 import numpy as np -import sys import time from mpi4py import MPI from scipy.stats.mstats import zscore @@ -292,8 +291,8 @@ def _correlationNormalization(self, corr): corr = np.nan_to_num(corr) time2 = time.time() logger.debug( - 'within-subject normalization for %d voxels using numpy zscore function,' - ' takes %.2f s' % + 'within-subject normalization for %d voxels ' + 'using numpy zscore function, takes %.2f s' % (sv, (time2 - time1)) ) return corr @@ -376,8 +375,8 @@ def _voxelScoring(self, task, clf): fcma_extension.normalization(corr, self.epochs_per_subj) time4 = time.time() logger.debug( - 'within-subject normalization for %d voxels using C++, ' - 'takes %.2f s' % + 'within-subject normalization for %d voxels ' + 'using C++, takes %.2f s' % ((task[1] - task[0]), (time4 - time3)) ) From 5bd0fa316c3228ad9db076af6619afa00f09b0ef Mon Sep 17 00:00:00 2001 From: yidawang Date: Thu, 8 Sep 2016 15:34:05 -0400 Subject: [PATCH 34/59] add extra link flag for cython compiling under MKL --- brainiak/fcma/cython_blas.pyxbld | 7 +++++++ setup.py | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) create mode 100644 brainiak/fcma/cython_blas.pyxbld diff --git a/brainiak/fcma/cython_blas.pyxbld b/brainiak/fcma/cython_blas.pyxbld new file mode 100644 index 000000000..b3127149d --- /dev/null +++ b/brainiak/fcma/cython_blas.pyxbld @@ -0,0 +1,7 @@ +def make_ext(modname, pyxfilename): + from distutils.extension import Extension + return Extension( + name=modname, + sources=[pyxfilename], + extra_link_args=['-lirc', '-lintlc'] + ) diff --git a/setup.py b/setup.py index d35262e6d..74726b446 100644 --- a/setup.py +++ b/setup.py @@ -179,7 +179,7 @@ def cmake_compiling(self, ext): cmdclass={'build_ext': BuildExt}, packages=find_packages(), package_data = { - 'brainiak.fcma': ['*.pyx'], + 'brainiak.fcma': ['*.pyx', '*.pyxbld'], }, zip_safe=False, ) From 7f822fa4f95f6a294ab0da5ffb0f496cc90516a4 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 8 Sep 2016 12:49:13 -0700 Subject: [PATCH 35/59] add compiler specific links to cython compiling --- brainiak/fcma/cython_blas.pyxbld | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/brainiak/fcma/cython_blas.pyxbld b/brainiak/fcma/cython_blas.pyxbld index b3127149d..1c50e3c85 100644 --- a/brainiak/fcma/cython_blas.pyxbld +++ b/brainiak/fcma/cython_blas.pyxbld @@ -1,7 +1,11 @@ def make_ext(modname, pyxfilename): from distutils.extension import Extension + import os + opts = [] + if 'icc' in os.environ["CC"]: + opts += ['-lirc', '-lintlc'] return Extension( - name=modname, - sources=[pyxfilename], - extra_link_args=['-lirc', '-lintlc'] + name = modname, + sources = [pyxfilename], + extra_link_args = opts, ) From 76992600d934828060a83ba686834435ded767fc Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 8 Sep 2016 22:27:27 -0700 Subject: [PATCH 36/59] make master rank selectable for better coverage in testing --- brainiak/fcma/voxelselector.py | 34 ++++++++++--- ...selection.py => test_voxel_selection_0.py} | 3 +- tests/fcma/test_voxel_selection_1.py | 51 +++++++++++++++++++ 3 files changed, 79 insertions(+), 9 deletions(-) rename tests/fcma/{test_voxel_selection.py => test_voxel_selection_0.py} (94%) create mode 100644 tests/fcma/test_voxel_selection_1.py diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index b9cf2d1aa..51058ad3f 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -78,21 +78,33 @@ class VoxelSelector: voxel_unit: int, default 100 The number of voxel assigned to a worker each time + + master_rank: int, default 0 + The process which serves as the master """ def __init__(self, raw_data, epochs_per_subj, labels, num_folds, - voxel_unit=100): + voxel_unit=100, + master_rank=0): self.raw_data = raw_data self.epochs_per_subj = epochs_per_subj self.num_voxels = raw_data[0].shape[1] self.labels = labels self.num_folds = num_folds self.voxel_unit = voxel_unit + self.master_rank = master_rank if self.num_voxels == 0: - raise ValueError("Zero processed voxels") + raise ValueError('Zero processed voxels') + if MPI.COMM_WORLD.Get_size() == 1: + raise RuntimeError('one process cannot run the ' + 'master-worker model') + if self.master_rank >= MPI.COMM_WORLD.Get_size(): + logger.warn('Master rank exceeds the number of ' + 'launched processes, set to 0') + self.master_rank = 0 # tags for MPI messages _WORKTAG = 0 @@ -116,7 +128,7 @@ def run(self, clf): the length of array equals the number of voxels """ rank = MPI.COMM_WORLD.Get_rank() - if rank == 0: + if rank == self.master_rank: results = self._master() # Sort the voxels results.sort(key=lambda tup: tup[1], reverse=True) @@ -153,7 +165,9 @@ def _master(self): # using_size is used when the number of tasks # is smaller than the number of workers using_size = size - for i in range(1, size): + for i in range(0, size): + if i == self.master_rank: + continue if current_task[1] == 0: using_size = i break @@ -182,11 +196,15 @@ def _master(self): else self.num_voxels - next_start current_task = (next_start, sending_voxels) - for i in range(1, using_size): + for i in range(0, using_size): + if i == self.master_rank: + continue result = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG) results += result - for i in range(1, size): + for i in range(0, size): + if i == self.master_rank: + continue comm.send(None, dest=i, tag=self._TERMINATETAG) @@ -210,13 +228,13 @@ def _worker(self, clf): comm = MPI.COMM_WORLD status = MPI.Status() while 1: - task = comm.recv(source=0, + task = comm.recv(source=self.master_rank, tag=MPI.ANY_TAG, status=status) if status.Get_tag(): break comm.send(self._voxelScoring(task, clf), - dest=0) + dest=self.master_rank) def _correlationComputation(self, task): """ use BLAS API to do correlation computation (matrix multiplication) diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection_0.py similarity index 94% rename from tests/fcma/test_voxel_selection.py rename to tests/fcma/test_voxel_selection_0.py index 92f2aebe2..373b82344 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection_0.py @@ -36,7 +36,8 @@ def test_voxel_selection(): fake_raw_data = [create_epoch(), create_epoch(), create_epoch(), create_epoch()] labels = [0, 1, 0, 1] - vs = VoxelSelector(fake_raw_data, 2, labels, 2) + # set master rank to be 0 + vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=0) # for cross validation, use SVM with precomputed kernel # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) diff --git a/tests/fcma/test_voxel_selection_1.py b/tests/fcma/test_voxel_selection_1.py new file mode 100644 index 000000000..e6f61d3b5 --- /dev/null +++ b/tests/fcma/test_voxel_selection_1.py @@ -0,0 +1,51 @@ +# 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. + +from brainiak.fcma.voxelselector import VoxelSelector +from scipy.stats.mstats import zscore +from sklearn import svm +import numpy as np +import math +from mpi4py import MPI +import time + + +def create_epoch(): + row = 12 + col = 1000 + mat = np.random.rand(row, col).astype(np.float32) + mat = zscore(mat, axis=0, ddof=0) + # if zscore fails (standard deviation is zero), + # set all values to be zero + mat = np.nan_to_num(mat) + mat = mat / math.sqrt(mat.shape[0]) + return mat + +def test_voxel_selection(): + fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch()] + labels = [0, 1, 0, 1] + # set master rank to be 0 + vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=1) + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=1 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + results = vs.run(clf) + # test scipy normalization + fake_corr = np.random.rand(1, 12, 100).astype(np.float32) + fake_corr = vs._correlationNormalization(fake_corr) + # make one process sleep a while to resolve file writing competition + if MPI.COMM_WORLD.Get_rank() == 0: + time.sleep(0.5) + return results, fake_corr From b1cdbd25aafd51dc98bb88cc44c049c6e192a7fe Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Thu, 8 Sep 2016 22:39:25 -0700 Subject: [PATCH 37/59] make the test example conciser, delete the example epoch data --- examples/fcma/data/fs_epoch_labels.npy | Bin 105504 -> 0 bytes ...selection_0.py => test_voxel_selection.py} | 7 ++- tests/fcma/test_voxel_selection_1.py | 51 ------------------ 3 files changed, 5 insertions(+), 53 deletions(-) delete mode 100644 examples/fcma/data/fs_epoch_labels.npy rename tests/fcma/{test_voxel_selection_0.py => test_voxel_selection.py} (89%) delete mode 100644 tests/fcma/test_voxel_selection_1.py diff --git a/examples/fcma/data/fs_epoch_labels.npy b/examples/fcma/data/fs_epoch_labels.npy deleted file mode 100644 index 0602190b2beabf7a5cee2df4dcaedce89d8ee062..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 105504 zcmeIwv5KB$6ot{G^(v+tg)|8%#0#)fY!s{{n23THN$k|C@Plw-IG5bU%z?(KNjGDgDk*9)-dzHKo;O3YnWRBPFk7GJrA4j-md1MmD$|$u=(!oY93md&3*IG z_M4Sl%>x5jfQPJM=7E7Mz(dwB^T5ER04JA9JKylI^VXM3dAL;C`G$v`x4vA;!==*B zH_b!a@04+A7T_Ukn0a6z3-FLN%seoV1$f9B=2n1{R%Ub0!{)oUt9fW;HupSizI(fx zhgN2D-#oPaW+hkiz(5w@A#0d*Ar950J?R?WbwEa#Qmu3MTvWA%l2C@JTS;Nc&16hEFtYK~iIB8`z_dIOA zd%K#4R%Ub0!{)oUt9fW;Huud#+izBKH4hAA0Uol3nFj{401sKi%mV|L0-Rhb?R>+- z&RbtD<>69k=Nle&-uiMW50^?i-!uEWkt7F!R7b7T_Ukm|Fo( zTA9r~51a4auI8bY+1&H6`R?s%9$J~ree=-vo0VM60|Qxrhpb`dfq^W*L)I|!z`&&d zCznb)-|(>W)|X3pxK!HthKHTEzFf+~rP9tf%|qMolyPYm;2~?6d0-$5@Q^jkJTQ<2 zc*q*&R)CXMW^>QO=DWA6d1z%e_dIOAd%K#4R%UbGJhc60C0FynKo;O3YnXXpAPex2 zHOxFPa4EparP9tfJnX#nQO=DWA6d1z%e_sv7wZ&q?O z4-8}h99oA2JP=Ao6@-1D&c?(J$GTA9s# z^U(I2m0ZmO16hEFtYPMXfh@p7)-dzHz@-2umr6U|@UZjNmrHrLRNDE5hn=^+T*||x K(#|*k=iy(7!Xh^S diff --git a/tests/fcma/test_voxel_selection_0.py b/tests/fcma/test_voxel_selection.py similarity index 89% rename from tests/fcma/test_voxel_selection_0.py rename to tests/fcma/test_voxel_selection.py index 373b82344..ccc72c04e 100644 --- a/tests/fcma/test_voxel_selection_0.py +++ b/tests/fcma/test_voxel_selection.py @@ -41,11 +41,14 @@ def test_voxel_selection(): # for cross validation, use SVM with precomputed kernel # no shrinking, set C=10 clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) - results = vs.run(clf) + results0 = vs.run(clf) + # set master rank to be 1 and do it again + vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=1) + results1 = vs.run(clf) # test scipy normalization fake_corr = np.random.rand(1, 12, 100).astype(np.float32) fake_corr = vs._correlationNormalization(fake_corr) # make one process sleep a while to resolve file writing competition if MPI.COMM_WORLD.Get_rank() == 0: time.sleep(0.5) - return results, fake_corr + return results0, results1, fake_corr diff --git a/tests/fcma/test_voxel_selection_1.py b/tests/fcma/test_voxel_selection_1.py deleted file mode 100644 index e6f61d3b5..000000000 --- a/tests/fcma/test_voxel_selection_1.py +++ /dev/null @@ -1,51 +0,0 @@ -# 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. - -from brainiak.fcma.voxelselector import VoxelSelector -from scipy.stats.mstats import zscore -from sklearn import svm -import numpy as np -import math -from mpi4py import MPI -import time - - -def create_epoch(): - row = 12 - col = 1000 - mat = np.random.rand(row, col).astype(np.float32) - mat = zscore(mat, axis=0, ddof=0) - # if zscore fails (standard deviation is zero), - # set all values to be zero - mat = np.nan_to_num(mat) - mat = mat / math.sqrt(mat.shape[0]) - return mat - -def test_voxel_selection(): - fake_raw_data = [create_epoch(), create_epoch(), - create_epoch(), create_epoch()] - labels = [0, 1, 0, 1] - # set master rank to be 0 - vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=1) - # for cross validation, use SVM with precomputed kernel - # no shrinking, set C=1 - clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) - results = vs.run(clf) - # test scipy normalization - fake_corr = np.random.rand(1, 12, 100).astype(np.float32) - fake_corr = vs._correlationNormalization(fake_corr) - # make one process sleep a while to resolve file writing competition - if MPI.COMM_WORLD.Get_rank() == 0: - time.sleep(0.5) - return results, fake_corr From 926cb66ca43bfb335b5c73caa8f2794bb36ad256 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 9 Sep 2016 21:41:01 -0700 Subject: [PATCH 38/59] better c compiler detection in the pyxbld file --- brainiak/fcma/cython_blas.pyxbld | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/brainiak/fcma/cython_blas.pyxbld b/brainiak/fcma/cython_blas.pyxbld index 1c50e3c85..3838602f9 100644 --- a/brainiak/fcma/cython_blas.pyxbld +++ b/brainiak/fcma/cython_blas.pyxbld @@ -1,8 +1,12 @@ def make_ext(modname, pyxfilename): from distutils.extension import Extension + from distutils import sysconfig import os opts = [] - if 'icc' in os.environ["CC"]: + if 'CC' in os.environ: + if 'icc' in os.environ['CC']: + opts += ['-lirc', '-lintlc'] + elif 'icc' in sysconfig.get_config_var('CC'): opts += ['-lirc', '-lintlc'] return Extension( name = modname, From 1365da858b68cea7dc0d938fe63d7a6d432c9db2 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 9 Sep 2016 23:36:30 -0700 Subject: [PATCH 39/59] compile cython_blas in __init__.py using only one process, no need to sleep in pr-check before py.test --- brainiak/fcma/__init__.py | 11 +++++++++++ brainiak/fcma/cython_blas.pyx | 5 +++++ brainiak/fcma/voxelselector.py | 2 ++ pr-check.sh | 1 - 4 files changed, 18 insertions(+), 1 deletion(-) diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py index 3880ca640..79b91e2ca 100644 --- a/brainiak/fcma/__init__.py +++ b/brainiak/fcma/__init__.py @@ -14,4 +14,15 @@ """Full correlation matrix analysis""" import pyximport +from mpi4py import MPI pyximport.install() + +a = 'uninstalled' +# make the cython building execute only once +if MPI.COMM_WORLD.Get_rank() == 0: + from . import cython_blas as blas + blas.installed() + a = 'installed' +# tell everybody the building is done +# functions as a barrier +a = MPI.COMM_WORLD.bcast(a, root=0) diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 4aa2c5eb6..84ccb6925 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -193,3 +193,8 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ for j in range(py_c.shape[0]): for k in range(j): py_c[k, j] = py_c[j, k] + +def installed(): + """ + This is an empty method for installing cython_blas library + """ diff --git a/brainiak/fcma/voxelselector.py b/brainiak/fcma/voxelselector.py index 51058ad3f..fc0a747e6 100644 --- a/brainiak/fcma/voxelselector.py +++ b/brainiak/fcma/voxelselector.py @@ -65,6 +65,8 @@ class VoxelSelector: in the list 5. voxel selection is always done in the auto-correlation, i.e. raw_data correlate with themselves + 6. if MPI jobs are running on multiple nodes, the user + home directory is shared by all nodes epochs_per_subj: int The number of epochs of each subject diff --git a/pr-check.sh b/pr-check.sh index 0ca5826a1..168d161af 100755 --- a/pr-check.sh +++ b/pr-check.sh @@ -116,7 +116,6 @@ pip3 install $ignore_installed -U -e . || \ exit_with_error_and_venv "pip3 failed to install BrainIAK" # run tests -sleep 10 # allow the run-time build of cython to complete ./run-tests.sh || \ exit_with_error_and_venv "run-tests failed" From 1260a42207cef771e99fd8c009390416a02f1821 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 09:27:19 -0700 Subject: [PATCH 40/59] add run-tests.sh back --- run-tests.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run-tests.sh b/run-tests.sh index d8ae83d33..f20c25f41 100644 --- a/run-tests.sh +++ b/run-tests.sh @@ -33,4 +33,4 @@ $mpi_command -n 2 coverage run -m pytest coverage combine coverage report coverage html -coverage xml \ No newline at end of file +coverage xml From b5f2c0e75ff2c1f340aa195a4b986c42deb1b07d Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 09:30:28 -0700 Subject: [PATCH 41/59] chmod of run-tests.sh --- run-tests.sh | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 run-tests.sh diff --git a/run-tests.sh b/run-tests.sh old mode 100644 new mode 100755 From 9e47e15b6e23121b4286405b797a3db170518e1f Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 14:23:22 -0700 Subject: [PATCH 42/59] remove cmake related info for this PR --- brainiak/fcma/CMakeLists.txt | 38 ------------------------------------ setup.py | 37 ++--------------------------------- 2 files changed, 2 insertions(+), 73 deletions(-) delete mode 100644 brainiak/fcma/CMakeLists.txt diff --git a/brainiak/fcma/CMakeLists.txt b/brainiak/fcma/CMakeLists.txt deleted file mode 100644 index e675bdf1d..000000000 --- a/brainiak/fcma/CMakeLists.txt +++ /dev/null @@ -1,38 +0,0 @@ -cmake_minimum_required (VERSION 2.6) -project (FCMA_EXT) - -if ($ENV{VERBOSE}) - set (CMAKE_VERBOSE_MAKEFILE ON) -endif ($ENV{VERBOSE}) - -# add the binary tree to the search path for include files -include_directories("${PROJECT_BINARY_DIR}") - -#set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") -#set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") -# set (CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") -set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fopenmp -std=c++14 -fPIC -flto") - -#if ($ENV{USE_MKL}) -# set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DUSE_MKL") -# set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -mkl=parallel") -#else ($ENV{USE_MKL}) -# set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lopenblas -lomp") -# set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lopenblas -lomp") -#endif ($ENV{USE_MKL}) - -set(CMAKE_C_COMPILER ${CC}) -set(CMAKE_CXX_COMPILER ${CXX}) - -if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "clang") - set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lomp") -elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") - set (CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -lgomp") -endif() - -# get all *.cc files recursively -file(GLOB_RECURSE SRC_FILES src/*.cc) - -# add the library -add_subdirectory(pybind11) -pybind11_add_module(fcma_extension ${SRC_FILES}) \ No newline at end of file diff --git a/setup.py b/setup.py index 74726b446..efd2e1798 100644 --- a/setup.py +++ b/setup.py @@ -33,11 +33,6 @@ def __str__(self): import pybind11 return pybind11.get_include(self.user) -class CMakeExtension(Extension): - def __init__(self, name, sourcedir=''): - Extension.__init__(self, name, sources=[]) - self.sourcedir = os.path.abspath(sourcedir) - ext_modules = [ Extension( 'brainiak.factoranalysis.tfa_extension', @@ -48,7 +43,6 @@ def __init__(self, name, sourcedir=''): get_pybind_include(user=True) ], ), - #CMakeExtension('brainiak.fcma.fcma_extension', 'brainiak/fcma/'), Extension( 'brainiak.fcma.fcma_extension', ['brainiak/fcma/src/fcma_extension.cc'], @@ -93,7 +87,7 @@ def cpp_flag(compiler): class BuildExt(build_ext): """A custom build extension for adding compiler-specific options. - Support both regular compiling and CMake compiling. + Support regular compiling. Support Unix/Linux and MacOS. """ if sys.platform == 'Windows': @@ -114,9 +108,7 @@ def build_extensions(self): # First, sanity-check the 'extensions' list, i.e. ext_module self.check_extensions_list(self.extensions) for ext in self.extensions: - if isinstance(ext, CMakeExtension): - self.cmake_compiling(ext) - elif isinstance(ext, Extension): + if isinstance(ext, Extension): self.regular_compiling(ext) def regular_compiling(self, ext): @@ -132,31 +124,6 @@ def regular_compiling(self, ext): ext.extra_link_args = opts self.build_extension(ext) - def cmake_compiling(self, ext): - # configuration for CMake - try: - out = subprocess.check_output(['cmake', '--version']) - except OSError: - raise RuntimeError("CMake must be installed to build the following extensions: " + - ", ".join(e.name for e in self.extensions)) - extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) - cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, - '-DPYTHON_EXECUTABLE=' + sys.executable] - - cfg = 'Debug' if self.debug else 'Release' - build_args = ['--config', cfg] - - cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] - build_args += ['--', '-j'] - - env = os.environ.copy() - env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format(env.get('CXXFLAGS', ''), - self.distribution.get_version()) - if not os.path.exists(self.build_temp): - os.makedirs(self.build_temp) - subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, cwd=self.build_temp, env=env) - subprocess.check_call(['cmake', '--build', '.'] + build_args, cwd=self.build_temp) - setup( name='brainiak', version=__version__, From a77b51fede089208b629ea67498490c0dc5c4b4f Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 19:14:29 -0700 Subject: [PATCH 43/59] restore original setup.py and address other review comments --- .gitignore | 2 +- brainiak/fcma/cython_blas.pyxbld | 6 ++--- setup.py | 43 ++++++-------------------------- 3 files changed, 10 insertions(+), 41 deletions(-) diff --git a/.gitignore b/.gitignore index d80d83f96..39739d1ae 100644 --- a/.gitignore +++ b/.gitignore @@ -20,4 +20,4 @@ htmlcov .idea .ipynb_checkpoints junit-*.xml -brainiak/fcma/__pycache__ +__pycache__ diff --git a/brainiak/fcma/cython_blas.pyxbld b/brainiak/fcma/cython_blas.pyxbld index 3838602f9..45ab31910 100644 --- a/brainiak/fcma/cython_blas.pyxbld +++ b/brainiak/fcma/cython_blas.pyxbld @@ -3,11 +3,9 @@ def make_ext(modname, pyxfilename): from distutils import sysconfig import os opts = [] - if 'CC' in os.environ: - if 'icc' in os.environ['CC']: + if ('CC' in os.environ and 'icc' in os.environ['CC']) or \ + 'icc' in sysconfig.get_config_var('CC'): opts += ['-lirc', '-lintlc'] - elif 'icc' in sysconfig.get_config_var('CC'): - opts += ['-lirc', '-lintlc'] return Extension( name = modname, sources = [pyxfilename], diff --git a/setup.py b/setup.py index efd2e1798..8ad6aeb99 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,6 @@ import os import sys import setuptools -import subprocess __version__ = '0.2' @@ -21,7 +20,6 @@ class get_pybind_include(object): """Helper class to determine the pybind11 include path - The purpose of this class is to postpone importing pybind11 until it is actually installed, so that the ``get_include()`` method can be invoked. """ @@ -33,6 +31,7 @@ def __str__(self): import pybind11 return pybind11.get_include(self.user) + ext_modules = [ Extension( 'brainiak.factoranalysis.tfa_extension', @@ -43,15 +42,6 @@ def __str__(self): get_pybind_include(user=True) ], ), - Extension( - 'brainiak.fcma.fcma_extension', - ['brainiak/fcma/src/fcma_extension.cc'], - include_dirs=[ - # Path to pybind11 headers - get_pybind_include(), - get_pybind_include(user=True) - ], - ), ] @@ -73,7 +63,6 @@ def has_flag(compiler, flagname): def cpp_flag(compiler): """Return the -std=c++[11/14] compiler flag. - The c++14 is prefered over c++11 (when it is available). """ if has_flag(compiler, '-std=c++14'): @@ -86,32 +75,16 @@ def cpp_flag(compiler): class BuildExt(build_ext): - """A custom build extension for adding compiler-specific options. - Support regular compiling. - Support Unix/Linux and MacOS. - """ - if sys.platform == 'Windows': - raise RuntimeError("BrainIAK cannot be built on Windows") - - # configuration for regular compiling + """A custom build extension for adding compiler-specific options.""" c_opts = { 'unix': ['-g0', '-fopenmp'], } + if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', '-ftemplate-depth-1024'] def build_extensions(self): - """the system will execute the run functions of the base class (distutils.command.build_ext) - and get to this function which is override in the sub class - """ - # First, sanity-check the 'extensions' list, i.e. ext_module - self.check_extensions_list(self.extensions) - for ext in self.extensions: - if isinstance(ext, Extension): - self.regular_compiling(ext) - - def regular_compiling(self, ext): ct = self.compiler.compiler_type opts = self.c_opts.get(ct, []) if ct == 'unix': @@ -120,9 +93,10 @@ def regular_compiling(self, ext): opts.append(cpp_flag(self.compiler)) if has_flag(self.compiler, '-fvisibility=hidden'): opts.append('-fvisibility=hidden') - ext.extra_compile_args = opts - ext.extra_link_args = opts - self.build_extension(ext) + for ext in self.extensions: + ext.extra_compile_args = opts + ext.extra_link_args = opts + build_ext.build_extensions(self) setup( name='brainiak', @@ -145,8 +119,5 @@ def regular_compiling(self, ext): ext_modules=ext_modules, cmdclass={'build_ext': BuildExt}, packages=find_packages(), - package_data = { - 'brainiak.fcma': ['*.pyx', '*.pyxbld'], - }, zip_safe=False, ) From 8f510f74896c4c9e24346defb2c9ecb47e63a155 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 19:27:09 -0700 Subject: [PATCH 44/59] add fcma extension back --- setup.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/setup.py b/setup.py index 8ad6aeb99..b74402663 100644 --- a/setup.py +++ b/setup.py @@ -42,6 +42,15 @@ def __str__(self): get_pybind_include(user=True) ], ), + Extension( + 'brainiak.fcma.fcma_extension', + ['brainiak/fcma/src/fcma_extension.cc'], + include_dirs=[ + # Path to pybind11 headers + get_pybind_include(), + get_pybind_include(user=True) + ], + ), ] From 44225ebd069d6ea15964e27cdc219f8a3e59edff Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Fri, 16 Sep 2016 22:04:12 -0700 Subject: [PATCH 45/59] fortmatting issue --- setup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/setup.py b/setup.py index b74402663..ca3bfeaad 100644 --- a/setup.py +++ b/setup.py @@ -20,6 +20,7 @@ class get_pybind_include(object): """Helper class to determine the pybind11 include path + The purpose of this class is to postpone importing pybind11 until it is actually installed, so that the ``get_include()`` method can be invoked. """ @@ -72,6 +73,7 @@ def has_flag(compiler, flagname): def cpp_flag(compiler): """Return the -std=c++[11/14] compiler flag. + The c++14 is prefered over c++11 (when it is available). """ if has_flag(compiler, '-std=c++14'): From f65f6c859168c4567dc42d40633de1a03ab2e283 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Mon, 19 Sep 2016 13:58:22 -0700 Subject: [PATCH 46/59] copy *pyxbld files when installing --- setup.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/setup.py b/setup.py index ca3bfeaad..21cd010bb 100644 --- a/setup.py +++ b/setup.py @@ -130,5 +130,8 @@ def build_extensions(self): ext_modules=ext_modules, cmdclass={'build_ext': BuildExt}, packages=find_packages(), + package_data = { + 'brainiak.fcma': ['*.pyx', '*.pyxbld'], + }, zip_safe=False, ) From ef554a2a86daaf287bdcaefb4b5d72d119947b74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mihai=20Capot=C4=83?= Date: Tue, 20 Sep 2016 14:50:07 -0700 Subject: [PATCH 47/59] Compile Cython code during install, not at runtime Add more string-like methods required by Cython to pybind11 workaround class. Remove unused `pytest-cython`. Remove `-fvisibility=hidden` which conflicts with Cython. --- .gitignore | 2 ++ brainiak/fcma/__init__.py | 14 -------------- brainiak/fcma/cython_blas.pyxbld | 13 ------------- requirements-dev.txt | 1 - setup.cfg | 4 ---- setup.py | 26 ++++++++++++++++++++++++-- 6 files changed, 26 insertions(+), 34 deletions(-) delete mode 100644 brainiak/fcma/cython_blas.pyxbld diff --git a/.gitignore b/.gitignore index 39739d1ae..afd089b49 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,5 @@ htmlcov .ipynb_checkpoints junit-*.xml __pycache__ + +brainiak/fcma/cython_blas.c diff --git a/brainiak/fcma/__init__.py b/brainiak/fcma/__init__.py index 79b91e2ca..5f9a9900f 100644 --- a/brainiak/fcma/__init__.py +++ b/brainiak/fcma/__init__.py @@ -12,17 +12,3 @@ # See the License for the specific language governing permissions and # limitations under the License. """Full correlation matrix analysis""" - -import pyximport -from mpi4py import MPI -pyximport.install() - -a = 'uninstalled' -# make the cython building execute only once -if MPI.COMM_WORLD.Get_rank() == 0: - from . import cython_blas as blas - blas.installed() - a = 'installed' -# tell everybody the building is done -# functions as a barrier -a = MPI.COMM_WORLD.bcast(a, root=0) diff --git a/brainiak/fcma/cython_blas.pyxbld b/brainiak/fcma/cython_blas.pyxbld deleted file mode 100644 index 45ab31910..000000000 --- a/brainiak/fcma/cython_blas.pyxbld +++ /dev/null @@ -1,13 +0,0 @@ -def make_ext(modname, pyxfilename): - from distutils.extension import Extension - from distutils import sysconfig - import os - opts = [] - if ('CC' in os.environ and 'icc' in os.environ['CC']) or \ - 'icc' in sysconfig.get_config_var('CC'): - opts += ['-lirc', '-lintlc'] - return Extension( - name = modname, - sources = [pyxfilename], - extra_link_args = opts, - ) diff --git a/requirements-dev.txt b/requirements-dev.txt index b7cd60cf5..37097bf57 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,7 +2,6 @@ coverage flake8 flake8-print pytest -pytest-cython restructuredtext-lint sphinx sphinx_rtd_theme diff --git a/setup.cfg b/setup.cfg index eebc7c44b..e095400f4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,3 @@ -[tool:pytest] -addopts = - --doctest-cython - [coverage:run] source = brainiak branch = True diff --git a/setup.py b/setup.py index ca3bfeaad..494c281ff 100644 --- a/setup.py +++ b/setup.py @@ -1,3 +1,5 @@ +from distutils import sysconfig + from setuptools import setup, Extension, find_packages from setuptools.command.build_ext import build_ext import os @@ -28,10 +30,22 @@ class get_pybind_include(object): def __init__(self, user=False): self.user = user + # Required by Cython + def __add__(self, x): + import pybind11 + return pybind11.get_include(self.user) + x + def __str__(self): import pybind11 return pybind11.get_include(self.user) + def endswith(self, x): + import pybind11 + return pybind11.get_include(self.user).endswith(x) + + def startswith(self, x): + import pybind11 + return pybind11.get_include(self.user).startswith(x) ext_modules = [ Extension( @@ -52,6 +66,10 @@ def __str__(self): get_pybind_include(user=True) ], ), + Extension( + 'brainiak.fcma.cython_blas', + ['brainiak/fcma/cython_blas.pyx'], + ), ] @@ -91,6 +109,12 @@ class BuildExt(build_ext): 'unix': ['-g0', '-fopenmp'], } + # FIXME Workaround for using the Intel compiler by setting the CC env var + # Other uses of ICC (e.g., cc binary linked to icc) are not supported + if (('CC' in os.environ and 'icc' in os.environ['CC']) + or 'icc' in sysconfig.get_config_var('CC')): + c_opts += ['-lirc', '-lintlc'] + if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', '-ftemplate-depth-1024'] @@ -102,8 +126,6 @@ def build_extensions(self): opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version()) opts.append(cpp_flag(self.compiler)) - if has_flag(self.compiler, '-fvisibility=hidden'): - opts.append('-fvisibility=hidden') for ext in self.extensions: ext.extra_compile_args = opts ext.extra_link_args = opts From 0b7b078c8c3fff195b5bf2096e64368ab75f07ce Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 18:04:02 -0700 Subject: [PATCH 48/59] only add c++11/14 to C++ extension compiling; check correctness in FCMA tests --- setup.py | 7 ++--- tests/fcma/test_voxel_selection.py | 46 +++++++++++++++++++----------- 2 files changed, 33 insertions(+), 20 deletions(-) diff --git a/setup.py b/setup.py index f958b06b4..30bc7e791 100644 --- a/setup.py +++ b/setup.py @@ -125,10 +125,12 @@ def build_extensions(self): if ct == 'unix': opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version()) - opts.append(cpp_flag(self.compiler)) for ext in self.extensions: ext.extra_compile_args = opts ext.extra_link_args = opts + if ext.language is 'c++': + ext.extra_compile_args.append(cpp_flag(self.compiler)) + ext.extra_link_args.append(cpp_flag(self.compiler)) build_ext.build_extensions(self) setup( @@ -152,8 +154,5 @@ def build_extensions(self): ext_modules=ext_modules, cmdclass={'build_ext': BuildExt}, packages=find_packages(), - package_data = { - 'brainiak.fcma': ['*.pyx', '*.pyxbld'], - }, zip_safe=False, ) diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index ccc72c04e..2da8e701b 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -18,13 +18,15 @@ import numpy as np import math from mpi4py import MPI -import time +from numpy.random import RandomState +# specify the random state to fix the random numbers +prng = RandomState(1234567890) def create_epoch(): row = 12 - col = 1000 - mat = np.random.rand(row, col).astype(np.float32) + col = 5 + mat = prng.rand(row, col).astype(np.float32) mat = zscore(mat, axis=0, ddof=0) # if zscore fails (standard deviation is zero), # set all values to be zero @@ -34,21 +36,33 @@ def create_epoch(): def test_voxel_selection(): fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), create_epoch(), create_epoch()] - labels = [0, 1, 0, 1] - # set master rank to be 0 - vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=0) + labels = [0, 1, 0, 1, 0, 1, 0, 1] + # 2 subjects, 4 epochs per subject + vs = VoxelSelector(fake_raw_data, 4, labels, 2) # for cross validation, use SVM with precomputed kernel - # no shrinking, set C=10 - clf = svm.SVC(kernel='precomputed', shrinking=False, C=10) - results0 = vs.run(clf) - # set master rank to be 1 and do it again - vs = VoxelSelector(fake_raw_data, 2, labels, 2, master_rank=1) - results1 = vs.run(clf) + # no shrinking, set C=1 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + results = vs.run(clf) + if MPI.COMM_WORLD.Get_rank() == 0: + expected_results = [(0, 0.75), (2, 0.625), (1, 0.5), + (3, 0.5), (4, 0.375)] + assert np.array_equal(results, expected_results), \ + 'voxel selection does not provide correct results' # test scipy normalization - fake_corr = np.random.rand(1, 12, 100).astype(np.float32) + fake_corr = prng.rand(1, 4, 5).astype(np.float32) fake_corr = vs._correlationNormalization(fake_corr) - # make one process sleep a while to resolve file writing competition if MPI.COMM_WORLD.Get_rank() == 0: - time.sleep(0.5) - return results0, results1, fake_corr + expected_fake_corr = [[[1.19203866, 0.18862808, -0.54350245, + -1.18334889, -0.16860008], + [-1.06594729, -1.08742261, -1.19447124, + 1.14114654, -0.67860204], + [0.7839641, 1.53981364, 0.24948341, + 0.82626557, 1.67902875], + [-0.91005552, -0.64101928, 1.48848987, + -0.78406328, -0.83182675]]] + assert np.allclose(fake_corr, expected_fake_corr), \ + 'within-subject normalization does not provide correct results' + return results, fake_corr \ No newline at end of file From 5eadef11e9b3aa23149e7a2ea2da3ee043604ea8 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 18:11:37 -0700 Subject: [PATCH 49/59] explicitly specify language for each extension --- setup.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/setup.py b/setup.py index 30bc7e791..bff716a2c 100644 --- a/setup.py +++ b/setup.py @@ -56,6 +56,7 @@ def startswith(self, x): get_pybind_include(), get_pybind_include(user=True) ], + language='c++', ), Extension( 'brainiak.fcma.fcma_extension', @@ -65,10 +66,12 @@ def startswith(self, x): get_pybind_include(), get_pybind_include(user=True) ], + language='c++', ), Extension( 'brainiak.fcma.cython_blas', ['brainiak/fcma/cython_blas.pyx'], + language='c', ), ] From 39f9f079e8f553b65164e4d4e939e7fcc3e47284 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 18:23:43 -0700 Subject: [PATCH 50/59] use == to check if two strings equal --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index bff716a2c..d7d60a3fb 100644 --- a/setup.py +++ b/setup.py @@ -131,7 +131,7 @@ def build_extensions(self): for ext in self.extensions: ext.extra_compile_args = opts ext.extra_link_args = opts - if ext.language is 'c++': + if ext.language == 'c++': ext.extra_compile_args.append(cpp_flag(self.compiler)) ext.extra_link_args.append(cpp_flag(self.compiler)) build_ext.build_extensions(self) From 031d5fb1bb4c02aeb3c55ce075f6b2e4e26e25e1 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 18:33:41 -0700 Subject: [PATCH 51/59] set icc flags to c_opts correctly --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index d7d60a3fb..48f90a0c4 100644 --- a/setup.py +++ b/setup.py @@ -116,7 +116,7 @@ class BuildExt(build_ext): # Other uses of ICC (e.g., cc binary linked to icc) are not supported if (('CC' in os.environ and 'icc' in os.environ['CC']) or 'icc' in sysconfig.get_config_var('CC')): - c_opts += ['-lirc', '-lintlc'] + c_opts['unix'] += ['-lirc', '-lintlc'] if sys.platform == 'darwin': c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7', From 3d86e86a0f9319e4e374745195ccc0ec51f8291b Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 19:00:41 -0700 Subject: [PATCH 52/59] detect extension language in a smarter way; let fcma test code read in data files --- setup.py | 6 ++---- tests/fcma/corr.npytest | Bin 0 -> 160 bytes tests/fcma/raw_data.npytest | Bin 0 -> 2000 bytes tests/fcma/test_voxel_selection.py | 17 +++++++++++------ 4 files changed, 13 insertions(+), 10 deletions(-) create mode 100644 tests/fcma/corr.npytest create mode 100644 tests/fcma/raw_data.npytest diff --git a/setup.py b/setup.py index 48f90a0c4..6a7902c5e 100644 --- a/setup.py +++ b/setup.py @@ -56,7 +56,6 @@ def startswith(self, x): get_pybind_include(), get_pybind_include(user=True) ], - language='c++', ), Extension( 'brainiak.fcma.fcma_extension', @@ -66,12 +65,10 @@ def startswith(self, x): get_pybind_include(), get_pybind_include(user=True) ], - language='c++', ), Extension( 'brainiak.fcma.cython_blas', ['brainiak/fcma/cython_blas.pyx'], - language='c', ), ] @@ -131,7 +128,8 @@ def build_extensions(self): for ext in self.extensions: ext.extra_compile_args = opts ext.extra_link_args = opts - if ext.language == 'c++': + lang = ext.language or self.compiler.detect_language(ext.sources) + if lang == 'c++': ext.extra_compile_args.append(cpp_flag(self.compiler)) ext.extra_link_args.append(cpp_flag(self.compiler)) build_ext.build_extensions(self) diff --git a/tests/fcma/corr.npytest b/tests/fcma/corr.npytest new file mode 100644 index 0000000000000000000000000000000000000000..dc4acb76012190e9fcb0eb3c2a888953666e9230 GIT binary patch literal 160 zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlV+l>qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$-ItnH_3Z|Mm3bhIlz;*U%s{M^Tcl)KAIPDi+pKCWW(%Jr1K(1YyP^H}} z)mL`A=6Kp)`7L1o;A4XQAHxLueTp*ncg{Q8PY%hnS6}?Xu1QkVe$zXBd%;P6?6khz GwgUi1S2`a6 literal 0 HcmV?d00001 diff --git a/tests/fcma/raw_data.npytest b/tests/fcma/raw_data.npytest new file mode 100644 index 0000000000000000000000000000000000000000..552d5ff9ea3f6e8f0deabd8aa0ee329d878756ac GIT binary patch literal 2000 zcmW;Ac{~(&AHeY$HbN&|TU$g-j<#dzGQaO^WfdE#hmNhu=1T6y*(xNLF&l%1 zs2t6fW`5rxD$-U+YNfJ8#g@dpqG#yYeLnwv-dT<=4o(M{3|q!YmWOAsJCC*5fW`K* zWLX%nyaIV4ydwdwfjke-f8F+oU$Eyt9_)Q2$n&3HxzWO4?K%qsD^m-DD1-l}ed=ow zS*oohbFI?BDj^3RZ$3y<-=u^8ZVwE4bRT>fk+DX?T*Njjc*K;t7_#n#F+h?|v2oj(n*Z`9im|Hv_U`HXm! zx(IR8)MABaaxupI;R5OJZi6!8B|27ik_`uwVAo&@lqToVxuW&Z?-U5jnO|U|z?XEN z5@Dw>6KX@!AkBM_n(JpmK1YaWE84(#A_@45biq0@7K18EpvWc{_W$Hdn;pa7Vr)cXk6P0z9UXBqap!VMoIY0QbynI<3i=%H(?GV z04)kyP|bFLZgD80F(>@7>6MC%{E`YP2R^*psKYVebrmO5?|@2p4ns?1)PAT8rM_+W zux2{Q*M-9KyL_f`_pGeU#r`+%-R z57l^RN_y`Ng0SN$xa5tKaIM*B{C*h!-Xp^2%n@vMI*ZoJs^ECaHmIIoj$Ho~7`!FJ zf>Q~gF5Zc)yjUVF>c;-Hx}awxKuJs{vLj-_r|BR%_())QQ!nlbD~Fb8pP=qU37m1w z$Lk^45U%?X3YLe{_I^jS8_S^I9P}oxT>-#tzKhY%uVNq=QmZh1TvTj|Qr8?hRx<^^ z_6*_D<(i6+HVKy6r<3SEn`qgxKxhycqKS4ZZONZNfAnFN-uWUSdVmVw2&BA_NAEmtdjADHup?z~1C?;H{YpDQ1N@S$GSC z>8s&>-WV!9YJhQRSe|~NRQ{nx7e`icAKWRnW_%j zm2`!U-JOmB<_ymFKfi^x?Kg4p)vKr)F2orZUs6$D2|KIjk$mf)aZvP&0#to`yfAeU z{cnb#plu;%GUp9-|DB5g!9B38Oi16KR>xIg#TXjliyGPSD9GCaN2L$&Yt=`n6qUj4 zOfL+dvlEWIxWXdA*=oSb<8nYq0ZP1l}wCjFk(f z!P0mNEuZVrwK53DQxn)dww8F^k`3$TB@p-UcDOnE51IVlo@Qj4T9>uJew?i7Xhl12K58C9 zkn}xAU7u&*<9Ath;LvVTVgC*%a^v85xCpe@zNTkFhiEa|6THvG<9q8NS_{QsmcvF9 zZW{{MC(@?9KpW#@pebk=`#;x_y1kjOWYSa-WEl?@XA$0CYG+6CE>rDKTj5LF3i`pg z0cRe40?}@Fh}gUY7wX)C*KQX)^L5+4z9?upg}B!|A?em{t1|=pz&j1dKg(Q>L2>+QIB&e;xA6dbxu zX8~R@KY&jz3qfgq1_D3EQSo6Z+6?;Oi+f#iqx$n;xLSlie(=NOs+n+{*$xBF0{Qxe z4Y0rc4JaqH6t^2?G^un4*t3MNV*hHAGs&g;Vi|sX(1f!t2I9QbQF?#LY@%kIWF1wN zOJ4mD0(|c%+#c{3e6u?lZFslwu-hUDGlHG5%ApQBbkaa|IR(15 Date: Tue, 20 Sep 2016 20:14:11 -0700 Subject: [PATCH 53/59] use deep copy to assign compiling flags --- setup.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 6a7902c5e..a0eb1b366 100644 --- a/setup.py +++ b/setup.py @@ -5,6 +5,7 @@ import os import sys import setuptools +from copy import deepcopy __version__ = '0.2' @@ -126,8 +127,8 @@ def build_extensions(self): opts.append('-DVERSION_INFO="%s"' % self.distribution.get_version()) for ext in self.extensions: - ext.extra_compile_args = opts - ext.extra_link_args = opts + ext.extra_compile_args = deepcopy(opts) + ext.extra_link_args = deepcopy(opts) lang = ext.language or self.compiler.detect_language(ext.sources) if lang == 'c++': ext.extra_compile_args.append(cpp_flag(self.compiler)) From f75d9f76c803d7939c00d0009f065005d596de81 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 20:28:30 -0700 Subject: [PATCH 54/59] switch the test sequence --- tests/fcma/test_voxel_selection.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index bd5f72563..2ac0ac99d 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -43,15 +43,6 @@ def test_voxel_selection(): labels = [0, 1, 0, 1, 0, 1, 0, 1] # 2 subjects, 4 epochs per subject vs = VoxelSelector(fake_raw_data, 4, labels, 2) - # for cross validation, use SVM with precomputed kernel - # no shrinking, set C=1 - clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) - results = vs.run(clf) - if MPI.COMM_WORLD.Get_rank() == 0: - expected_results = [(0, 0.75), (2, 0.625), (1, 0.5), - (3, 0.5), (4, 0.375)] - assert np.array_equal(results, expected_results), \ - 'voxel selection does not provide correct results' # test scipy normalization #fake_corr = prng.rand(1, 4, 5).astype(np.float32) fake_corr = np.load('tests/fcma/corr.npytest') @@ -67,6 +58,15 @@ def test_voxel_selection(): -0.78406328, -0.83182675]]] assert np.allclose(fake_corr, expected_fake_corr), \ 'within-subject normalization does not provide correct results' + # for cross validation, use SVM with precomputed kernel + # no shrinking, set C=1 + clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) + results = vs.run(clf) + if MPI.COMM_WORLD.Get_rank() == 0: + expected_results = [(0, 0.75), (2, 0.625), (1, 0.5), + (3, 0.5), (4, 0.375)] + assert np.array_equal(results, expected_results), \ + 'voxel selection does not provide correct results' return results, fake_corr if __name__ == '__main__': From 257882d91415dfb3ec39b5e8191d2a2b04a21504 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 21:04:57 -0700 Subject: [PATCH 55/59] make kernel matrix scale setting adjustable to the initial values; remove file reading in fcma test --- brainiak/fcma/cython_blas.pyx | 5 ++++- tests/fcma/test_voxel_selection.py | 16 +++++++--------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/brainiak/fcma/cython_blas.pyx b/brainiak/fcma/cython_blas.pyx index 84ccb6925..6f4e2ad3f 100644 --- a/brainiak/fcma/cython_blas.pyx +++ b/brainiak/fcma/cython_blas.pyx @@ -183,7 +183,10 @@ def compute_kernel_matrix(py_uplo, py_trans, py_n, py_k, py_alpha, py_a, int py_ 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 - py_c *= .001 + 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/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index 2ac0ac99d..cecef9f68 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -35,17 +35,15 @@ def create_epoch(): return mat def test_voxel_selection(): - #fake_raw_data = [create_epoch(), create_epoch(), - # create_epoch(), create_epoch(), - # create_epoch(), create_epoch(), - # create_epoch(), create_epoch()] - fake_raw_data = np.load('tests/fcma/raw_data.npytest') + fake_raw_data = [create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch(), + create_epoch(), create_epoch()] labels = [0, 1, 0, 1, 0, 1, 0, 1] # 2 subjects, 4 epochs per subject vs = VoxelSelector(fake_raw_data, 4, labels, 2) # test scipy normalization - #fake_corr = prng.rand(1, 4, 5).astype(np.float32) - fake_corr = np.load('tests/fcma/corr.npytest') + fake_corr = prng.rand(1, 4, 5).astype(np.float32) fake_corr = vs._correlationNormalization(fake_corr) if MPI.COMM_WORLD.Get_rank() == 0: expected_fake_corr = [[[1.19203866, 0.18862808, -0.54350245, @@ -63,8 +61,8 @@ def test_voxel_selection(): clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) results = vs.run(clf) if MPI.COMM_WORLD.Get_rank() == 0: - expected_results = [(0, 0.75), (2, 0.625), (1, 0.5), - (3, 0.5), (4, 0.375)] + expected_results = [(0, 0.875), (2, 0.75), (1, 0.5), + (3, 0.5), (4, 0.5)] assert np.array_equal(results, expected_results), \ 'voxel selection does not provide correct results' return results, fake_corr From 3f14b8062cdf31cf19dda3564fbadbf4f25f0a69 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 21:21:57 -0700 Subject: [PATCH 56/59] tolerate one sample error in cross validation --- tests/fcma/test_voxel_selection.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index cecef9f68..54e182d10 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -61,9 +61,11 @@ def test_voxel_selection(): clf = svm.SVC(kernel='precomputed', shrinking=False, C=1) results = vs.run(clf) if MPI.COMM_WORLD.Get_rank() == 0: - expected_results = [(0, 0.875), (2, 0.75), (1, 0.5), - (3, 0.5), (4, 0.5)] - assert np.array_equal(results, expected_results), \ + output = [None] * len(results) + for tuple in results: + output[tuple[0]] = int(8*tuple[1]) + expected_output = [7, 4, 6, 4, 4] + assert np.allclose(output, expected_output, atol=1), \ 'voxel selection does not provide correct results' return results, fake_corr From 9a3f00bd74044a4c34a7979cad751d8c18ec3762 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Tue, 20 Sep 2016 21:34:10 -0700 Subject: [PATCH 57/59] delete unnecessary data file; improve fcma test coverage --- tests/fcma/corr.npytest | Bin 160 -> 0 bytes tests/fcma/raw_data.npytest | Bin 2000 -> 0 bytes tests/fcma/test_voxel_selection.py | 2 +- 3 files changed, 1 insertion(+), 1 deletion(-) delete mode 100644 tests/fcma/corr.npytest delete mode 100644 tests/fcma/raw_data.npytest diff --git a/tests/fcma/corr.npytest b/tests/fcma/corr.npytest deleted file mode 100644 index dc4acb76012190e9fcb0eb3c2a888953666e9230..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 160 zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlV+l>qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$-ItnH_3Z|Mm3bhIlz;*U%s{M^Tcl)KAIPDi+pKCWW(%Jr1K(1YyP^H}} z)mL`A=6Kp)`7L1o;A4XQAHxLueTp*ncg{Q8PY%hnS6}?Xu1QkVe$zXBd%;P6?6khz GwgUi1S2`a6 diff --git a/tests/fcma/raw_data.npytest b/tests/fcma/raw_data.npytest deleted file mode 100644 index 552d5ff9ea3f6e8f0deabd8aa0ee329d878756ac..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2000 zcmW;Ac{~(&AHeY$HbN&|TU$g-j<#dzGQaO^WfdE#hmNhu=1T6y*(xNLF&l%1 zs2t6fW`5rxD$-U+YNfJ8#g@dpqG#yYeLnwv-dT<=4o(M{3|q!YmWOAsJCC*5fW`K* zWLX%nyaIV4ydwdwfjke-f8F+oU$Eyt9_)Q2$n&3HxzWO4?K%qsD^m-DD1-l}ed=ow zS*oohbFI?BDj^3RZ$3y<-=u^8ZVwE4bRT>fk+DX?T*Njjc*K;t7_#n#F+h?|v2oj(n*Z`9im|Hv_U`HXm! zx(IR8)MABaaxupI;R5OJZi6!8B|27ik_`uwVAo&@lqToVxuW&Z?-U5jnO|U|z?XEN z5@Dw>6KX@!AkBM_n(JpmK1YaWE84(#A_@45biq0@7K18EpvWc{_W$Hdn;pa7Vr)cXk6P0z9UXBqap!VMoIY0QbynI<3i=%H(?GV z04)kyP|bFLZgD80F(>@7>6MC%{E`YP2R^*psKYVebrmO5?|@2p4ns?1)PAT8rM_+W zux2{Q*M-9KyL_f`_pGeU#r`+%-R z57l^RN_y`Ng0SN$xa5tKaIM*B{C*h!-Xp^2%n@vMI*ZoJs^ECaHmIIoj$Ho~7`!FJ zf>Q~gF5Zc)yjUVF>c;-Hx}awxKuJs{vLj-_r|BR%_())QQ!nlbD~Fb8pP=qU37m1w z$Lk^45U%?X3YLe{_I^jS8_S^I9P}oxT>-#tzKhY%uVNq=QmZh1TvTj|Qr8?hRx<^^ z_6*_D<(i6+HVKy6r<3SEn`qgxKxhycqKS4ZZONZNfAnFN-uWUSdVmVw2&BA_NAEmtdjADHup?z~1C?;H{YpDQ1N@S$GSC z>8s&>-WV!9YJhQRSe|~NRQ{nx7e`icAKWRnW_%j zm2`!U-JOmB<_ymFKfi^x?Kg4p)vKr)F2orZUs6$D2|KIjk$mf)aZvP&0#to`yfAeU z{cnb#plu;%GUp9-|DB5g!9B38Oi16KR>xIg#TXjliyGPSD9GCaN2L$&Yt=`n6qUj4 zOfL+dvlEWIxWXdA*=oSb<8nYq0ZP1l}wCjFk(f z!P0mNEuZVrwK53DQxn)dww8F^k`3$TB@p-UcDOnE51IVlo@Qj4T9>uJew?i7Xhl12K58C9 zkn}xAU7u&*<9Ath;LvVTVgC*%a^v85xCpe@zNTkFhiEa|6THvG<9q8NS_{QsmcvF9 zZW{{MC(@?9KpW#@pebk=`#;x_y1kjOWYSa-WEl?@XA$0CYG+6CE>rDKTj5LF3i`pg z0cRe40?}@Fh}gUY7wX)C*KQX)^L5+4z9?upg}B!|A?em{t1|=pz&j1dKg(Q>L2>+QIB&e;xA6dbxu zX8~R@KY&jz3qfgq1_D3EQSo6Z+6?;Oi+f#iqx$n;xLSlie(=NOs+n+{*$xBF0{Qxe z4Y0rc4JaqH6t^2?G^un4*t3MNV*hHAGs&g;Vi|sX(1f!t2I9QbQF?#LY@%kIWF1wN zOJ4mD0(|c%+#c{3e6u?lZFslwu-hUDGlHG5%ApQBbkaa|IR(15 Date: Tue, 20 Sep 2016 23:18:47 -0700 Subject: [PATCH 58/59] add logistic regression to fcma test --- tests/fcma/test_voxel_selection.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tests/fcma/test_voxel_selection.py b/tests/fcma/test_voxel_selection.py index 263dea497..cd9d1c6c6 100644 --- a/tests/fcma/test_voxel_selection.py +++ b/tests/fcma/test_voxel_selection.py @@ -15,6 +15,7 @@ from brainiak.fcma.voxelselector import VoxelSelector from scipy.stats.mstats import zscore from sklearn import svm +from sklearn.linear_model import LogisticRegression import numpy as np import math from mpi4py import MPI @@ -66,8 +67,17 @@ def test_voxel_selection(): output[tuple[0]] = int(8*tuple[1]) expected_output = [7, 4, 6, 4, 4] assert np.allclose(output, expected_output, atol=1), \ - 'voxel selection does not provide correct results' - return results, fake_corr + 'voxel selection via SVM does not provide correct results' + # for cross validation, use logistic regression + clf = LogisticRegression() + results = vs.run(clf) + if MPI.COMM_WORLD.Get_rank() == 0: + output = [None] * len(results) + for tuple in results: + output[tuple[0]] = int(8*tuple[1]) + expected_output = [6, 3, 6, 4, 4] + assert np.allclose(output, expected_output, atol=1), \ + 'voxel selection via SVM does not provide correct results' if __name__ == '__main__': - test_voxel_selection() \ No newline at end of file + test_voxel_selection() From e4823b18486f4cc4025f65c710dd88536ff5e2c1 Mon Sep 17 00:00:00 2001 From: Yida Wang Date: Wed, 21 Sep 2016 09:49:44 -0700 Subject: [PATCH 59/59] rename the code file in fcma example --- examples/fcma/{run.py => voxel_selection.py} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename examples/fcma/{run.py => voxel_selection.py} (98%) diff --git a/examples/fcma/run.py b/examples/fcma/voxel_selection.py similarity index 98% rename from examples/fcma/run.py rename to examples/fcma/voxel_selection.py index 3ced3c8a0..afcf29963 100644 --- a/examples/fcma/run.py +++ b/examples/fcma/voxel_selection.py @@ -188,7 +188,7 @@ def prepareData(data_dir, extension, mask_file, epoch_file): """ example running command: -mpirun -np 2 python run.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz +mpirun -np 2 python voxel_selection.py /Users/yidawang/data/face_scene/raw nii.gz /Users/yidawang/data/face_scene/mask.nii.gz data/fs_epoch_labels.npy 12 18 """ if __name__ == '__main__':