diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..8fae155 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required (VERSION 2.8) + +include(cmake/LookUp-GreatCMakeCookOff.cmake) +set(LIBRARY_OUTPUT_PATH ${CMAKE_SOURCE_DIR}/build/lib) +set(EXECUTABLE_OUTPUT_PATH ${CMAKE_SOURCE_DIR}/tests) +string(ASCII 27 Esc) + +option(double "Enable testing" off) + +if(double) + add_definitions(-DDOUBLE__PRECISION) +endif() + + +set(ColourReset "${Esc}[m") +set(Red "${Esc}[1;31m") +message("${Red}export PYTHONPATH=$PYTHONPATH:${LIBRARY_OUTPUT_PATH} ${ColourReset}") + +find_package(Numpy REQUIRED) +find_package(CoherentPython REQUIRED) +find_package(CUDA QUIET REQUIRED) + +find_library(FOUND_CUBLAS cublas) +include_directories(${PYTHON_INCLUDE_DIRS}) +include_directories(${NUMPY_INCLUDE_DIRS}/numpy) + +set(CUBLAS_TARGET_LINK ${FOUND_CUBLAS}) + +subdirs(gp_emulator/gpu) diff --git a/MANIFEST.in b/MANIFEST.in index 4cad124..61f6ad6 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,3 @@ include setup.py README.rst -recursive-include gp_emulator/GaussianProcess.py __init__.py lhd.py multivariate_gp.py save_emulators.py +#recursive-include gp_emulator/ GaussianProcess.py __init__.py lhd.py multivariate_gp.py save_emulators.py +#include build/lib/_gpu_predict.so diff --git a/README.rst b/README.rst index d640933..bdd96bf 100644 --- a/README.rst +++ b/README.rst @@ -8,8 +8,38 @@ GP emulators This repository contains an implementation of GPs for emulation in Python. Although many different implementations exist, this particular one deals with fast GP predictions for large number of input vectors, where the training data sets are typically modest (e.g. less than 300 samples). Access to the emulation's partial derivatives and Hessian matrix is calculated, and training is also taken care of. -You can install with +Requirements: +--------- +* python ( 2.7 or later ) +* numpy +* scipy +* GPU predict module: + * cmake + * CUDA 5.0 or later + * CUnit - python setup.py install +Install with GPU predict: +--------- +1. decide precision by modifying CMakeList.txt + `add_definition(-DDOUBLE_PRECISION) # double` +2. install + `python setup.py install` + +Tests (for GPU predict): +---------- -The only requirements are (if memory serves) numpy and scipy. +test: `python setup.py test` + +benchmark:`python setup.py benchmark` + +1. Unit testing ([unit_tests.py](https://github.com/UCL/gp_emulator/blob/master/tests/unit_tests.py)): + * random inputs generated by python predict + * operate unit testings of GPU functions. + * compare GPU and python outputs +2. benchmark [benchmark.py](https://github.com/UCL/gp_emulator/blob/master/tests/benchmark.py) + * obtain speedup of GPU predict + * random inputs +3. testing emulator ([testing_emulator.py](https://github.com/UCL/gp_emulator/blob/master/tests/test_perband_emulator.py)) + * read real data from .zpn file + * run emulator with and without GPU + diff --git a/build/.gitignore b/build/.gitignore new file mode 100644 index 0000000..5e7d273 --- /dev/null +++ b/build/.gitignore @@ -0,0 +1,4 @@ +# Ignore everything in this directory +* +# Except this file +!.gitignore diff --git a/cmake/LookUp-GreatCMakeCookOff.cmake b/cmake/LookUp-GreatCMakeCookOff.cmake new file mode 100644 index 0000000..fae65ca --- /dev/null +++ b/cmake/LookUp-GreatCMakeCookOff.cmake @@ -0,0 +1,59 @@ +# This file is part of the GreatCMakeCookOff package and distributed under the MIT Licences. +# Upon inclusion in a cmake file, it will download the GreatCMakeCookOff and make itself known to +# CMake. It should be added explicitely to build systems that make use of recipes from the cook-off. +# And it should be included prior to using cook-off recipes: +# +# ```{CMake} +# include(LookUp-GreatCMakeCookOff) +# ``` + + +# First attempts to find the package +set(COOKOFF_DOWNLOAD_DIR "${PROJECT_BINARY_DIR}/external/src/GreatCMakeCookOff") +find_package(GreatCMakeCookOff NO_MODULE PATHS "${COOKOFF_DOWNLOAD_DIR}" QUIET) + +# Otherwise attempts to download it. +# Does not use ExternalProject_Add to avoid doing a recursive cmake step. +if(NOT GreatCMakeCookOff_FOUND) + message(STATUS "[GreatCMakeCookOff] not found. Will attempt to clone it.") + + # Need git for cloning. + find_package(Git) + if(NOT GIT_FOUND) + message(FATAL_ERROR "[Git] not found. Cannot download GreatCMakeCookOff") + endif() + + # Remove GreatCMakeCookOff directory if it exists + if(EXISTS "${COOKOFF_DOWNLOAD_DIR}") + execute_process( + COMMAND ${CMAKE_COMMAND} -E remove_directory "${COOKOFF_DOWNLOAD_DIR}" + OUTPUT_QUIET + ) + endif() + if(NOT COOKOFF_GITREPO) + set(COOKOFF_GITREPO https://github.com/UCL/GreatCMakeCookOff.git) + endif() + execute_process( + COMMAND ${GIT_EXECUTABLE} clone "${COOKOFF_GITREPO}" + "${COOKOFF_DOWNLOAD_DIR}" + RESULT_VARIABLE CLONING_COOKOFF + OUTPUT_QUIET + ERROR_VARIABLE CLONING_ERROR + ) + + if(NOT ${CLONING_COOKOFF} EQUAL 0) + message(STATUS "${CLONING_ERROR}") + message(FATAL_ERROR "[GreatCMakeCookOff] git cloning failed.") + else() + message(STATUS "[GreatCMakeCookOff] downloaded to ${COOKOFF_DOWNLOAD_DIR}") + find_package(GreatCMakeCookOff NO_MODULE PATHS "${COOKOFF_DOWNLOAD_DIR}" QUIET) + endif() + + set(GreatCMakeCookOff_DIR "${COOKOFF_DOWNLOAD_DIR}/cmake") + set(GreatCMakeCookOff_FOUND TRUE) +endif() +unset(COOKOFF_DOWNLOAD_DIR) + +# Adds GreatCMakeCookOff to module paths +initialize_cookoff() + diff --git a/data/.gitignore b/data/.gitignore new file mode 100644 index 0000000..5e7d273 --- /dev/null +++ b/data/.gitignore @@ -0,0 +1,4 @@ +# Ignore everything in this directory +* +# Except this file +!.gitignore diff --git a/doc/Gomez_Lewis.pdf b/doc/Gomez_Lewis.pdf new file mode 100644 index 0000000..ace2bae Binary files /dev/null and b/doc/Gomez_Lewis.pdf differ diff --git a/doc/build.sh b/doc/build.sh new file mode 100755 index 0000000..e24c8de --- /dev/null +++ b/doc/build.sh @@ -0,0 +1 @@ +pandoc report.md -o report.pdf -H report.latex diff --git a/doc/kernel.png b/doc/kernel.png new file mode 100644 index 0000000..ac0094b Binary files /dev/null and b/doc/kernel.png differ diff --git a/doc/report.latex b/doc/report.latex new file mode 100644 index 0000000..daa08fc --- /dev/null +++ b/doc/report.latex @@ -0,0 +1,48 @@ +\let\oldlongtable\longtable +\def\longtable{\tiny \oldlongtable} +\usepackage[svgnames]{xcolor} +\usepackage{libertine} +\usepackage{framed} +\let\oldquote\quote +\let\endoldquote\endquote +\definecolor{verypalegreen}{RGB}{252,255,250} +\colorlet{shadecolor}{verypalegreen} + +\makeatletter +\def\shadequote{\begin{snugshade}\begin{oldquote}} +\def\endshadequote{% + \end{oldquote}\end{snugshade}} +\makeatother +\renewenvironment{quote}{\begin{shadequote}}{\end{shadequote}} +\usepackage{lettrine} % Package to accentuate the first letter of the text +\usepackage{fix-cm} % Custom font sizes - used for the initial letter in the document + +\definecolor{uclmidgreen}{RGB}{130,141,55} + +\newcommand{\initial}[1]{ % Defines the command and style for the first letter +\lettrine[lines=3,lhang=0.3,nindent=0em]{ +\color{DarkGreen} +{\textsf{#1}}}{}} + +\usepackage{sectsty} % Enables custom section titles +\sectionfont{\color{uclmidgreen} \usefont{OT1}{phv}{b}{n}} % Change the font of all section commands +\subsectionfont{\color{DarkSeaGreen} \usefont{OT1}{phv}{b}{n}} % Change the font of all section commands + +\usepackage{titling} % Allows custom title configuration +\newcommand{\HorRule}{\color{DarkSeaGreen} \rule{\linewidth}{1pt}} % Defines the gold horizontal rule around the title + +\pretitle{ +\vspace{-30pt} \begin{flushleft} \HorRule \fontsize{38}{38} \usefont{OT1}{phv}{b}{n} \color{uclmidgreen} \selectfont} % Horizontal rule before the title + +\posttitle{\par\end{flushleft}\begin{flushleft}\fontsize{25}{25} \usefont{OT1}{phv}{b}{n} \color{Black} GPU acceleration of Gaussian process emulator\selectfont \end{flushleft}\vskip 0.5em } % Whitespace under the title and subtitle + +\preauthor{\begin{flushleft}\large \lineskip 0.5em \usefont{OT1}{phv}{b}{sl} \color{uclmidgreen} Sinan Shi}% Author font configuration + +\postauthor{\footnotesize \usefont{OT1}{phv}{m}{sl} \color{Black} % Configuration for the institution name + + +University College London% Your institution + +\par\end{flushleft}\HorRule} % Horizontal rule after the title + +\hypersetup{linkcolor=black} diff --git a/doc/report.md b/doc/report.md new file mode 100644 index 0000000..5dfba6e --- /dev/null +++ b/doc/report.md @@ -0,0 +1,95 @@ +% Project Report +# Introduction +Land surface coverage data from satellite observations must be interpreted using physical models, which are computationally expansive. Using a surrogate Gaussian process (GP) emulator instead of complex physical models has been proved much more efficient. The GP emulator is written in python, which is part of a python-based Earth Observation Land Data Assimilation System (EO-LDAS). However the current performance of the GP emulator is still far away from ready to process vast data from the satellites. + +The performance bottleneck of the emulator is the _predict_ function. The objective of this project is to provide an appropriate parallel solution for this predict function based on an investigate the performance of both CPU and GPU platforms. A prototype GPU CUDA/C++ implementation of the _predict_ function has been integrated to the python emulator, which has achieved 65x speedup on K20 GPU comparing to the python CPU implementation. + + +# Python predict function +## Code explanation + +The emulator calculates a hyper-parameters set first and uses these hyper-parameters to calculate the approximation of the parameters, uncertainty and the associated gradient by the _predict_ function. The _predict_ function has two inputs, first is the 2D array consists of Nparam (~10) parameter vectors (Npredict), such as biophysical parameters or state variables (e.g leaf area index) by satellites, and second the inverted hyper-parameter covariance array (Ntrain x Ntrain)(~250) provided by the emulator. The length of the input vectors can reach 1e8. + +The python _predict_ function is wrapped in GaussianProcess class, which has around 30 lines with six major numpy or scipy functions/operations of either matrix operations or array element wise operations. There are two features of _predict_ that can benefit the GPU implementation, which is first, there is no branching; second, inside each parameter vector (~ 1e8), elements can be calculated independently, i.e. we can arbitrarily decide the size based on the performance tuning and the capability of GPU. We will go through this in the later sections. + +## Performance analysis +Tests have been carried out on Legion with Xeon 5500 CPU. Interestingly, according to the performance table below, the most computational intensive functions are not the most time consuming ones, e.g. the matrix multiplication costs the same amount of time per call as the element wise array multiplication where the previous one has one magnitude higher the computational complexity than the later one. This inconsistency is caused by a different level of optimisation of numpy/scipy. If we take a close look at the numpy source code, its matrix multiplication function has a C extension with CBLAS. That’s the reason why the performance numpy matrix operations is almost comparable to the C++ implementation. The element wise operations take 88% of total CPU wall time. On the contrary, due to a good optimisation under numpy and less number of calling, the matrix operations take only 12%. + +| Rank | Operation | Call | Percentage of total wall time (%) | +|:----: |-------------------------------------------------- |:----: |:---------------------------------: | +| 1 | Array element wise multiplication (Matmul - Ele) | 10 | 58.01 | +| 2 | Matrix vector subtraction (MvSub) | 10 | 15.08 | +| 3 | Array update: exp (Exp) | 3 | 9.71 | +| 4 | Matrix multiplication | 1 | 8.47 | +| 5 | Euclidean distance (Dist) | 1 | 5.54 | +| 6 | Matrix vector multiplication (Mv) | 10 | 3.00 | +| | Total | | 99.81 | + +## Comparison with C-BLAS implementation +![Kernel Performance with Vector Size 2e6](kernel.png) + +In general, C/C++ is regarded as one of the most efficient languages on CPU platforms. Before implementing the GPU _predict_, we would like to know how much potential does a CPU has on C/C++ with reasonable amount of effort. A very primitive C++ with CBLAS implementation has been carried out for this purpose. + +All major functions except distance calculation have a performance improvement. The performance of the element wise array multiplication has a significant improvement. However what we didn’t expect is that the distance function of C++ is even slower the python scipy _cdist_. +The reason of that is maybe the bad memory access pattern obtained in the our distance function. We will just keep it in mind, and not do any further optimisation on it since this implementation is just to give us an idea about how fast can CPU achieve. The scipy distance performance can more or less represent the highest level of CPU performance with reasonable amount of effort. + +The result of the performance comparison shows the C++ predict function can be 1.3x faster than its python counterpart. If we ignore impact of distance function, C++ implementation is around 2.5x – 5x faster than original python implementation. The performance can vary from platforms, compilers and python versions. + + +# GPU implementation of GP-Emulator +## CUDA/C++ implementation +CUDA with C/C++ extension is one of the most common and straightforward ways of heterogeneous programming on CUDA (NVidia) device. The CUDA APIs enable programs to control the CUDA device and memory explicitly. + +In spite that there is no standard python CUDA extension provided by NVidia, there are many third party APIs or wrappers that allow users to do heterogeneous programming in python, such as numbapro, theano, and pycuda. They are usually easier to use, like theano, wrapped entirely the CUDA derivatives in python functions, leaving nothing for users to control the device explicitly. Some other APIs, such as Numbapro, which introduces not only wrapped functions, but also a CUDA dialect. In general, the third party APIs are easier to implement and easier to maintain but more difficult to control and their behaviour is more difficult to understand. + +Good software should find a balance between the software sustainability and the performance. However in this project we are more interested in investigating the behaviour of GPU and exploiting the all the resource of it. CUDA can give us more freedom in tuning. The prototype has been made in C++/CUDA, integrated to python emulator with C-Python API. With the CUDA implementation, we can set a performance reference of future GPU development. + +## Memory issues +In the emulator, we have to at least hold two Npredict x Ntrain (~ 1e8 x 250) matrix (93G) at the same time. GPU memory (~5G) cannot cope with such a large demand. Fortunately, elements inside each parameter vector are independent to each other; in other word, we can truncate the large input vectors (Npredict) into multiple pieces and feed them slice by slice to the _predict_ function. This feature of the emulator solves the memory scarcity issue and moreover enables multiple GPUs implementation in the future. We will go through this in the later sessions. The size of input vector we feed every time to the GPU predict function is 2e5(~400M). + +## Performance of GPU Kernels +| Rank | Operation | Call | Percentage of total wall time (%) | +|:----: |------------ |:----: |:---------------------------------: | +| 1 | MvSub | 10 | 21.1 | +| 2 | Matmul-Ele | 10 | 17.1 | +| 3 | Dist | 1 | 16.3 | +| 4 | Mv | 10 | 6.5 | +| 5 | Matmul | 2 | 3.8 | +| 6 | Exp | 1 | 0.1 | +| | Total | | 99.81 | + +To avoid GPU/CPU memory transferring overhead, the whole _predict_ has been migrated on GPU, rather than only computational intensive operations. The GPU _predict_ is formed by multiple kernels, which correspond to its CPU counterparts, i.e. six major operations. + +On the one hand we made array element wise operations kernels on GPU. GPU is very good at dealing with these problems due to its embarrassingly parallel feature. Meanwhile, cuBLAS has been implemented for the matrix operation component of the _predict_. + +There are three major indicators to measure the GPU optimisation. First is occupancy (active warp/ maximum warp), which indicates if there are sufficient warp (32 threads) of GPU has been launched so that the latency may be covered; the rest two are the indicators of utilisation GPU resource including memory bandwidth (achieved bandwidth/device maximum bandwidth) and computation (flop achieved/maximum flops). + +Array element wise operations still are the most time consuming part of the _predict_. Let us examine the two most expansive kernels (MvSub & Matmul-Ele) here. These kernels are memory bound, which means the memory operations are heavier than computation, i.e. each floating-point operation needs two times memory access. The memory bandwidth of both kernels has been largely saturated, i.e. utilisation of memory bandwidth is 85% respectively, while the computation is relatively low, 20% - 35%, due to its memory bound nature. It means the computation units of GPU are not fully exploited, due to the latency. The occupancy achieved by both kernels are very high, around 75%, which means the GPU has almost tried its best to overlap the latency by launching more concurrent threads. At this point, the further optimisation can only lie upon the improvement to memory access pattern so that the access of memory can be more efficient, which requires much more effort. + +## Predict function performance +Based on the performance benchmark on Emerald K20 GPU. GPU doesn’t introduce a very large overhead in the _predict_. Memory copy between device and host takes only 0.1% of total run time. The entire GPU wall time is only 65.8% of total wall time, which means these once very expansive kernels are no longer the major performance bottleneck anymore, instead the python part of the _predict_ (34.2%), (mainly truncating and distributing data to feed GPU) becomes the new bottleneck. To optimise this part will be the first priority of future optimisation. + +Conventionally, speedup is used as a metric to measure the performance gain of both serial and parallel optimisation. But strictly taken, one can only use it to compare the performance on the same platform. In heterogeneous programming, the speedup is just for giving an idea about the performance rather than a strict measurement, since the chose of different platform as a base of comparison can change entirely the final result. The speedup in our project will be based on the performance of Intel Xeon 5500 CPU on Legion, which is the fastest CPU platform we get in this project. + +The tests have been carried out on 3 nvidia GPU cards, M2070 on Legion (tests are not ready due to the temporary outage of legion GPU), M2090 and K20 on Emerald. K20 is a Kepler architecture GPU, which is the latest architecture of NVidia GPU, while M2070/90 is Fermi. The most distinctive difference between them is that K20 has more than 4 times thread processors than M series. _predict_ function on K20 shows 3 times faster than of M series. It takes K20 875 second to finish the _predict_ function with the vector size of 1e8. + +Since GPU is optimised on its single precision operations, the single precision _predict_ function is much faster than the double precision one. Both implementations have been applied in the source code which can be switched by the users in light of the the performance and the errors. + +![Kernel Performance with Vector Size 2e6](speedup.png) + + +## Software development details + +The main objective of this project is to accelerate the GPU Gaussian Process _predict_ function, which is a very small part of the code. There is no major change of the structure. The predict function C++/CUDA extension compilation is managed by a cross platform compilation tool cmake which has also been integrated in the python setup. + +Unit tests have been introduced to the program for the purpose of verifying the correctness of GPU implementation especially by comparing the GPU results with original python results. Tests are carried out recursively with the changing the lengths input vectors under the CUnit testing framework. + +The code is compiled and tested every overnightly on RSDT continuous integration system Jenkins. + +# Summary and the future work + +In this project, we have investigated the performance _predict_ in the GP emulator. The most computational intensive operations, such as matrix multiplication and distance calculation are not necessary to be the most time consuming operation on CPU, due to its optimisation by scipy and numpy library. The major bottleneck of _predict_ function is the array element wise operations which can be easily implemented on GPU. A GPU implementation prototype has also been made, by porting element wise operations as kernels on GPU and using cuBLAS library to deal with matrix operations. It achieves 65x speedup when running on K20 GPU with single precision. + +The further single GPU optimisation of _predict_ should focus on data distribution part in python. Kernels on GPU still have some space to improve, if the memory access pattern can be improved. It will also be interesting to checkout the performance of some third party APIs such as Numbapro. Third party APIs can reduce the code complexity. In the future, developers can balance code sustainability and performance with the performance reference of the CUDA implementation. + +With the base of single GPU, it seems to be very promising to have multiple GPUs implementation. First, the single GPU implementation has already achieved quite a good performance. Second, the predict function is embarrassingly paralleled which means all inputs can be divided and distributed arbitrarily, as what we did in the single GPU case, so to scale to multiple GPUs should not be difficult. Third, the communications among GPU nodes will be relatively small, which occurs only at the beginning and in the end. Multiple GPUs implementation can fundamentally change the performance to meet the final scientific target. diff --git a/doc/report.pdf b/doc/report.pdf new file mode 100644 index 0000000..770ad02 Binary files /dev/null and b/doc/report.pdf differ diff --git a/doc/speedup.png b/doc/speedup.png new file mode 100644 index 0000000..1e7ae76 Binary files /dev/null and b/doc/speedup.png differ diff --git a/gp_emulator/GaussianProcess.py b/gp_emulator/GaussianProcess.py index f70c38d..47129a2 100644 --- a/gp_emulator/GaussianProcess.py +++ b/gp_emulator/GaussianProcess.py @@ -3,6 +3,8 @@ import numpy as np import scipy.spatial.distance as dist import random +import sys +import _gpu_predict def k_fold_cross_validation(X, K, randomise = False): """ @@ -206,7 +208,7 @@ def learn_hyperparameters ( self, n_tries=15, verbose=False ): self._set_params ( params[idx]) return (log_like[idx], params[idx] ) - def predict ( self, testing, do_unc=True ): + def cpu_predict ( self, testing, do_unc=True ): """Make a prediction for a set of input vectors, as well as calculate the partial derivatives of the emulated model, and optionally, the "emulation uncertainty". @@ -221,34 +223,124 @@ def predict ( self, testing, do_unc=True ): do_unc: flag, optional Calculate the uncertainty (if you don't set this flag, it can shave a few us""" - + ( nn, D ) = testing.shape assert D == self.D - expX = np.exp ( self.theta ) a = dist.cdist ( np.sqrt(expX[:(self.D)])*self.inputs, \ np.sqrt(expX[:(self.D)])*testing, 'sqeuclidean') - a = expX[self.D]*np.exp(-0.5*a) b = expX[self.D] mu = np.dot( a.T, self.invQt) + if do_unc: - var = b - np.sum ( a * np.dot(self.invQ,a), axis=0) + var = b - np.sum ( a * np.dot(self.invQ,a), axis=0) # Derivative and partial derivatives of the function deriv = np.zeros ( ( nn, self.D ) ) for d in xrange ( self.D ): aa = self.inputs[:,d].flatten()[None,:] - testing[:,d].flatten()[:,None] c = a*aa.T - deriv[:, d] = expX[d]*np.dot(c.T, self.invQt) if do_unc: return mu, var, deriv else: return mu, deriv + + def get_gpu_block( self, size, block_size ): + ''' + Distribute a size long vector block_size, and return the start index + and end index of each block. To ensure size of the last block + is not too small, the last two blocks will have a equal size. + ''' + ind_start = range( np.int(0), np.int(size), np.int(block_size) ) + ind_end = np.append( ind_start[1:], size ) + nblocks = len(ind_start) + + # compute the size of last two blocks. + if nblocks > 1: + last_two_block_size = ( ind_end[ nblocks - 1 ] - ind_start[ nblocks - 2 ] ) / 2 + ind_end[ nblocks - 2 ] = ind_start[ nblocks - 2 ] + last_two_block_size + ind_start[ nblocks - 1 ] = ind_end[ nblocks - 2 ] + + assert np.all(ind_end - ind_start <= block_size ) + return ind_start, ind_end + + + def gpu_predict ( self, testing, precision, threshold): + ''' + Parameters: + -------------- + testing: 2D array n_predict * n_inputs + precision: np.float32 / np.float64 + threshold: see predict() threshold. + ''' + import _gpu_predict + n_predict, n_inputs = testing.shape + n_train = self.inputs.shape[0] + theta_size=self.theta.size + + assert n_inputs == self.D + + #_predict_wrap() has to be fed by one dimentional array + inputs = precision(self.inputs.reshape(self.inputs.shape[0] * self.inputs.shape[1])) + invQt = precision(self.invQt) + invQ = precision(self.invQ.reshape(self.invQ.shape[0] * self.invQ.shape[1])) + expX = precision(np.exp(self.theta)) + + result = [] + error = [] + deriv = np.array([]).reshape((0, n_inputs)) + ind_start, ind_end = self.get_gpu_block( n_predict, threshold ) + + for block_start, block_end in zip(ind_start, ind_end): + testing_block = testing[block_start:block_end,:] + testing_block = testing_block.reshape( testing_block.shape[0] * testing_block.shape[1] ) + + n_predict_block = np.int(block_end - block_start) + result_block = np.zeros( n_predict_block ) + error_block = np.zeros( n_predict_block ) + deriv_block = np.zeros( n_predict_block * n_inputs ) + + testing_block = precision(testing_block) + result_block = precision(result_block) + error_block = precision(error_block) + deriv_block = precision(deriv_block) + + _gpu_predict.predict_wrap( + expX, inputs, invQt, invQ, testing_block, + result_block, error_block, deriv_block, + n_predict_block, n_train, n_inputs, theta_size) + + result = np.append(result, result_block) + error = np.append(error, error_block) + #deriv produced by gpu is transposed, so here we need transpose them back. + deriv = np.append(deriv, deriv_block.reshape( n_inputs, n_predict_block ).T, axis = 0) + + return result, error, deriv + + + + def predict(self, testing, do_unc = True, is_gpu = False, precision = np.float64, threshold = 2e5): + ''' + Parameters: + -------------- + testing: 2D array n_predict * n_inputs + precision: np.float32 / np.float64 + do_unc: tag to switch on or off calculation of the uncertainty. It can only affect cpu_predict() + threshold: is the maximum number of n_predict that gpu_predict() can deal with. + If the n_predict is larger than the threshold, data will be truncated, + and gpu_predict will be excuted for multiple time. + ''' + if is_gpu == True: + return self.gpu_predict(testing, precision, threshold = threshold) + else: + return self.cpu_predict(testing, do_unc) + + def hessian ( self, testing ): '''calculates the hessian of the GP for the testing sample. @@ -272,52 +364,3 @@ def hessian ( self, testing ): cc = a*(aa.T) hess[:, d,d2] = np.dot(cc.T, self.invQt) return hess - -if __name__ == "__main__": - import matplotlib.pyplot as plt - np.set_printoptions(precision=2, suppress=True) - #input_obs = np.array ( [[-4, -3, -1, 0, 2]]).T - #target1 = np.array ([-2, 0, 1., 2., -1]) - #data = np.loadtxt ("mvreg.dat", delimiter="," ) - #target1 = data[ :,0] - #target2 = data[ :,1] - #target3 = data[ :,2] - #input_obs = data[ :, 3: ] - wheat = np.loadtxt("argentine_wheat.dat")[:, 1:] - yields = wheat[:,0] - mu = yields.mean() - sigma = yields.std() - yields = (yields - mu ) / sigma - wheat [:, 0] = yields - rmse = [] - for ( train,validate) in k_fold_cross_validation ( wheat, 5, randomise=True): - train = np.array ( train ) - validate = np.array ( validate ) - yields_t = train [ :, 0] - inputs_t = train [ :, 1:] - yields_v = validate [ :, 0] - inputs_v = validate [ :, 1:] - gp = GaussianProcess ( inputs_t, yields_t ) - theta_min= gp.learn_hyperparameters (n_tries=2) - pred_mu, pred_var, par_dev = gp.predict ( inputs_v ) - print "TEST" - print inputs_v - print pred_mu, pred_var, par_dev - r = ( yields_v - pred_mu )**2#/pred_var - rmse.append ( [ np.sqrt(r.mean()), theta_min[1] ]) - - - #gp.theta[2] = 0. - ######theta = gp.theta - ######print theta - ######gp.predict ( input_obs ) - - ######x = np.linspace(-5,5,100) - - ######(mu, var ) = gp.predict ( x[:, np.newaxis]) - ######plt.plot ( x, mu, '-r' ) - ######plt.fill_between ( x, mu+np.sqrt(var), mu-np.sqrt(var), color='0.8') - #######plt.errorbar ( x, mu, yerr=np.sqrt(var)*0.5 ) - ######plt.plot ( input_obs, target1, 'gs' ) - ######plt.title("theta: %s" % theta ) - ######plt.show() diff --git a/gp_emulator/__init__.py b/gp_emulator/__init__.py index 9d486dc..0411b6a 100644 --- a/gp_emulator/__init__.py +++ b/gp_emulator/__init__.py @@ -1,4 +1,3 @@ - from GaussianProcess import GaussianProcess, k_fold_cross_validation from multivariate_gp import MultivariateEmulator from lhd import lhd diff --git a/gp_emulator/gpu/CMakeLists.txt b/gp_emulator/gpu/CMakeLists.txt new file mode 100644 index 0000000..b4246d7 --- /dev/null +++ b/gp_emulator/gpu/CMakeLists.txt @@ -0,0 +1,17 @@ +cmake_minimum_required (VERSION 2.8) + +cuda_add_library(_gpu_predict _gpu_predict.cpp + predict.cu + computeTranspose.cu + kernel_vectorTimesMatrix.cu + kernel_init_array.cu + kernel_cdist.cu + kernel_matrixExp.cu + kernel_elementwiseMult.cu + kernel_scalarMinusVec.cu + kernel_rowSum.cu + kernel_crossMinus.cu MODULE) +set_target_properties(_gpu_predict PROPERTIES PREFIX "") +cuda_add_cublas_to_target(_gpu_predict) +target_link_libraries(_gpu_predict ${PYTHON_LIBRARIES}) +subdirs(tests) diff --git a/gp_emulator/gpu/_gpu_predict.cpp b/gp_emulator/gpu/_gpu_predict.cpp new file mode 100644 index 0000000..6f5cd39 --- /dev/null +++ b/gp_emulator/gpu/_gpu_predict.cpp @@ -0,0 +1,181 @@ +/**********************************//** + * source code of share library _gpu_predict.so + * importing inputs from python + * calling GPU predict function by predict_wrap() + * python utility functions are based on the following template: + * http://wiki.scipy.org/Cookbook/C_Extensions/NumPy_arrays + **********************************/ + +#include "gpu_predict.h" +#include + +/*********************************//** + * Set up the methods table + *********************************/ +static PyMethodDef gpuMethods[] = +{ + {"predict_wrap", predict_wrap, METH_VARARGS}, + {NULL, NULL} /* Sentinel - marks the end of this structure */ +}; + +/**********************************//** + * Initialize the C_test functions + * Module name must be _C_arraytest in compile and linked + **********************************/ +extern "C"{ +void init_gpu_predict() +{ + (void) Py_InitModule("_gpu_predict", gpuMethods); + import_array(); // Must be present for NumPy. Called first after above line. +} +} + + + +/**********************************//** + * Check the dimension and data type + * of input python arrays. + **********************************/ +void checkRealVector(PyArrayObject *vec) +{ + if( sizeof( real ) == sizeof( double ) ) + { + if (vec->descr->type_num != NPY_DOUBLE || vec->nd != 1) + { + printf( "In checkRealVector: array must be of type Double vector\n"); + exit(EXIT_FAILURE); + } + + } + if( sizeof( real ) == sizeof( float ) ) + { + if (vec->descr->type_num != NPY_FLOAT32 || vec->nd != 1) + { + printf("In checkRealVector: array must be of type Float32 vector\n"); + exit(EXIT_FAILURE); + } + } +} + + + +/**********************************//** + * Convert data (vector) from PyArrayObject to real + **********************************/ +real *pyvector_to_Carrayptrs(PyArrayObject *arrayin) +{ + checkRealVector(arrayin); + return (real *) arrayin->data; /* pointer to arrayin data as real */ +} + + + +/**********************************//** + * getPredictDataFromPython + * Assign multiple input python arrays, which are all in 1D + * to C array. + **********************************/ +void getPredictDataFromPython(PyObject *args, real **c_theta_exp, real **c_invQt, real **c_invQ, + real **c_predict, real **c_train, + real **c_result, real **c_error, real **c_deriv, + int *Npredict, int *Ntrain, int *Ninputs, int *theta_size) +{ + PyArrayObject *py_theta_exp, *py_train, *py_invQt, *py_invQ, *py_predict; + PyArrayObject *py_result, *py_error, *py_deriv; + + PyArg_ParseTuple ( args, "O!O!O!O!O!O!O!O!iiii", + &PyArray_Type, &py_theta_exp, + &PyArray_Type, &py_train, + &PyArray_Type, &py_invQt, + &PyArray_Type, &py_invQ, + &PyArray_Type, &py_predict, + &PyArray_Type, &py_result, + &PyArray_Type, &py_error, + &PyArray_Type, &py_deriv, + Npredict, Ntrain, Ninputs, theta_size); + + *c_theta_exp = pyvector_to_Carrayptrs( py_theta_exp); + *c_train = pyvector_to_Carrayptrs( py_train ); + *c_invQt = pyvector_to_Carrayptrs( py_invQt ); + *c_invQ = pyvector_to_Carrayptrs( py_invQ ); + *c_predict = pyvector_to_Carrayptrs( py_predict ); + *c_result = pyvector_to_Carrayptrs( py_result ); + *c_error = pyvector_to_Carrayptrs( py_error ); + *c_deriv = pyvector_to_Carrayptrs( py_deriv ); +} + + + +/**********************************//** + * predict_wrap + * 1) call getPredictDataFromPython; + * 2) transpose 2D arrays (data arranged in 1D) to column major; + * 3) calling GPU predict function; + **********************************/ +PyObject *predict_wrap ( PyObject *self, PyObject *args ) +{ + int i; + int Npredict, Ntrain, Ninputs, theta_size; + real *c_theta_exp; + real *c_train, *c_invQt, *c_invQ, *c_predict; + real *c_result, *c_error, *c_deriv; + + getPredictDataFromPython(args, &c_theta_exp, &c_invQt, &c_invQ, + &c_predict, &c_train, + &c_result, &c_error, &c_deriv, + &Npredict, &Ntrain, &Ninputs, &theta_size); + + //transpose 2D array to column major to cope with cublas + real *c_invQ_T, *c_train_T, *c_predict_T; + c_invQ_T = computeTranspose( c_invQ, Ntrain, Ntrain ); + c_train_T = computeTranspose( c_train, Ninputs, Ntrain ); + c_predict_T = computeTranspose( c_predict, Ninputs, Npredict); + + //get c_theta_exp_sqrt + real *c_theta_exp_sqrt; + c_theta_exp_sqrt = (real *)malloc( sizeof(real) * theta_size ); + for( i = 0; i < theta_size; i++ ) + { + c_theta_exp_sqrt[i] = sqrt( c_theta_exp[i] ); + } + + //predict + gpuPredict gpu_predict( + c_theta_exp, c_theta_exp_sqrt, + c_invQt, c_invQ_T, + c_predict_T,c_train_T, + c_result, c_error, c_deriv, + Npredict, Ntrain, Ninputs, theta_size ); + gpu_predict.predict(); + + + free(c_invQ_T); + free(c_train_T); + free(c_predict_T); + free(c_theta_exp_sqrt); + + Py_INCREF(Py_None); + return Py_None; +} + + + + + + + + + + + + + + + + + + + + + + diff --git a/gp_emulator/gpu/computeTranspose.cu b/gp_emulator/gpu/computeTranspose.cu new file mode 100644 index 0000000..1041385 --- /dev/null +++ b/gp_emulator/gpu/computeTranspose.cu @@ -0,0 +1,18 @@ +#include "gpu_predict.h" + +real *computeTranspose(const real *matrix, const int lead_dim_in, const int lead_dim_out) +{ + int x, y; + real * temp; + temp = ( real *)malloc(sizeof(real) * lead_dim_in * lead_dim_out); + + for ( y = 0; y < lead_dim_out; ++y ) + { + for ( x = 0; x < lead_dim_in; ++x ) + { + temp[(x * lead_dim_out) + y] = matrix[(y * lead_dim_in) + x]; + + } + } + return(temp); +} diff --git a/gp_emulator/gpu/gpu_predict.h b/gp_emulator/gpu/gpu_predict.h new file mode 100644 index 0000000..e31aabd --- /dev/null +++ b/gp_emulator/gpu/gpu_predict.h @@ -0,0 +1,155 @@ +/* A file to test importing C modules for handling arrays to Python */ + +#include "Python.h" +#include "arrayobject.h" +#include +#include +#include +#include +#include +/* Includes, cuda */ +#include +#include +#include + +#define CUDA_BLOCK 2 +#define MIN_NPREDICT 1000 //ensure there will be enough amount of threads to launch. +#define MAX_NUM_THREAD 1024 //sm_20 +#define MAX_NUM_BLOCK 65536 + +#ifdef DOUBLE__PRECISION + #define real double +#else + #define real float +#endif + +#ifdef DOUBLE__PRECISION + #define CUBLAS_GEMV cublasDgemv + #define CUBLAS_GEMM cublasDgemm + #define CUBLAS_GEAM cublasDgeam +#else + #define CUBLAS_GEMV cublasSgemv + #define CUBLAS_GEMM cublasSgemm + #define CUBLAS_GEAM cublasSgeam +#endif + +__forceinline__ __host__ __device__ +int IDX2D(int row, int col, int lead_dim) +{ + return(((col)*(lead_dim))+(row)); +} + +PyArrayObject *pyvector(PyObject *objin); +real*pyvector_to_Carrayptrs(PyArrayObject *arrayin); +real **pymatrix_to_Carrayptrs(PyArrayObject *arrayin); +real **ptrvector(long n); +PyObject *predict_wrap ( PyObject *self, PyObject *args ); + + +void getPredictDataFromPython(PyObject *args, real **c_theta_exp, real **c_invQt, real **c_invQ, + real **c_testing, real **c_inputs, + real **c_mu, real **c_var, real **c_deriv, + int *N, int *M, int *D, int *theta_size); + +#ifdef __cplusplus +extern "C"{ +#endif +real *computeTranspose(const real *matrix, const int size_in, const int size_out); +void gpu_vectorTimesMatrix(const real *A, const real *v, real *res, int nrows, int ncols); +void gpu_init_array(real *vec, const int init_val, const int vec_len); +void gpu_cdist(const real *input1, const real *input2, real *output, const int nrow1, const int ncol1, + const int nrow2, const int ncol2); +void gpu_matrixExp( real *matrix,const real alpha,const real beta, const int size ); +void gpu_elementwiseMult( const real *v1, real *v2, const int size ); +void gpu_scalarMinusVec( real *matrix, const real scalar, const int size ); +real* gpu_rowSum(const real *A, const int A_nrows,const int A_ncols); +void gpu_crossMinus(const real *v1, const real *v2, real *mat_res, const int v1_len, const int v2_len); +#ifdef __cplusplus +} +#endif + +class gpuPredict +{ + real *c_result, *c_error, *c_deriv; + const real *c_theta_exp, *c_theta_exp_sqrt; + const real *c_invQt, *c_invQ, *c_predict, *c_train, *c_inputs; + const int Npredict, Ntrain, Ninputs; + const int theta_size; + + real *d_result, *d_error, *d_deriv; + real *d_theta_exp, *d_theta_exp_sqrt; + real *d_invQt, *d_invQ, *d_predict, *d_train, *d_inputs; + real *d_dist_matrix; real *d_dist_matrix_T; + cublasHandle_t handle; + + public: + gpuPredict(real *ctheta_exp, + real *ctheta_exp_sqrt, + real *cinvQt, + real *cinvQ, + real *cpredict, + real *ctrain, + real *cresult, + real *cerror, + real *cderiv, + int npredict, + int ntrain, + int ninputs, + int thetasize): + c_theta_exp(ctheta_exp), + c_theta_exp_sqrt(ctheta_exp_sqrt), + c_invQt(cinvQt), + c_invQ(cinvQ), + c_predict(cpredict), + c_train(ctrain), + c_result(cresult), + c_error(cerror), + c_deriv(cderiv), + Npredict(npredict), + Ntrain(ntrain), + Ninputs(ninputs), + theta_size(thetasize){}; + + + real * gpu_transpose(real *, const int, const int); + void init_gpu(void); + void compute_distance(void); + void compute_result(void); + void compute_error(void); + void compute_deriv(void); + void predict(void); + void free_gpu(void); +}; + + + + + + + +// error check macros +#define cudaCheckErrors(msg) \ + do { \ + cudaError_t __err = cudaGetLastError(); \ + if (__err != cudaSuccess) { \ + fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \ + msg, cudaGetErrorString(__err), \ + __FILE__, __LINE__); \ + fprintf(stderr, "*** FAILED - ABORTING\n"); \ + exit(1); \ + } \ + } while (0) + +// for CUBLAS V2 API +#define cublasCheckErrors(fn) \ + do { \ + cublasStatus_t __err = fn; \ + if (__err != CUBLAS_STATUS_SUCCESS) { \ + fprintf(stderr, "Fatal cublas error: %d (at %s:%d)\n", \ + (int)(__err), \ + __FILE__, __LINE__); \ + fprintf(stderr, "*** FAILED - ABORTING\n"); \ + exit(1);\ + } \ + } while (0) + diff --git a/gp_emulator/gpu/kernel_cdist.cu b/gp_emulator/gpu/kernel_cdist.cu new file mode 100644 index 0000000..f98a76b --- /dev/null +++ b/gp_emulator/gpu/kernel_cdist.cu @@ -0,0 +1,47 @@ +/*********************************************//** + * Squared Euclidiean distance function: + * - Equivalent to scipy cdist() function + * - In1_ld and In2_ld, are leading dimention of the input1 and input2, + * in our case it should be always the column. + * + *********************************************/ +#include "gpu_predict.h" +#define min(a, b) (((a) < (b)) ? (a) : (b)) +#define CDIST_NTHREAD_X 256 +#define CDIST_NTHREAD_Y 4 +#define CDIST_NTHREAD_Z 1 + +__global__ +void kernel_cdist(const real *input1, const real *input2, real *output, const int nrow1, const int nrow2, const int ncol) +{ + int ix, iy, iz; + ix = blockIdx.x * blockDim.x + threadIdx.x;//N + iy = blockIdx.y * blockDim.y + threadIdx.y;//M + iz = blockIdx.z * blockDim.z + threadIdx.z; + if( ix < nrow2 && iy < nrow1 && iz < ncol) + output[IDX2D(ix, iy, nrow2)] += pow(input1[IDX2D(iy, iz, nrow1)] - input2[IDX2D(ix, iz, nrow2)],2); +} + + +void gpu_cdist(const real *input1, const real *input2, real *output, const int nrow1, const int ncol1, const int nrow2, const int ncol2) +{ + if( nrow2 < MIN_NPREDICT ) + { + printf("gpu_cdist: %d [nrow2(Npredict)] < %d\n",nrow2, MIN_NPREDICT); + exit(EXIT_FAILURE); + } + + dim3 nthread, nblock; + nthread.x = CDIST_NTHREAD_X; + nthread.y = CDIST_NTHREAD_Y; + nthread.z = CDIST_NTHREAD_Z; + + nblock.x = ceil( float(nrow2) / float(nthread.x) ); + nblock.y = ceil( float(nrow1) / float(nthread.y) ); + nblock.z = ceil( float(ncol1) / float(nthread.z) ); + + kernel_cdist<<>>(input1, input2, output, nrow1, nrow2, ncol1); +} + + + diff --git a/gp_emulator/gpu/kernel_crossMinus.cu b/gp_emulator/gpu/kernel_crossMinus.cu new file mode 100644 index 0000000..80c4cac --- /dev/null +++ b/gp_emulator/gpu/kernel_crossMinus.cu @@ -0,0 +1,35 @@ +#include "gpu_predict.h" +/*********************************************//** + * cross minus: + * aa_{ix, iy} = inputs_{ix} - testing_{iy} + *********************************************/ +__global__ +void kernel_crossMinus(const real *vec1,const real *vec2, real *matrix_result, const int vec1_len, const int vec2_len) +{ + int ix = blockIdx.x * blockDim.x + threadIdx.x; + int iy = blockIdx.y * blockDim.y + threadIdx.y; + + if(ix < vec1_len && iy < vec2_len) + matrix_result[IDX2D(ix, iy, vec1_len)] = vec1[ix] - vec2[iy]; +} + + + +void gpu_crossMinus(const real *vec1, const real *vec2, real *matrix_result, const int vec1_len, const int vec2_len) +{ + dim3 nthread, nblock; + if( vec2_len < MAX_NUM_THREAD ) + nthread.y = vec2_len; + else + nthread.y = MAX_NUM_THREAD; + + nthread.x = 1; + nblock.x = vec1_len; + nblock.y = ceil( float(vec2_len) / float(nthread.y) ); + + kernel_crossMinus<<< nblock, nthread >>>(vec1, vec2, matrix_result, vec1_len, vec2_len); +} + + + + diff --git a/gp_emulator/gpu/kernel_elementwiseMult.cu b/gp_emulator/gpu/kernel_elementwiseMult.cu new file mode 100644 index 0000000..1c397ac --- /dev/null +++ b/gp_emulator/gpu/kernel_elementwiseMult.cu @@ -0,0 +1,42 @@ +/*********************************************//** + * vector elementwise multiplication + * vector2_{i} = vector2_{i} * vector1_{i} + *********************************************/ +#include "gpu_predict.h" +__global__ +void kernel_elementwiseMult(const real *vector1, real *vector2, const int size) +{ + int i, index; + int ix = blockIdx.x * blockDim.x + threadIdx.x; + for( i = 0; i < CUDA_BLOCK; ++i ) + { + index = ix * CUDA_BLOCK + i; + if( index < size) + { + vector2[index] = vector2[index] * vector1[index]; + } + } +} + + +void gpu_elementwiseMult( const real *vector1, real *vector2, const int size ) +{ + int nthread, nblock; + + if( size < float(MAX_NUM_THREAD) / float(CUDA_BLOCK) ) + { + printf("gpu_elementwiseMult: size = %d [ < MAX_NUM_THREAD / CUDA_BLOCK ].\n", size); + exit(EXIT_FAILURE); + } + + nthread = MAX_NUM_THREAD; + nblock = ceil( float(size) / float(CUDA_BLOCK) / float(nthread) ); + + if( nblock > MAX_NUM_BLOCK ) + { + printf("gpu_elementwiseMult: nblock outside the range of [1, MAX_NUM_BLOCK]\n"); + exit(EXIT_FAILURE); + } + + kernel_elementwiseMult<<>>(vector1, vector2, size); +} diff --git a/gp_emulator/gpu/kernel_init_array.cu b/gp_emulator/gpu/kernel_init_array.cu new file mode 100644 index 0000000..39f6ce0 --- /dev/null +++ b/gp_emulator/gpu/kernel_init_array.cu @@ -0,0 +1,45 @@ +/*********************************************//** + * initialise an array [real] with an assigned + * value init_val. + *********************************************/ +#include "gpu_predict.h" +__global__ +void kernel_init_array(real *vec, const real init_val, const int vec_len) +{ + int i, index; + int ix = blockIdx.x * blockDim.x + threadIdx.x; + for( i = 0; i < CUDA_BLOCK; ++i ) + { + index = ix * CUDA_BLOCK + i; + if( index < vec_len) + { + vec[index] = init_val; + } + } +} + +void gpu_init_array(real *vec, const int init_val, const int vec_len) +{ + int nthread, nblock; + + if( CUDA_BLOCK * vec_len < MAX_NUM_THREAD) + { + nthread = ceil( vec_len / CUDA_BLOCK ); + nblock = 1; + } + else + { + if( vec_len / CUDA_BLOCK > MAX_NUM_BLOCK * MAX_NUM_THREAD ) // largest block size for sm_2.0 + { + printf("gpu_init_array: vector length / CUDA_BLOCK = %d > MAX_NUM_BLOCK * MAX_NUM_THREAD.\n", vec_len/CUDA_BLOCK); + exit(EXIT_FAILURE); + } + nthread = MAX_NUM_THREAD; + nblock = ceil( float(vec_len) / nthread / float(CUDA_BLOCK) ); + } + kernel_init_array<<< nblock, nthread >>>(vec, init_val, vec_len); +} + + + + diff --git a/gp_emulator/gpu/kernel_matrixExp.cu b/gp_emulator/gpu/kernel_matrixExp.cu new file mode 100644 index 0000000..95b56de --- /dev/null +++ b/gp_emulator/gpu/kernel_matrixExp.cu @@ -0,0 +1,38 @@ +/*********************************************//** + * Do the following operation: + * - matrix_{i} = beta * e^{alpha * matrix_{i}} + *********************************************/ +#include "gpu_predict.h" +__global__ +void kernel_matrixExp(real *matrix, const real alpha, const real beta, const int size) +{ + int i, index; + int ix = blockIdx.x * blockDim.x + threadIdx.x; + for( i = 0; i < CUDA_BLOCK; ++i ) + { + index = ix * CUDA_BLOCK + i; + if( index < size) + matrix[ index ] = beta * exp( alpha * matrix[ index ]); + } +} + +void gpu_matrixExp( real *matrix,const real alpha,const real beta, const int size ) +{ + int nthread, nblock; + + if( size < float( MAX_NUM_THREAD ) / float(CUDA_BLOCK) ) + { + printf("gpu_matrixExp: size = %d [ < MAX_NUM_THREAD / CUDA_BLOCK ].\n", size); + exit(EXIT_FAILURE); + } + + if( size > MAX_NUM_BLOCK * MAX_NUM_THREAD ) + { + printf("gpu_matrixExp: size = %d [ > MAX_NUM_BLOCK * MAX_NUM_THREAD / CUDA_BLOCK ].\n", size); + exit(EXIT_FAILURE); + } + + nthread = MAX_NUM_THREAD; + nblock = ceil( float(size) / float(CUDA_BLOCK) / float(nthread) ); + kernel_matrixExp<<>>(matrix, alpha, beta, size); +} diff --git a/gp_emulator/gpu/kernel_rowSum.cu b/gp_emulator/gpu/kernel_rowSum.cu new file mode 100644 index 0000000..49e5cf7 --- /dev/null +++ b/gp_emulator/gpu/kernel_rowSum.cu @@ -0,0 +1,28 @@ +/*********************************************//** + * row sum: + * return vector rowSum(matrix) + *********************************************/ +#include "gpu_predict.h" + +real* gpu_rowSum(const real *matrix, const int matrix_nrows,const int matrix_ncols) +{ + cublasHandle_t handle; + cublasCreate(&handle); + + real alpha = 1.f; + real beta = 0.f; + real *vec_one; + real *d_result; + + cudaMalloc((void **)&vec_one, sizeof(real) * matrix_ncols); + cudaMalloc((void **)&d_result, sizeof(real) * matrix_ncols); + + gpu_init_array(vec_one, 1.0, matrix_ncols); + gpu_init_array(d_result, 0.0, matrix_ncols); + + cublasCheckErrors(CUBLAS_GEMV(handle, CUBLAS_OP_T, matrix_nrows, matrix_ncols, &alpha, matrix, matrix_nrows, vec_one, 1, &beta, d_result, 1)); + + cudaFree(vec_one); + cublasDestroy(handle); + return d_result; +} diff --git a/gp_emulator/gpu/kernel_scalarMinusVec.cu b/gp_emulator/gpu/kernel_scalarMinusVec.cu new file mode 100644 index 0000000..c24b959 --- /dev/null +++ b/gp_emulator/gpu/kernel_scalarMinusVec.cu @@ -0,0 +1,38 @@ +/*********************************************//** + * vector scalar operation, update vec by, + * vec = scalar - vec + *********************************************/ +#include "gpu_predict.h" +__global__ +void kernel_scalarMinusVec(real *vec, const real scalar, const int size) +{ + int ix = blockIdx.x * blockDim.x + threadIdx.x; + if( ix < size ) + vec[ix] = scalar - vec[ix]; +} + + +void gpu_scalarMinusVec( real *matrix, const real scalar, const int size ) +{ + int nthread, nblock; + if( size > MAX_NUM_BLOCK * MAX_NUM_THREAD ) + { + printf("gpu_scalarMinusVec: size = %d [ > MAX_NUM_BLOCK * MAX_NUM_THREAD ].\n", size); + exit(EXIT_FAILURE); + } + + if( size < MAX_NUM_THREAD ) + { + nthread = size; + nblock = 1; + } + else + { + nthread = MAX_NUM_THREAD; + nblock = ceil( float(size) / float(nthread) ); + } + + kernel_scalarMinusVec<<>>(matrix, scalar, size); +} + + diff --git a/gp_emulator/gpu/kernel_vectorTimesMatrix.cu b/gp_emulator/gpu/kernel_vectorTimesMatrix.cu new file mode 100644 index 0000000..63283be --- /dev/null +++ b/gp_emulator/gpu/kernel_vectorTimesMatrix.cu @@ -0,0 +1,41 @@ +/*********************************************//** + * vector matrix elementwise multiplication + * res_{ix,iy} = matrix_{ix,iy} * v_{iy} + * matrix_lead_dim: leading dimension of matrix matrix + * matrix_second_dim: secondary dimension of matrix matrix + *********************************************/ +#include "gpu_predict.h" +#define VTM_THREADX 100 //best thread number in x dim for vectorTimesMatrix +__global__ +void kernel_vectorTimesMatrix(const real *matrix, const real * vector, real *res, int matrix_lead_dim, int matrix_second_dim) +{ + int ix, iy; + ix = blockIdx.x * blockDim.x + threadIdx.x; + iy = blockIdx.y * blockDim.y + threadIdx.y; + if( ix < matrix_lead_dim && iy < matrix_second_dim) + res[IDX2D(ix, iy, matrix_lead_dim)] = matrix[IDX2D(ix, iy, matrix_lead_dim)] * vector[iy]; +} + + +void gpu_vectorTimesMatrix(const real *matrix, const real *vector, real *res, int nrows, int ncols) +{ + dim3 nthread, nblock; + if( nrows > MIN_NPREDICT ) + { + nthread.x = VTM_THREADX; + nthread.y = ncols; + nblock.x = ceil(float(nrows)/VTM_THREADX); + nblock.y = 1; + } + else + { + nthread.x = 1; + nthread.y = ncols; + nblock.x = nrows; + nblock.y = 1; + } + kernel_vectorTimesMatrix<<>>(matrix, vector, res, nrows, ncols); +} + + + diff --git a/gp_emulator/gpu/predict.cu b/gp_emulator/gpu/predict.cu new file mode 100644 index 0000000..87737c3 --- /dev/null +++ b/gp_emulator/gpu/predict.cu @@ -0,0 +1,182 @@ + +/*********************************************//** + * CUDA implementation of derived from GaussianProcess + * in python code. + * + * Sinan Shi + *********************************************/ +#include "gpu_predict.h" +#include + +void gpuPredict::init_gpu(void) +{ + cublasCreate(&handle); + cudaMalloc( (void **)&d_theta_exp, sizeof(real) * theta_size ); + cudaMalloc( (void **)&d_theta_exp_sqrt, sizeof(real) * theta_size ); + cudaMalloc( (void **)&d_invQt, sizeof(real) * Ntrain); + cublasCheckErrors(cublasSetVector( + theta_size, sizeof(real), c_theta_exp_sqrt, 1, d_theta_exp_sqrt, 1 )); + cublasCheckErrors(cublasSetVector( + theta_size, sizeof(real), c_theta_exp, 1, d_theta_exp, 1 )); + cublasCheckErrors(cublasSetVector( + Ntrain, sizeof(real), c_invQt, 1, d_invQt, 1)); + + //allocate and copy matrix on device + cudaMalloc( (void **)&d_train, sizeof(real) * Ntrain * Ninputs ); + cudaMalloc( (void **)&d_invQ, sizeof(real) * Ntrain * Ntrain ); + cudaMalloc( (void **)&d_predict, sizeof(real) * Npredict * Ninputs ); + cublasCheckErrors(cublasSetMatrix( + Ntrain, Ninputs, sizeof(real), c_train, Ntrain, d_train, Ntrain )); + cublasCheckErrors(cublasSetMatrix( + Ntrain, Ntrain, sizeof(real), c_invQ, Ntrain, d_invQ, Ntrain )); + cublasCheckErrors(cublasSetMatrix( + Npredict, Ninputs, sizeof(real), c_predict, Npredict, d_predict, Npredict)); +} + +/*********************************//* + * Euclidian distance calculation + * 1) res_mv1 = theta_exp_sqrt_{,Ninputs} * train + * 2) res_mv2 = theta_exp_sqrt_{,Ninputs} * predict + * 3) dist_matrix = cidist(res_mv1, res_mv2) + * 4) dist_matrix = -0.5 * exp(expX_{Ninputs}) + * Notice: distance is equivalent to distance^T (a^T) in python + ********************************/ +void gpuPredict::compute_distance(void) +{ + real *d_res_mv1, *d_res_mv2; + cudaMalloc((void **)&d_res_mv1, sizeof(real) * Ntrain * Ninputs); + cudaMalloc((void **)&d_res_mv2, sizeof(real) * Npredict * Ninputs); + cudaMalloc((void **)&d_dist_matrix, sizeof(real) * Ntrain * Npredict); + + gpu_vectorTimesMatrix(d_train, d_theta_exp_sqrt, d_res_mv1, Ntrain, Ninputs); + gpu_vectorTimesMatrix(d_predict, d_theta_exp_sqrt, d_res_mv2, Npredict, Ninputs); + gpu_init_array( d_dist_matrix, 0.0, Npredict * Ntrain ); + gpu_cdist(d_res_mv1, d_res_mv2, d_dist_matrix, Ntrain, Ninputs, Npredict, Ninputs); + gpu_matrixExp(d_dist_matrix, -0.5, c_theta_exp[Ninputs], Ntrain * Npredict); + + cudaFree(d_res_mv1); + cudaFree(d_res_mv2); +} + +/*********************************//* + * compute result: + * c_result = dist_matrix * invQt (dot product) + ********************************/ +void gpuPredict::compute_result(void) +{ + real *d_result; + cudaMalloc((void **)&d_result, sizeof(real) * Npredict); + real alpha = 1.f; + real beta = 0.f; + cublasCheckErrors(CUBLAS_GEMV(handle, CUBLAS_OP_N, Npredict, Ntrain, &alpha, d_dist_matrix, + Npredict, d_invQt, 1, &beta, d_result, 1)); + cudaMemcpy(c_result, d_result, sizeof(real) * Npredict, cudaMemcpyDeviceToHost); + cudaFree(d_result); +} + + +real * gpuPredict::gpu_transpose(real *d_matrix, const int nrow, const int ncol) +{ + real *d_matrix_T; + real alpha = 1.f; + real beta = 0.f; + cudaMalloc((void **)&d_matrix_T, sizeof(real) * nrow * ncol ); + cublasCheckErrors(CUBLAS_GEAM(handle, CUBLAS_OP_T, CUBLAS_OP_N, + nrow, ncol, &alpha, d_matrix, ncol, &beta, + d_matrix, nrow, d_matrix_T, nrow)); + return( d_matrix_T ); +} + +/********************************* + * compute error: + * c_error = b - rowsum(a * dot(invQ, d_dist_matrix_T)) + * arguments d_invQ, d_dist_matrix have been freed. + ********************************/ +void gpuPredict::compute_error() +{ + real alpha = 1.f; + real beta = 0.f; + real *d_res_dot; + real *d_error; + + cudaMalloc((void **)&d_res_dot, sizeof(real) * Ntrain * Npredict); + cublasCheckErrors(CUBLAS_GEMM( + handle, CUBLAS_OP_N, CUBLAS_OP_T, + Ntrain, Npredict, Ntrain, + &alpha, d_invQ, Ntrain, + d_dist_matrix, Npredict, + &beta, d_res_dot, Ntrain)); // dot(invQ, d_dist_matrix_T) + d_dist_matrix_T = gpuPredict::gpu_transpose(d_dist_matrix, Ntrain, Npredict); + gpu_elementwiseMult(d_dist_matrix_T, d_res_dot, Ntrain * Npredict); + d_error = gpu_rowSum(d_res_dot, Ntrain, Npredict); + gpu_scalarMinusVec(d_error, c_theta_exp[Ninputs], Npredict ); + + cudaMemcpy(c_error, d_error, sizeof(real) * Npredict, cudaMemcpyDeviceToHost); + + cudaFree(d_dist_matrix); + cudaFree(d_invQ); + cudaFree(d_error); + cudaFree(d_res_dot); +} + + + +/*********************************//* + * compute deriv: + ********************************/ +void gpuPredict::compute_deriv( void ) +{ + int i; + real alpha; + real beta = 0.f; + + real *d_deriv, *d_aa; + cudaMalloc((void **)&d_deriv, sizeof(real) * Npredict ); + cudaMalloc((void **)&d_aa, sizeof(real) * Ntrain * Npredict); + real *ptr_train, *ptr_predict, *ptr_deriv; + ptr_train = d_train; + ptr_predict = d_predict; + ptr_deriv = c_deriv; + + for( i = 0; i < Ninputs; ++i){ + gpu_crossMinus(ptr_train, ptr_predict, d_aa, Ntrain, Npredict ); + ptr_train = ptr_train + Ntrain; + ptr_predict = ptr_predict + Npredict; + alpha = c_theta_exp[i]; + gpu_elementwiseMult(d_dist_matrix_T, d_aa, Ntrain * Npredict); + cublasCheckErrors(CUBLAS_GEMV(handle, CUBLAS_OP_T, Ntrain, Npredict, &alpha, d_aa, Ntrain, d_invQt, 1, &beta, d_deriv,1)); + + cudaMemcpy(ptr_deriv, d_deriv, sizeof(real) * Npredict, cudaMemcpyDeviceToHost); + ptr_deriv = ptr_deriv + Npredict; + } + + cudaFree(d_deriv); + cudaFree(d_aa); +} + +void gpuPredict::free_gpu(void) +{ + cublasDestroy(handle); + cudaFree(d_invQt); + cudaFree(d_dist_matrix_T); + cudaFree(d_train); + cudaFree(d_theta_exp); + cudaFree(d_predict); + cudaFree(d_theta_exp_sqrt); +} + +void gpuPredict::predict( void ) +{ + gpuPredict::init_gpu(); + gpuPredict::compute_distance(); + gpuPredict::compute_result(); + gpuPredict::compute_error(); + gpuPredict::compute_deriv(); + gpuPredict::free_gpu(); +} + + + + + + diff --git a/gp_emulator/gpu/tests/CMakeLists.txt b/gp_emulator/gpu/tests/CMakeLists.txt new file mode 100644 index 0000000..70e9af1 --- /dev/null +++ b/gp_emulator/gpu/tests/CMakeLists.txt @@ -0,0 +1,31 @@ +find_library(CUNIT_LIBRARY NAMES cunit libcunit cunitlib) +find_path(CUNIT_INCLUDE_DIR NAMES CUnit/CUnit.h) +#mark_as_advanced(CUNIT_INCLUDE_DIR) +#MARK_AS_ADVANCED(CUNIT_LIBRARY) + + + +include_directories(${CUNIT_INCLUDE_DIR}) + +cuda_add_executable(gpu_predict_test unit_test.cu + readTestData.cu + compare_result.cu + testCdist.cu + testCublasgemm.cu + testMatrixExp.cu + testPredict.cu + testVecTimesMat.cu + testInitArray.cu + ../predict.cu + ../computeTranspose.cu + ../kernel_vectorTimesMatrix.cu + ../kernel_init_array.cu + ../kernel_cdist.cu + ../kernel_matrixExp.cu + ../kernel_elementwiseMult.cu + ../kernel_scalarMinusVec.cu + ../kernel_rowSum.cu + ../kernel_crossMinus.cu) +cuda_add_cublas_to_target(gpu_predict_test) +target_link_libraries(gpu_predict_test ${CUNIT_LIBRARY}) + diff --git a/gp_emulator/gpu/tests/compare_result.cu b/gp_emulator/gpu/tests/compare_result.cu new file mode 100644 index 0000000..cc312f6 --- /dev/null +++ b/gp_emulator/gpu/tests/compare_result.cu @@ -0,0 +1,27 @@ +#include "gpu_predict_test.h" +//#define max(a,b) (((a)>(b))?(a):(b)) + +void compare_result( const real *test_val, const real *origin_val, const int len, + const real epsilon_average, const real epsilon_max, char var_name[]) +{ + int i; + real error = 0.0 ; + real error_all = 0.0; + real error_max = 0.0; + + for( i = 0; i < len; ++i ) + { + error = abs( test_val[i] - origin_val[i] ) / abs( origin_val[i] ); + error_all += error; + error_max = max( error_max, error ); + } + error_all = error_all / len; + + if(error_all > epsilon_average || error_max > epsilon_max) + printf(" [%s] Average error = %.1e, Max error = %.1e\n",var_name, error_all, error_max); + CU_ASSERT(error_all < epsilon_average); + CU_ASSERT(error_max < epsilon_max); + +} + + diff --git a/gp_emulator/gpu/tests/gpu_predict_test.h b/gp_emulator/gpu/tests/gpu_predict_test.h new file mode 100644 index 0000000..f689678 --- /dev/null +++ b/gp_emulator/gpu/tests/gpu_predict_test.h @@ -0,0 +1,44 @@ +#include +#include +#include "../gpu_predict.h" +#include "CUnit/Basic.h" +#include + +extern real *expX, *expXsqrt, *in_train, *in_predict, *invQ, *invQt; +extern real *cdist_a, *cdist_expa; +extern real *cdist_test_var1, *cdist_test_var2, *cdist_test_var3; +extern real *error_test1; +extern real *result_py, *error_py, *deriv_py; +extern int Npredict_t, Ntrain_t, Ninputs_t; +extern int theta_size_t; + +#ifdef DOUBLE__PRECISION + #define EPSILON_AVG 1e-6 + #define EPSILON_MAX 1e-2 +#else + #define EPSILON_AVG 1e-3 + #define EPSILON_MAX 1 +#endif + + +real *readTestData(char *file_name, int size); +void compare_result(const real *test_val, const real *origin_val, + const int len, const real epsilon_average, + const real epsilon_max, char var_name[]); +void testInitArray(void); + +void testVecTimesMat(const real *c_vec, + const real *c_mat, const real *c_res, + const int vec_len, const int mat_nrows, + const int mat_ncols); +void testCdist(const real *in1, const real *in2, + const real *res, const int in1_nrows, + const int in2_nrows, const int in_ncols); +void testMatrixExp(const real *mat, + const real *res, const real alpha, + const real beta,const int size); +void testCublasgemm(const real *c_mat1, + const real *c_mat2, const real *c_res, + const int mat1_nrows, const int mat1_ncols, + const int mat2_nrows, const int mat2_ncols); +void testPredict(void); diff --git a/gp_emulator/gpu/tests/readTestData.cu b/gp_emulator/gpu/tests/readTestData.cu new file mode 100644 index 0000000..6f1938f --- /dev/null +++ b/gp_emulator/gpu/tests/readTestData.cu @@ -0,0 +1,30 @@ +#include "gpu_predict_test.h" +#include + + +real *readTestData(char *file_name, const int size) +{ + int i; + char path[400]; + sprintf(path, "data/"); + strcat(path,file_name); + + double *data_raw; + real *data; + + data_raw = (double *)malloc(sizeof(double) * size); + data = (real *)malloc(sizeof(real) * size); + FILE *file_ptr; + file_ptr = fopen(path,"rb"); + fread(data_raw, sizeof(double), size, file_ptr); + + for( i = 0; i < size; ++i ) + data[i] = (real)data_raw[i]; + + fclose(file_ptr); + return(data); + +} + + + diff --git a/gp_emulator/gpu/tests/testCdist.cu b/gp_emulator/gpu/tests/testCdist.cu new file mode 100644 index 0000000..95527c1 --- /dev/null +++ b/gp_emulator/gpu/tests/testCdist.cu @@ -0,0 +1,37 @@ +#include "gpu_predict_test.h" + +void testCdist( const real *matrix1,const real *matrix2, const real *result, + const int matrix1_nrows, const int matrix2_nrows, const int matrix_ncols ) +{ + int i; + real *matrix1_T, *matrix2_T, *gpu_result; + real *d_matrix1, *d_matrix2, *d_result; + + gpu_result = (real *)malloc( sizeof(real) * matrix1_nrows * matrix2_nrows ); + matrix1_T = computeTranspose( matrix1, matrix_ncols, matrix1_nrows ); + matrix2_T = computeTranspose( matrix2, matrix_ncols, matrix2_nrows ); + + /*GPU part*/ + for( i = 0; i < matrix1_nrows * matrix2_nrows; i++ ) + gpu_result[i] = 0; + + + cudaMalloc((void **)&d_matrix1, sizeof(real) * matrix1_nrows * matrix_ncols ); + cudaMalloc((void **)&d_matrix2, sizeof(real) * matrix2_nrows * matrix_ncols ); + cudaMalloc((void **)&d_result, sizeof(real) * matrix2_nrows * matrix1_nrows ); + cublasCheckErrors(cublasSetMatrix( matrix1_nrows, matrix_ncols, sizeof(real), matrix1_T, matrix1_nrows, d_matrix1, matrix1_nrows) ); + cublasCheckErrors(cublasSetMatrix( matrix2_nrows, matrix_ncols, sizeof(real), matrix2_T, matrix2_nrows, d_matrix2, matrix2_nrows) ); + cublasCheckErrors(cublasSetMatrix( matrix2_nrows, matrix1_nrows, sizeof(real), gpu_result, matrix2_nrows, d_result,matrix2_nrows) ); + gpu_cdist(d_matrix1, d_matrix2, d_result, matrix1_nrows, matrix_ncols, matrix2_nrows, matrix_ncols); + cudaMemcpy(gpu_result, d_result, sizeof(real) * matrix2_nrows * matrix1_nrows, cudaMemcpyDeviceToHost); + + compare_result(gpu_result, result, matrix1_nrows * matrix2_nrows, EPSILON_AVG, EPSILON_MAX, "RESULT"); + + free(matrix1_T); + free(matrix2_T); + free(gpu_result); + + cudaFree(d_matrix1); + cudaFree(d_matrix2); + cudaFree(d_result); +} diff --git a/gp_emulator/gpu/tests/testCublasgemm.cu b/gp_emulator/gpu/tests/testCublasgemm.cu new file mode 100644 index 0000000..6fc6e4b --- /dev/null +++ b/gp_emulator/gpu/tests/testCublasgemm.cu @@ -0,0 +1,48 @@ +#include "gpu_predict_test.h" + +void testCublasgemm(const real *c_matrix1, const real *c_matrix2, const real *c_result, + const int matrix1_nrows, const int matrix1_ncols, const int matrix2_nrows, + const int matrix2_ncols) +{ + cublasHandle_t handle; + cublasCreate(&handle); + + + real *d_matrix1, *d_matrix2, *d_result; + real *c_gpu_result; + c_gpu_result = (real *)malloc( sizeof(real) * matrix2_nrows * matrix2_ncols); + + cudaMalloc( (void **)&d_matrix1, sizeof(real) * matrix1_nrows * matrix1_ncols ); + cudaMalloc( (void **)&d_matrix2, sizeof(real) * matrix2_nrows * matrix2_ncols ); + cudaMalloc( (void **)&d_result, sizeof(real) * matrix1_nrows * matrix2_ncols ); + + + real *c_matrix1_T; + + + c_matrix1_T = computeTranspose(c_matrix1, matrix1_nrows, matrix1_ncols); + + cudaMemcpy( d_matrix1, c_matrix1_T, sizeof(real) * matrix1_nrows * matrix1_ncols, cudaMemcpyHostToDevice); + cudaMemcpy( d_matrix2, c_matrix2, sizeof(real) * matrix1_nrows * matrix2_ncols, cudaMemcpyHostToDevice); + + + real alpha = 1.f; + real beta = 0.f; + CUBLAS_GEMM( handle, CUBLAS_OP_N, CUBLAS_OP_T, + matrix1_nrows, matrix2_ncols, matrix1_ncols, + &alpha, + d_matrix1, matrix1_nrows, + d_matrix2, matrix2_ncols, + &beta, + d_result, matrix1_ncols ); + + cudaMemcpy( c_gpu_result, d_result, sizeof(real) * matrix1_nrows * matrix2_ncols, cudaMemcpyDeviceToHost); + c_gpu_result = computeTranspose(c_gpu_result, matrix1_nrows, matrix2_ncols); + compare_result( c_gpu_result, c_result, matrix1_nrows * matrix2_ncols, EPSILON_AVG, EPSILON_MAX, "RESULTS"); + + free(c_matrix1_T); + cudaFree(d_matrix1); + cudaFree(d_matrix2); + cudaFree(d_result); + free(c_gpu_result); +} diff --git a/gp_emulator/gpu/tests/testInitArray.cu b/gp_emulator/gpu/tests/testInitArray.cu new file mode 100644 index 0000000..dc00e62 --- /dev/null +++ b/gp_emulator/gpu/tests/testInitArray.cu @@ -0,0 +1,30 @@ +#include"gpu_predict_test.h" + +void testInitArray(void) +{ + real *c_vec, *d_vec; + real val; + int len; + int i, test, error; + c_vec = (real *)malloc( int(sizeof(real) * 1e7) ); + cudaMalloc( (void **)&d_vec, int(sizeof(real) * 1e7) ); + + for( test = 0; test < 6; ++test ) + { + val = pow( 10.0, test-5); + len = pow( 10.0, test ); + gpu_init_array(d_vec, val, len); + cudaMemcpy(d_vec, c_vec, sizeof(real) * len, cudaMemcpyDeviceToHost); + + error = 0; + for( i = 0; i < len; ++i ) + { + if( ( c_vec[i] - val ) / val > 1e-15 ) + error++; + } + CU_ASSERT(error == 0); + } + + free(c_vec); + cudaFree(d_vec); +} diff --git a/gp_emulator/gpu/tests/testMatrixExp.cu b/gp_emulator/gpu/tests/testMatrixExp.cu new file mode 100644 index 0000000..8019b79 --- /dev/null +++ b/gp_emulator/gpu/tests/testMatrixExp.cu @@ -0,0 +1,20 @@ +#include "gpu_predict_test.h" + + +void testMatrixExp(const real *mat, const real *res, const real alpha,const real beta,const int size) +{ + real *d_mat; + real *gpu_res; + + gpu_res = (real *)malloc( sizeof(real) * size ); + cudaMalloc( (void **)&d_mat, sizeof(real) * size ); + + cudaMemcpy( d_mat, mat, sizeof(real) * size, cudaMemcpyHostToDevice ); + gpu_matrixExp( d_mat, alpha, beta, size ); + cudaMemcpy( gpu_res, d_mat, sizeof(real) * size, cudaMemcpyDeviceToHost); + compare_result( gpu_res, res, size, EPSILON_AVG, EPSILON_MAX, "RESULT"); + + cudaFree(d_mat); + free(gpu_res); +} + diff --git a/gp_emulator/gpu/tests/testPredict.cu b/gp_emulator/gpu/tests/testPredict.cu new file mode 100644 index 0000000..d95ebb2 --- /dev/null +++ b/gp_emulator/gpu/tests/testPredict.cu @@ -0,0 +1,39 @@ +#include "gpu_predict_test.h" +#include + +void testPredict(void) +{ + real *gpu_result, *gpu_error, *gpu_deriv; + gpu_result = (real *)malloc(sizeof(real) * Npredict_t); + gpu_error = (real *)malloc(sizeof(real) * Npredict_t); + gpu_deriv = (real *)malloc(sizeof(real) * Npredict_t * Ninputs_t); + + real *invQ_T, *train_T, *predict_T; + invQ_T = (real *)malloc(sizeof(real) * Ntrain_t * Ntrain_t); + train_T = (real *)malloc(sizeof(real) * Npredict_t * Ninputs_t); + predict_T = (real *)malloc(sizeof(real) * Npredict_t * Ninputs_t); + + + invQ_T = computeTranspose( invQ, Ntrain_t, Ntrain_t ); + train_T = computeTranspose( in_train, Ninputs_t, Ntrain_t ); + predict_T = computeTranspose( in_predict, Ninputs_t, Npredict_t ); + + gpuPredict gpu_predict(expX, expXsqrt, invQt, invQ_T, predict_T, train_T, + gpu_result, gpu_error, gpu_deriv, Npredict_t, Ntrain_t, Ninputs_t, theta_size_t); + gpu_predict.predict(); + + gpu_deriv = computeTranspose( gpu_deriv, Npredict_t, Ninputs_t); + + compare_result( gpu_result, result_py, Npredict_t, EPSILON_AVG, EPSILON_MAX, "result"); + compare_result( gpu_error, error_py, Npredict_t, EPSILON_AVG, EPSILON_MAX, "error"); + compare_result( gpu_deriv, deriv_py, Npredict_t * Ninputs_t, EPSILON_AVG, EPSILON_MAX, "deriv"); + + free(invQ_T); + free(train_T); + free(predict_T); + cudaFree(gpu_result); + cudaFree(gpu_error); + cudaFree(gpu_deriv); +} + + diff --git a/gp_emulator/gpu/tests/testVecTimesMat.cu b/gp_emulator/gpu/tests/testVecTimesMat.cu new file mode 100644 index 0000000..ca39dbf --- /dev/null +++ b/gp_emulator/gpu/tests/testVecTimesMat.cu @@ -0,0 +1,36 @@ +#include "gpu_predict_test.h" + +void testVecTimesMat(const real *c_vec,const real *c_matrix, const real *c_res,const int vec_len, const int matrix_nrows, const int matrix_ncols) +{ + real *d_vec, *d_matrix, *d_res; + real *c_res_gpu; + real *c_matrix_T; + + CU_ASSERT (vec_len == matrix_ncols); + cudaMalloc((void **)&d_vec, sizeof(real) * vec_len ); + cudaMalloc((void **)&d_matrix, sizeof(real) * matrix_nrows * matrix_ncols ); + cudaMalloc((void **)&d_res, sizeof(real) * matrix_nrows * matrix_ncols ); + + c_matrix_T = (real *)malloc( sizeof(real) * matrix_nrows * matrix_ncols); + c_res_gpu = (real *)malloc( sizeof(real) * matrix_nrows * matrix_ncols ); + + c_matrix_T = computeTranspose(c_matrix, matrix_ncols, matrix_nrows); + + cublasCheckErrors(cublasSetVector( vec_len, sizeof(real), c_vec, 1, d_vec, 1 )); + cublasCheckErrors(cublasSetMatrix( matrix_nrows, matrix_ncols, sizeof(real), c_matrix_T, matrix_nrows, d_matrix, matrix_nrows)); + + gpu_vectorTimesMatrix( d_matrix, d_vec, d_res, matrix_nrows, matrix_ncols); + + cudaMemcpy(c_res_gpu, d_res, sizeof(real) * matrix_nrows * matrix_ncols, cudaMemcpyDeviceToHost); + c_res_gpu = computeTranspose(c_res_gpu, matrix_nrows, matrix_ncols); + compare_result( c_res_gpu, c_res, matrix_nrows * matrix_ncols, EPSILON_AVG, EPSILON_MAX, "RESULT"); + + free(c_res_gpu); + free(c_matrix_T); + + cudaFree(d_vec); + cudaFree(d_matrix); + cudaFree(d_res); + +} + diff --git a/gp_emulator/gpu/tests/unit_test.cu b/gp_emulator/gpu/tests/unit_test.cu new file mode 100644 index 0000000..07cf933 --- /dev/null +++ b/gp_emulator/gpu/tests/unit_test.cu @@ -0,0 +1,121 @@ +#include +#include +#include "gpu_predict_test.h" +#include "cuda.h" + +real *expX, *expXsqrt, *in_train, *in_predict, *invQ, *invQt; +real *cdist_a, *cdist_expa; +real *cdist_test_var1, *cdist_test_var2, *cdist_test_var3; +real *error_test1; +real *result_py, *error_py, *deriv_py; +int Npredict_t, Ntrain_t, Ninputs_t; +int theta_size_t; + +int init_suite1(void) +{ + invQ = readTestData( "invQ.bin", Ntrain_t * Ntrain_t); + invQt = readTestData( "invQt.bin", Ntrain_t ); + + expX = readTestData( "expX.bin", Ninputs_t+2); + expXsqrt = readTestData( "expXsqrt.bin", Ninputs_t ); + in_train = readTestData( "in_train.bin", Ntrain_t * Ninputs_t ); + in_predict = readTestData( "in_predict.bin", Npredict_t * Ninputs_t ); + cdist_test_var1 = readTestData( "cdist_test_var1.bin", Ntrain_t * Ninputs_t); + cdist_test_var2 = readTestData( "cdist_test_var2.bin", Npredict_t * Ninputs_t); + cdist_a = readTestData( "cdist_a.bin", Ntrain_t * Npredict_t); + cdist_expa = readTestData( "cdist_expa.bin", Ntrain_t * Npredict_t ); + + error_test1 = readTestData( "error_test1.bin", Ntrain_t * Npredict_t ); + result_py = readTestData( "result.bin", Ntrain_t * Npredict_t); + error_py = readTestData( "error.bin", Npredict_t); + deriv_py = readTestData( "deriv.bin", Ntrain_t * Npredict_t ); + + return 0; +} + +int clean_suite(void) +{ + free(expXsqrt); + free(in_train); + free(in_predict); + free(cdist_a); + free(cdist_test_var1); + free(cdist_test_var2);// have all variable been fully cleaned?? + return 0; +} + +void tests_VecTimesMat(void) +{ + testVecTimesMat(expXsqrt, in_train, cdist_test_var1, Ninputs_t, Ntrain_t, Ninputs_t ); + testVecTimesMat(expXsqrt, in_predict, cdist_test_var2, Ninputs_t, Npredict_t, Ninputs_t ); +} + +void tests_cdist(void) +{ + testCdist(cdist_test_var1,cdist_test_var2, cdist_a, Ntrain_t, Npredict_t, Ninputs_t); +} + + +void tests_matrixExp(void) +{ + testMatrixExp( cdist_a, cdist_expa, -0.5, expX[Ninputs_t], Ntrain_t * Npredict_t ); +} + +void tests_cublasgemm(void) +{ + testCublasgemm(invQ, cdist_expa, error_test1, Ntrain_t, Ntrain_t, Ntrain_t, Npredict_t); +} + + +int main(int argc, char *argv[]) +{ + if( argc != 4) + { + printf("ERROR: number of arguments is wrong (Ntrain_t, Npredict_t, Ninputs_t)\n"); + } + + Ntrain_t = atoi(argv[1]); + Npredict_t = atoi(argv[2]); + Ninputs_t = atoi(argv[3]); + theta_size_t = Ninputs_t + 2; + printf("===============================\n"); + printf("Testing with problem size:\n( npredict = %d, ntrain = %d, ninputs = %d )\n", Npredict_t, Ntrain_t, Ninputs_t); + printf("==============================="); + + CU_pSuite pSuite = NULL; + + /* initialize the CUnit test registry */ + if (CUE_SUCCESS != CU_initialize_registry()) + return CU_get_error(); + + /* add a suite to the registry */ + pSuite = CU_add_suite("Unit Tests", init_suite1, clean_suite); + + if (NULL == pSuite) { + CU_cleanup_registry(); + return CU_get_error(); + } + + + /* add the tests to the suite */ + if ((NULL == CU_add_test(pSuite, "test of gpu_vectorTimesMatrix", tests_VecTimesMat)) || + (NULL == CU_add_test(pSuite, "test of gpu_init_array", testInitArray))|| + (NULL == CU_add_test(pSuite, "test of gpu_cdist", tests_cdist))|| + (NULL == CU_add_test(pSuite, "test of gpu_MatrixExp", tests_matrixExp))|| + (NULL == CU_add_test(pSuite, "test of cublasgemm", tests_cublasgemm))|| + (NULL == CU_add_test(pSuite, "test of gpu_predict", testPredict)) + ) + { + CU_cleanup_registry(); + return CU_get_error(); + } + + //Run all tests using the CUnit Basic interface + CU_basic_set_mode(CU_BRM_VERBOSE); + CU_basic_run_tests(); + CU_cleanup_registry(); + printf("\n"); + return CU_get_error(); +} + + diff --git a/gp_emulator/multivariate_gp.py b/gp_emulator/multivariate_gp.py index fb3203b..379741c 100644 --- a/gp_emulator/multivariate_gp.py +++ b/gp_emulator/multivariate_gp.py @@ -30,7 +30,7 @@ import shutil import numpy as np -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt from GaussianProcess import GaussianProcess @@ -192,7 +192,7 @@ def compress ( self, X ): """Project full-rank vector into PC basis""" return X.dot ( self.basis_functions.T ).T - def predict ( self, y, do_deriv=True ): + def predict ( self, y, do_deriv=True, is_gpu=False ): """Prediction of input vector The individual GPs predict the PC weights, and these are used to @@ -212,7 +212,7 @@ def predict ( self, y, do_deriv=True ): if do_deriv: deriv = np.zeros ( ( y.shape[1], self.basis_functions.shape[1] ) ) for i in xrange ( self.n_pcs ): - pred_mu, pred_var, grad = self.emulators[i].predict ( y ) + pred_mu, pred_var, grad = self.emulators[i].predict ( y , is_gpu=is_gpu ) fwd += pred_mu * self.basis_functions[i] if do_deriv: deriv += np.matrix(grad).T * np.matrix(self.basis_functions[i]) diff --git a/setup.py b/setup.py index 97d52fb..4e053e8 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,69 @@ #!/usr/bin/env python - from distutils.core import setup +from setuptools import setup +import distutils.command.build as _build +import setuptools.command.install as _install + +import sys +import os +import os.path as op +import distutils.spawn as ds +import distutils.dir_util as dd + + + +def run_cmake(): + if ds.find_executable('cmake') is None: + print "CMake is required" + print "Please install cmake version >= 2.6 and re-run setup" + sys.exit(-1) + + build_dir = op.join(op.split(__file__)[0], 'build') + dd.mkpath(build_dir) + os.chdir(build_dir) + + try: + ds.spawn(['cmake','../']) + except ds.DistutilsExecError: + print "Error while running cmake" + sys.exit(-1) + try: + ds.spawn(['make','-j']) + except ds.DistutilsExecError: + print "Error while compiling" + sys.exit(-1) + + +class install(_install.install): + def run(self): + cwd = os.getcwd() + run_cmake() + os.chdir(cwd) + _install.install.run(self) + +class build(_build.build): + def run(self): + cwd = os.getcwd() + run_cmake() + os.chdir(cwd) + _build.build.run(self) + +class benchmark(_install.install): + def run(self): + os.system("python tests/benchmark.py") + +class test(_install.install): + def run(self): + cwd = os.getcwd() + testdir = op.join(op.split(__file__)[0], 'tests') + os.chdir(testdir) + os.system("python unit_tests.py") + os.chdir(cwd) + + + + setup(name='gp_emulator', version='1.4.3', description='A Python GaussianProcess emulator software package', @@ -9,4 +71,5 @@ author_email='j.gomez-dans@ucl.ac.uk', url='http://bitbucket.org/gomezdansj/gp_emulator', packages=['gp_emulator'], + cmdclass={'build':build, 'install':install, 'test':test, 'benchmark':benchmark}, ) diff --git a/tests/benchmark.py b/tests/benchmark.py new file mode 100644 index 0000000..72b5a1b --- /dev/null +++ b/tests/benchmark.py @@ -0,0 +1,62 @@ +#!/usr/local/bin/python + +import numpy as np +import scipy.spatial.distance as dist +import _gpu_predict +import time +from gp_emulator import GaussianProcess +from types import MethodType +import sys + +def set_testing_val (self, Ninputs, Npredict, Ntrain): + self.D = Ninputs + self.theta = np.random.random((Ninputs+2)) + self.invQ = np.random.random(( Ntrain, Ntrain )) + self.invQt = np.random.random ((Ntrain)) + + +if __name__ == '__main__': + + GaussianProcess.set_testing_val = MethodType(set_testing_val, None, GaussianProcess) + print 'Problem_size\tCPU time\tGPU time\tSpeedup\tStatus' + print '-----------------------------' + + for Npredict in xrange(np.int(1e5), np.int(1e6), np.int(1e5)): + Ntrain = 250 + Ninputs = 10 + + inputs = np.random.random(( Ntrain, Ninputs)) + testing = np.random.random(( Npredict, Ninputs)) + + gp = GaussianProcess(inputs, []) + gp.set_testing_val(Ninputs, Npredict, Ntrain) + + #CPU predict + start = time.time() + [mu_c, var_c, deriv_c] = gp.predict(testing, is_gpu=False ) + end = time.time() + cputime = end -start + + #GPU predict + start = time.time() + [mu_g, var_g, deriv_g] = gp.predict(testing, is_gpu = True, precision = np.float32, threshold = 1e5) + end =time.time() + gputime = end - start + print "%d\t%.2fs\t%.2fs\t%.2fx\t" % (Npredict, cputime, gputime, cputime/gputime), + + + + # checking results + try: + e_mu = max(abs(mu_c - mu_g) ) / np.max(abs(mu_c)) + e_var = max(abs(var_c - var_g)) / np.max(abs(var_c)) + e_deriv = np.max(abs(deriv_c - deriv_g)) / np.max(abs(deriv_c)) + except ValueError: + print 'Results have invalid data type or dimension.' + if e_mu > 1e-5 or e_var > 1e-5 or e_deriv > 1e-5: + print 'FAILED\t', + else: + print 'Pass\t', + print 'e_mu=%.2g\te_var=%.2g\te_deriv=%.2g\t'%(e_mu, e_var, e_deriv) + + diff --git a/tests/benchmark2.py b/tests/benchmark2.py new file mode 100644 index 0000000..18cc4a0 --- /dev/null +++ b/tests/benchmark2.py @@ -0,0 +1,123 @@ +#!/usr/local/bin/python + +import numpy as np +import scipy.spatial.distance as dist +import _gpu_predict +import time +from gp_emulator import GaussianProcess +from types import MethodType +import sys + +def set_testing_val (self, Ninputs, Npredict, Ntrain): + self.D = Ninputs + self.theta = np.random.random((Ninputs+2)) + self.invQ = np.random.random(( Ntrain, Ntrain )) + self.invQt = np.random.random ((Ntrain)) + + +def predict_benchmark( self, testing, do_unc=True ): + start_predict = time.time() + ( nn, D ) = testing.shape + assert D == self.D + expX = np.exp ( self.theta ) + + start_mvel = time.time() + in_dist1 = np.sqrt(expX[:(self.D)])*self.inputs + in_dist2 = np.sqrt(expX[:(self.D)])*testing + end_mvel = time.time() + print 'mvel',2, end_mvel - start_mvel + + start_dist = time.time() + a = dist.cdist ( in_dist1, in_dist2, 'sqeuclidean') + end_dist = time.time() + print 'dist',1, end_dist - start_dist + + start_exp = time.time() + a = expX[self.D]*np.exp(-0.5*a) + end_exp = time.time() + print 'exp-ele',3, end_exp - start_exp + + b = expX[self.D] + + start_mv =time.time() + mu = np.dot( a.T, self.invQt) + end_mv = time.time() + print 'mv',1,end_mv - start_mv + + if do_unc: + start_mm = time.time() + k = np.dot(self.invQ,a) + end_mm = time.time() + print 'mm',1,end_mm - start_mm + + start_mulele = time.time() + k = a * k + end_mulele = time.time() + print 'mul-ele',1, end_mulele - start_mulele + + start_sum = time.time() + var = b - np.sum (b, axis=0) + end_sum = time.time() + print 'sum',1, end_sum - start_sum + # Derivative and partial derivatives of the function + deriv = np.zeros ( ( nn, self.D ) ) + + crossminus = 0 + mul_ele2 = 0 + mv2 = 0 + + for d in xrange ( self.D ): + s_crossminus = time.time() + aa = self.inputs[:,d].flatten()[None,:] - testing[:,d].flatten()[:,None] + e_crossminus = time.time() + crossminus = crossminus + e_crossminus - s_crossminus + + s_mul_ele2 = time.time() + c = a*aa.T + e_mul_ele2 = time.time() + mul_ele2 = mul_ele2 + e_mul_ele2 - s_mul_ele2 + + s_mv2 = time.time() + deriv[:, d] = expX[d]*np.dot(c.T, self.invQt) + e_mv2 = time.time() + mv2 = mv2 + e_mv2 - s_mv2 + end_predict = time.time() + + print 'crossminus',10, crossminus + print 'mul_ele2',10,mul_ele2 + print 'mv2',10, mv2 + print 'predict',1,end_predict - start_predict + + if do_unc: + return mu, var, deriv + else: + return mu, deriv + + +if __name__ == '__main__': + + GaussianProcess.set_testing_val = MethodType(set_testing_val, None, GaussianProcess) + GaussianProcess.predict_benchmark = MethodType(predict_benchmark, None, GaussianProcess) + print 'Problem_size\tCPU time\tGPU time\tSpeedup\tStatus' + print '-----------------------------' + + for i in xrange(10): + + Npredict = np.int(1e6) + Ntrain = 250 + Ninputs = 10 + + inputs = np.random.random(( Ntrain, Ninputs)) + testing = np.random.random(( Npredict, Ninputs)) + + gp = GaussianProcess(inputs, []) + gp.set_testing_val(Ninputs, Npredict, Ntrain) + + #CPU predict + start = time.time() + [mu_c, var_c, deriv_c] = gp.predict_benchmark(testing ) + end = time.time() + cputime = end -start + print 'cputime', cputime + print '==================' + diff --git a/tests/data/.gitignore b/tests/data/.gitignore new file mode 100644 index 0000000..5e7d273 --- /dev/null +++ b/tests/data/.gitignore @@ -0,0 +1,4 @@ +# Ignore everything in this directory +* +# Except this file +!.gitignore diff --git a/tests/test_emulator.py b/tests/test_emulator.py new file mode 100644 index 0000000..16a7a49 --- /dev/null +++ b/tests/test_emulator.py @@ -0,0 +1,67 @@ +import numpy as np +#import matplotlib.pyplot as plt +from gp_emulator import MultivariateEmulator, GaussianProcess +from types import MethodType +import scipy.spatial.distance as dist +import os + + +#d = np.loadtxt("../data/test_00100.txt", delimiter=",") +d = np.loadtxt("../data/validation_set_100000.txt",delimiter=",") +gp = MultivariateEmulator(dump="../data/prosail_vza30_sza0_saa0_vaa0_n250.npz") + + +wv = np.arange (400, 2501) +#for isample in [0, 10, 20, 50]: +# plt.plot ( wv, d[isample, 10:], '-', lw=2) +# r = gp.predict(d[isample,:10])[0] +# plt.plot ( wv, d[isample, 10:], '--', lw=2) + +def perband_emulators ( emulators, band_pass ): + """This function creates per band emulators from the full-spectrum + emulator. Should be faster in many cases""" + + n_bands = band_pass.shape[0] + x_train_pband = [ emulators.X_train[:,band_pass[i,:]].mean(axis=1) \ + for i in xrange( n_bands ) ] + x_train_pband = np.array ( x_train_pband ) + emus = [] + # add get_testing_val + #GaussianProcess.get_testing_val = MethodType(get_testing_val, None, GaussianProcess) + + for i in xrange( n_bands ): + gp = GaussianProcess ( emulators.y_train[:]*1, \ + x_train_pband[i,:] ) + print 'p1', (emulators.y_train[:]*1).shape, 'p2', x_train_pband[i,:].shape + gp.learn_hyperparameters ( n_tries=5 ) + emus.append ( gp ) + return emus + + +precision = "float32" +new_output = d[:,10:][:,452:486].mean(axis=1) +band_pass = np.zeros( (1,2101), dtype=np.bool) +band_pass[:, 442:476] = 1 +gp_single = perband_emulators ( gp, band_pass ) + + +print isinstance ( gp_single[0], GaussianProcess) +X = d[:, :10] +gpur = gp_single[0].predict ( X ,is_gpu=True, prec = precision) +cpur = gp_single[0].predict ( X, is_gpu=False) + + +for i in xrange(3): + val_absolute = np.mean(np.abs(cpur[i])) + error_absolute = np.abs( cpur[i] - gpur[i] ) + error_relative = error_absolute / np.abs( cpur[i] ) + print '\n--------------' + print 'output', i, precision + print 'absolute_val\t%.2g'%val_absolute + print 'error_absolute\t', '%.2g [average]\t%.2g [max]'% (np.mean(error_absolute), np.max(error_absolute)) + print 'error_relative\t', '%.2g [average]\t%.2g [max]'% (np.mean(error_relative), np.max(error_relative)) + + + + + diff --git a/tests/test_perband_emulator.py b/tests/test_perband_emulator.py index 5baeed4..64279ea 100644 --- a/tests/test_perband_emulator.py +++ b/tests/test_perband_emulator.py @@ -4,16 +4,18 @@ -#d = np.loadtxt("test_00100.txt", delimiter=",") -d = np.loadtxt("data/validation_set_001000.txt",delimiter=",") -gp = MultivariateEmulator(dump="data/prosail_30_0_30_0.npz") +#d = np.loadtxt("../data/test_00100.txt", delimiter=",") +d = np.loadtxt("../data/validation_set_001000.txt",delimiter=",") +gp = MultivariateEmulator(dump="../data/prosail_30_0_30_0.npz") wv = np.arange (400, 2501) -for isample in [0, 10, 20, 50]: - plt.plot ( wv, d[isample, 10:], '-', lw=2) - r = gp.predict(d[isample,:10])[0] - plt.plot ( wv, d[isample, 10:], '--', lw=2) +#for isample in [0, 10, 20, 50]: +# plt.plot ( wv, d[isample, 10:], '-', lw=2) +# r = gp.predict(d[isample,:10])[0] +# plt.plot ( wv, d[isample, 10:], '--', lw=2) + + @@ -25,13 +27,11 @@ def perband_emulators ( emulators, band_pass ): x_train_pband = [ emulators.X_train[:,band_pass[i,:]].mean(axis=1) \ for i in xrange( n_bands ) ] x_train_pband = np.array ( x_train_pband ) - print x_train_pband.shape - print emulators.y_train.shape - print n_bands emus = [] for i in xrange( n_bands ): gp = GaussianProcess ( emulators.y_train[:]*1, \ x_train_pband[i,:] ) + print 'p1', (emulators.y_train[:]*1).shape, 'p2', x_train_pband[i,:].shape gp.learn_hyperparameters ( n_tries=5 ) emus.append ( gp ) return emus @@ -46,9 +46,15 @@ def perband_emulators ( emulators, band_pass ): print isinstance ( gp_single[0], GaussianProcess) X = d[:, :10] -r = gp_single[0].predict ( X )[0] -plt.plot ( new_output, r,'o') -plt.plot([0, 0.6], [0, 0.6], 'k--') +gpur = gp_single[0].predict ( X ,is_gpu=True)[0] +cpur = gp_single[0].predict ( X, is_gpu=False)[0] +print gpur - cpur +print 'predict finish' +#plt.plot ( new_output, r,'o') +#plt.plot([0, 0.6], [0, 0.6], 'k--') + +#plt.show() +#%timeit r = gp_single[0].predict ( X )[0] + -#%timeit r = gp_single[0].predict ( X )[0] \ No newline at end of file diff --git a/tests/unit_tests.py b/tests/unit_tests.py new file mode 100755 index 0000000..9933a9b --- /dev/null +++ b/tests/unit_tests.py @@ -0,0 +1,91 @@ +import numpy as np +import scipy.spatial.distance as dist +import time +import os as os +class GP: + def __init__ (self, P, N, M): + self.D = P + D=P + self.theta = np.random.random((P+2)) + self.inputs = np.random.random ((M,P)) + self.invQ = np.random.random((M,M)) + self.invQt = np.random.random ((M)) + + def get_testing_val (self, testing, do_unc = True): + + ( nn, D ) = testing.shape + assert D == self.D + expX=np.exp(self.theta) + + t_expX = expX + t_expXsqrt = np.sqrt(expX[:(self.D)]) + t_inputs = self.inputs + t_testing = testing + t_invQ = self.invQ + + cdist_test_var_1= t_expXsqrt * t_inputs + cdist_test_var_2 = np.sqrt(expX[:(self.D)])*testing + cdist_test_var_3 = dist.cdist ( np.sqrt(expX[:(self.D)])*self.inputs, np.sqrt(expX[:(self.D)])*testing, 'sqeuclidean') + cdist_test_var_4 = expX[self.D]*np.exp(-0.5*cdist_test_var_3) + + mu = np.dot( cdist_test_var_4.T, self.invQt) + var_test_1 = np.dot(self.invQ, cdist_test_var_4) + + + + t_expX.tofile("./data/expX.bin") + t_expXsqrt.tofile("./data/expXsqrt.bin") + t_inputs.tofile("./data/in_train.bin") + t_testing.tofile("./data/in_predict.bin") + cdist_test_var_1.tofile("./data/cdist_test_var1.bin") + cdist_test_var_2.tofile("./data/cdist_test_var2.bin") + cdist_test_var_3.tofile("./data/cdist_a.bin") + cdist_test_var_4.tofile("./data/cdist_expa.bin") + + t_invQ.tofile("./data/invQ.bin") + mu.tofile("./data/result.bin") + var_test_1.tofile("./data/error_test1.bin") + self.invQt.tofile("./data/invQt.bin") + + + ( nn, D ) = testing.shape + assert D == self.D + + expX = np.exp ( self.theta ) + + a = dist.cdist ( np.sqrt(expX[:(self.D)])*self.inputs, np.sqrt(expX[:(self.D)])*testing, 'sqeuclidean') + a = expX[self.D]*np.exp(-0.5*a) + b = expX[self.D] + + mu = np.dot( a.T, self.invQt) + if do_unc: + var = b - np.sum ( a * np.dot(self.invQ,a), axis=0) + # Derivative and partial derivatives of the function + deriv = np.zeros ( ( nn, self.D ) ) + + for d in xrange ( self.D ): + aa = self.inputs[:,d].flatten()[None,:] - testing[:,d].flatten()[:,None] + c = a*aa.T + deriv[:, d] = expX[d]*np.dot(c.T, self.invQt) + var.tofile("./data/error.bin") + deriv.tofile("./data/deriv.bin") + + if do_unc: + return mu, var, deriv + else: + return mu, deriv + + +if __name__ == '__main__': + for size in xrange(int(12345), int(1e5),int(6e4)): + N=size + M=250 + P=10 + + testing=np.random.random((N,P)) + gp=GP(P,N,M) + + gp.get_testing_val(testing) + command = "./gpu_predict_test %d %d %d" % (M, N, P) + os.system(command) +