Python筆記_第四篇_高階編程_進程、線程、協程_5.GPU加速

Numba:高性能計算的高生產率
　　在這篇文章中，筆者將向你介紹一個來自Anaconda的Python編譯器Numba，它能夠在CUDA-capable GPU或多核cpu上編譯Python代碼。Python一般不是一種編譯語言，你可能想知道爲何要使用Python編譯器。答案固然是：運行本地編譯的代碼要比運行動態的、解譯的代碼快不少倍。Numba容許你爲Python函數指定類型簽名，從而在運行時啓用編譯(這就是「Just-in-Time」，即時，也能夠說JIT編譯)。Numba動態編譯代碼的能力意味着你不會所以而拋棄Python的靈活性。這是向提供高生產率編程和高性能計算的完美結合邁出的一大步。

　　使用Numba能夠編寫標準的Python函數，並在CUDA-capable GPU上運行它們。Numba是爲面向數組的計算任務而設計的，很像你們經常使用的NumPy庫。在面向數組的計算任務中，數據並行性對於像GPU這樣的加速器是很天然的。Numba瞭解NumPy數組類型，並使用它們生成高效的編譯代碼，用於在GPU或多核CPU上執行。所需的編程工做能夠很簡單，就像添加一個函數修飾器來指示Numba爲GPU編譯同樣。例如，在下面的代碼中，@ vectorize decorator會生成一個編譯的、矢量化的標量函數在運行時添加的版本，這樣它就能夠用於在GPU上並行處理數據數組。　　

　　要在CPU上編譯和運行相同的函數，咱們只需將目標更改成「CPU」，它將在編譯水平上帶來性能，在CPU上向量化C代碼。這種靈活性能夠幫助你生成更可重用的代碼，並容許你在沒有GPU的機器上開發。

import numpy as np

from numba import vectorize

@vectorize(['float32(float32, float32)'], target='cuda')

def Add(a, b):

return a + b

# Initialize arrays

N = 100000
A = np.ones(N, dtype=np.float32)
B = np.ones(A.shape, dtype=A.dtype)
C = np.empty_like(A, dtype=A.dtype)
# Add arrays on GPU
C = Add(A, B)

　　關於Python 的GPU-Accelerated庫

　　CUDA並行計算平臺的優點之一是其可用的GPU加速庫的闊度。Numba團隊的另外一個項目叫作pyculib，它提供了一個Python接口，用於CUDA cuBLAS(dense linear algebra，稠密線性代數)，cuFFT(Fast Fourier Transform，快速傅里葉變換)，和cuRAND(random number generation，隨機數生成)庫。許多應用程序都可以經過使用這些庫得到顯著的加速效果，而不須要編寫任何特定於GPU的代碼。例如，下面的代碼使用「XORWOW」僞隨機數生成器在GPU上生成100萬個均勻分佈的隨機數。

　　import numpy as np

　　from pyculib import rand as curand

　　prng = curand.PRNG(rndtype=curand.PRNG.XORWOW)

　　rand = np.empty(100000)

　　prng.uniform(rand)

　　print rand[:10]
　　CUDA Python的高並行性

　　Anaconda(原名Continuum Analytics)認識到，在某些計算上實現大的速度須要一個更具表現力的編程接口，它比庫和自動循環矢量化更詳細地控制並行性。所以，Numba有另外一組重要的特性，構成了其非正式名稱「CUDA Python」。Numba公開了CUDA編程模型，正如CUDA C/ C++，可是使用純python語法，這樣程序員就能夠建立自定義、調優的並行內核，而不會放棄python帶來的便捷和優點。Numba的CUDA JIT(經過decorator或函數調用可用)在運行時編譯CUDA Python函數，專門針對你所使用的類型，它的CUDA Python API提供了對數據傳輸和CUDA流的顯式控制，以及其餘特性。

　　下面的代碼示例演示了一個簡單的Mandelbrot設置內核。請注意，mandel_kernel函數使用Numba提供的cuda.threadIdx，cuda.blockIdx，cuda.blockDim和cuda.gridDim架構來計算當前線程的全局X和Y像素索引。與其餘CUDA語言同樣，咱們經過插入在括號內一個「執行配置」(CUDA-speak用於線程數和線程塊啓動內核)，在函數名和參數列表之間中: mandel_kernel[griddim, blockdim](-2.0, 1.0, -1.0, 1.0, d_image, 20)。你還能夠看到使用to_host和to_device API函數來從GPU中複製數據。

　　Mandelbrot的例子將在Github上持續更新。　　

@cuda.jit(device=True)

　　def mandel(x, y, max_iters):

　　"""
　　Given the real and imaginary parts of a complex number,
　　determine if it is a candidate for membership in the Mandelbrot
　　set given a fixed number of iterations.
　　"""

　　c = complex(x, y)

　　z = 0.0j

　　for i in range(max_iters):

　　z = z*z + c

　　if (z.real*z.real + z.imag*z.imag) >= 4:

　　return i

　　return max_iters

　　@cuda.jit

　　def mandel_kernel(min_x, max_x, min_y, max_y, image, iters):

　　height = image.shape[0]

　　width = image.shape[1]

　　pixel_size_x = (max_x - min_x) / width

　　pixel_size_y = (max_y - min_y) / height

　　startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x

　　startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y

　　gridX = cuda.gridDim.x * cuda.blockDim.x;

　　gridY = cuda.gridDim.y * cuda.blockDim.y;

　　for x in range(startX, width, gridX):

　　real = min_x + x * pixel_size_x

　　for y in range(startY, height, gridY):

　　imag = min_y + y * pixel_size_y

　　image[y, x] = mandel(real, imag, iters)

　　gimage = np.zeros((1024, 1536), dtype = np.uint8)

　　blockdim = (32, 8)

　　griddim = (32,16)

　　start = timer()

　　d_image = cuda.to_device(gimage)

　　mandel_kernel[griddim, blockdim](-2.0, 1.0, -1.0, 1.0, d_image, 20)

　　d_image.to_host()

　　dt = timer() - start

　　print "Mandelbrot created on GPU in %f s" % dt

　　imshow(gimage)
　　在一臺帶有NVIDIA Tesla P100 GPU和Intel Xeon E5-2698 v3 CPU的服務器上，這個CUDA Python Mandelbrot代碼運行的速度比純Python版本快1700倍。1700倍彷佛有些不切實際，但請記住，咱們正在比較編譯的、並行的、GPU加速的Python代碼來解釋CPU上的單線程Python代碼。

 

　　今天開始使用Numba吧

　　Numba爲Python開發人員提供了一種簡單的進入GPU加速計算的方法，並提供了一種使用愈來愈複雜的CUDA代碼的方法，其中至少有新語法和術語。你能夠從簡單的函數decorator開始實現自動編譯函數，或者使用pyculib的強大的CUDA庫。當你提升對並行編程概念的理解時，當你須要對於並行線程的有表現力且靈活的控制時，CUDA能夠在不須要你第一天就徹底瞭解的狀況下使用。

　　Numba是一個BSD認證的開源項目，它自己嚴重依賴於LLVM編譯器的功能。Numba的GPU後端使用了基於LLVM的NVIDIA編譯器SDK。CUDA庫周圍的pyculib包裝器也是開源且通過BSD認證的。

　　要開始使用Numba，第一步是下載並安裝Anaconda Python發行版，這是一個「徹底免費的、用於大規模數據處理、預測分析和科學計算的Python發行版」，其中包括許多流行的軟件包(Numpy、Scipy、Matplotlib、iPython等)和「conda」，這是一個強大的包管理器。一旦您安裝了Anaconda，經過鍵入conda安裝numba cudatoolkit pyculib，安裝所需的CUDA包。而後在ContinuumIO github存儲庫中查看CUDA的Numba教程。筆者建議你在Anaconda的博客上查看Numba的帖子。

 

Nvidia的CUDA 架構爲咱們提供了一種便捷的方式來直接操縱GPU 並進行編程，可是基於 
C語言的CUDA實現較爲複雜，開發週期較長。而python 做爲一門普遍使用的語言，具備 
簡單易學、語法簡單、開發迅速等優勢。做爲第四種CUDA支持語言，相信python必定會 
在高性能計算上有傑出的貢獻–pyCUDA。

pyCUDA特色

pyCUDA工做流程

調用的基本例子

包含內容

pyCUDA特色
CUDA徹底的python實現

編碼更爲靈活、迅速、自適應調節代碼

更好的魯棒性，自動管理目標生命週期和錯誤檢測

包含易用的工具包，包括基於GPU的線性代數庫、reduction和scan，添加了快速傅里葉變換包和線性代數包LAPACK

完整的幫助文檔Wiki

pyCUDA的工做流程
具體的調用流程以下：開始
編寫python程序
python程序檢查？
調用pyCUDA編譯CUDA 源碼並上傳GPU
編譯正確？
PyCUDA’s numpy進行數據讀入處理
數據讀入處理成功？
輸出GPU 加速處理結果
結束

調用基本例子

import pycuda.autoinit
import pycuda.driver as drv
import numpy

from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void multiply_them(float *dest, float *a, float *b)
{
  const int i = threadIdx.x;
  dest[i] = a[i] * b[i];
}
""")

multiply_them = mod.get_function("multiply_them")

a = numpy.random.randn(400).astype(numpy.float32)
b = numpy.random.randn(400).astype(numpy.float32)

dest = numpy.zeros_like(a)
multiply_them(
        drv.Out(dest), drv.In(a), drv.In(b),
        block=(400,1,1), grid=(1,1))

print dest-a*b
補充內容：對於GPU 加速python還有功能包，例如處理圖像的pythonGPU加速包—— pyGPU

以及專門的GPU 加速python機器學習包—— scikitCUDA

GPU入門
如今GPU編程正變得愈來愈流行，因爲CPU串行執行的侷限性，程序若是在GPU上運行，則可真正作到多線程，併發執行，極大減小運行時間，這對於分秒必爭的科學計算，以及新興的人工智能領域都帶來了極大的便利。 目前，GPU編程以NVIDIA的CUDA平臺爲主，支持四種語言C、C++、Fortran(PGI)以及Python。目前CUDA的最新版本已經達到7.5，具體配置能夠看官方指導和其它教程，這裏不作具體介紹。 下面咱們具體來看Python的GPU編程。 個人顯卡是GeForce GT 740M,安裝CUDA7.5,使用Python2.7搭配相關庫。 首先咱們要引入一些必要的包 from numbapro import cuda 是cuda包是必須導入的，不然不能使用GPU。 引入以後就能夠調用cuda對象了。例如，建立一個一維網格

tx = cuda.threadIdx.x
bx = cuda.blockIdx.x
bw = cuda.blockDim.x
i = tx + bx * bw
array[i] = something(i)
也能夠簡化成

i = cuda.grid(1)
array[i] = something(i)
上面兩段代碼實現的功能是同樣的。

接下來咱們瞭解一下CUDA Stream的操做 CUDA流是對CUDA設備的命令隊列，經過特殊的流，CUDA API能夠變爲異步，這也意味着請求肯在隊列結束前返回。存儲器傳送指令和內核調用均可以使用CUDA流。 下面咱們來看示例：

stream = cuda.stream()
devary = cuda.to_device(an_array, stream=stream)
a_cuda_kernel[griddim, blockdim, stream](devary)
cuda.copy_to_host(an_array, stream=stream)
# 在an_array中的數據可能還沒有就緒
stream.synchronize()
# an_array中的數據已經就緒
另外一種語法是使用Python環境

stream = cuda.stream()
with stream.auto_synchronize():
devary = cuda.to_device(an_array, stream=stream)
a_cuda_kernel[griddim, blockdim, stream](devary)
devary.copy_to_host(an_array, stream=stream)
# an_array中的數據已經就緒
接下來是關於共享內存： 爲了達到最大性能，CUDA內核須要使用共享內存用於緩存數據，CUDA編譯器支持使用cuda.shared.array(shape, dtype)方法用來指定使用內核中的對象。 下面看一個例子

bpg = 50
tpb = 32
n = bpg * tpb

@jit(argtypes=[float32[:,:], float32[:,:], float32[:,:]], target='gpu')
def cu_square_matrix_mul(A, B, C):
sA = cuda.shared.array(shape=(tpb, tpb), dtype=float32)
sB = cuda.shared.array(shape=(tpb, tpb), dtype=float32)

tx = cuda.threadIdx.x
ty = cuda.threadIdx.y
bx = cuda.blockIdx.x
by = cuda.blockIdx.y
bw = cuda.blockDim.x
bh = cuda.blockDim.y

x = tx + bx * bw
y = ty + by * bh

acc = 0.
for i in range(bpg):
if x < n and y < n:
sA[ty, tx] = A[y, tx + i * tpb]
sB[ty, tx] = B[ty + i * tpb, x]

cuda.syncthreads()

if x < n and y < n:
for j in range(tpb):
acc += sA[ty, j] * sB[j, tx]

cuda.syncthreads()

if x < n and y < n:
C[y, x] = acc
以上就是python GPU編程的入門介紹。

咱們以FFT 爲例，
看看究竟如何利用GPU進行加速。 先看示例代碼，而後進行講解。

import sys
import numpy as np
from scipy.signal import fftconvolve
from scipy import misc, ndimage
from matplotlib import pyplot as plt
from numbapro.cudalib import cufft
from numbapro import cuda, vectorize
from timeit import default_timer as timer

@vectorize(['complex64(complex64, complex64)'], target='gpu')
#目標平臺是64位機器且擁有GPU
def vmult(a, b):
return a * b

def best_grid_size(size, tpb):
bpg = np.ceil(np.array(size, dtype=np.float) / tpb).astype(np.int).tolist()
return tuple(bpg)

def main():
# 構建過濾器
laplacian_pts = '''
-4 -1 0 -1 -4
-1 2 3 2 -1
0 3 4 3 0
-1 2 3 2 -1
-4 -1 0 -1 -4
'''.split()

laplacian = np.array(laplacian_pts, dtype=np.float32).reshape(5, 5)

# 構建圖像
try:
filename = sys.argv[1]
image = ndimage.imread(filename, flatten=True).astype(np.float32)
except IndexError:
image = misc.lena().astype(np.float32)

print("Image size: %s" % (image.shape,))

response = np.zeros_like(image)
response[:5, :5] = laplacian

# CPU
ts = timer()
cvimage_cpu = fftconvolve(image, laplacian, mode='same')
te = timer()
print('CPU: %.2fs' % (te - ts))

# GPU
threadperblock = 32, 8
blockpergrid = best_grid_size(tuple(reversed(image.shape)), threadperblock)
print('kernel config: %s x %s' % (blockpergrid, threadperblock))

# cuFFT系統，觸發器初始化.
# 對於小數據集來講，效果更明顯.
# 不該該計算這裏浪費的時間
cufft.FFTPlan(shape=image.shape, itype=np.complex64, otype=np.complex64)

# 開始GPU運行計時
ts = timer()
image_complex = image.astype(np.complex64)
response_complex = response.astype(np.complex64)

d_image_complex = cuda.to_device(image_complex)
d_response_complex = cuda.to_device(response_complex)

cufft.fft_inplace(d_image_complex)
cufft.fft_inplace(d_response_complex)

vmult(d_image_complex, d_response_complex, out=d_image_complex)

cufft.ifft_inplace(d_image_complex)

cvimage_gpu = d_image_complex.copy_to_host().real / np.prod(image.shape)

te = timer()
print('GPU: %.2fs' % (te - ts))

# 繪製結果
plt.subplot(1, 2, 1)
plt.title('CPU')
plt.imshow(cvimage_cpu, cmap=plt.cm.gray)
plt.axis('off')

plt.subplot(1, 2, 2)
plt.title('GPU')
plt.imshow(cvimage_gpu, cmap=plt.cm.gray)
plt.axis('off')

plt.show()

if __name__ == '__main__':
main()
看看運行結果：

 

時間對比：

Image size: (512L, 512L)
CPU: 0.66s
kernel config: (16, 64) x (32, 8)
GPU: 0.09s
[Finished in 61.4s]
1.導入的包

Numpy提供常見的數學函數，包含許多有用的數學庫

Scipy是python下的數值計算庫，和Numpy同樣是科學計算不可或缺的庫

Matplotlib是用以繪製二維圖形的Python模塊

Numbapro是CUDA提供的專用庫

timeit是計時工具

2.首先程序提供了一個測試數據集並進行轉化，對於CPU部分，直接調用Scipy中的fftconvolve函數計算出結果，而GPU則主要調用了numbapro中cufft庫，具體使用參考官方文檔。

3.計算出結果後繪圖展現，主要使用matplotlib的方法，只需設置plt對象並傳入參數便可。

以上在代碼中有詳細註釋，不同展開。總的來講，使用GPU仍是得參考官方的文檔，並結合Python編程，才能解決複雜的問題並達到良好的效果。

第三篇
相信若是你使用過Python Numpy包，必定了解NumPy（Numeric Python）提供了許多高級的數值編程工具，如：矩陣數據類型、矢量處理，以及精密的運算庫。它專爲進行嚴格的數字處理而產生。多爲不少大型金融公司使用，以及核心的科學計算組織如：Lawrence Livermore，NASA用其處理一些原本使用C++，Fortran或Matlab等所作的任務。 　　可是因爲複雜的計算，Numpy的計算效率不免受到影響，所以咱們對它進行了許多優化，用於優化的包有PyPy、Numba 與 Cython，而NumbaPro就是創建在Numba和cuda基礎上的高級優化方法。 　　下面咱們一塊兒來看。 　　使用NumbaPro,咱們能夠對Numpy中的方法進行優化，使Python代碼能夠動態編譯爲機器碼，並在運行中加載，使得GPU充分發揮多線程的優點。針對GPU，Numbapro也能夠自動完成工做，並優化GPU體系結構的代碼。另外，基於CUDA API編寫的Python代碼也能夠有效地利用硬件。 　　說了這麼多，下面就讓咱們從簡單的示例開始學習。

from numbapro import vectorize
@vectorize(['float32(float32, float32)'], target='cpu')
def sum(a, b):
return a + b
　　若是須要使用GPU來運行，只須要將第二行改爲@vectorize(['float32(float32, float32)'], target='gpu')

　　對於更復雜的操做，可使用Just-In-Time (JIT)來編譯。

from numbapro import cuda

@cuda.jit('void(float32[:], float32[:], float32[:])')
def sum(a, b, result):
i = cuda.grid(1) # 等價於threadIdx.x + blockIdx.x * blockDim.x
result[i] = a[i] + b[i]

# 調用: sum[grid_dim, block_dim](big_input_1, big_input_2, result_array)
　　下面繼續看一個具體的應用：

import numpy as np
import math
import time
from numba import *
from numbapro import cuda
from blackscholes_numba import black_scholes, black_scholes_numba
#import logging; logging.getLogger().setLevel(0)

RISKFREE = 0.02
VOLATILITY = 0.30

A1 = 0.31938153
A2 = -0.356563782
A3 = 1.781477937
A4 = -1.821255978
A5 = 1.330274429
RSQRT2PI = 0.39894228040143267793994605993438

@cuda.jit(argtypes=(double,), restype=double, device=True, inline=True)
def cnd_cuda(d):
K = 1.0 / (1.0 + 0.2316419 * math.fabs(d))
ret_val = (RSQRT2PI * math.exp(-0.5 * d * d) *
(K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5))))))
if d > 0:
ret_val = 1.0 - ret_val
return ret_val

@cuda.jit(argtypes=(double[:], double[:], double[:], double[:], double[:],
double, double))
def black_scholes_cuda(callResult, putResult, S, X,
T, R, V):
# S = stockPrice
# X = optionStrike
# T = optionYears
# R = Riskfree
# V = Volatility
i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x
if i >= S.shape[0]:
return
sqrtT = math.sqrt(T[i])
d1 = (math.log(S[i] / X[i]) + (R + 0.5 * V * V) * T[i]) / (V * sqrtT)
d2 = d1 - V * sqrtT
cndd1 = cnd_cuda(d1)
cndd2 = cnd_cuda(d2)

expRT = math.exp((-1. * R) * T[i])
callResult[i] = (S[i] * cndd1 - X[i] * expRT * cndd2)
putResult[i] = (X[i] * expRT * (1.0 - cndd2) - S[i] * (1.0 - cndd1))

def randfloat(rand_var, low, high):
return (1.0 - rand_var) * low + rand_var * high

def main (*args):
OPT_N = 4000000
iterations = 10
if len(args) >= 2:
iterations = int(args[0])

callResultNumpy = np.zeros(OPT_N)
putResultNumpy = -np.ones(OPT_N)
stockPrice = randfloat(np.random.random(OPT_N), 5.0, 30.0)
optionStrike = randfloat(np.random.random(OPT_N), 1.0, 100.0)
optionYears = randfloat(np.random.random(OPT_N), 0.25, 10.0)
callResultNumba = np.zeros(OPT_N)
putResultNumba = -np.ones(OPT_N)
callResultNumbapro = np.zeros(OPT_N)
putResultNumbapro = -np.ones(OPT_N)

time0 = time.time()
for i in range(iterations):
black_scholes(callResultNumpy, putResultNumpy, stockPrice,
optionStrike, optionYears, RISKFREE, VOLATILITY)
time1 = time.time()
print("Numpy Time: %f msec" %
((1000 * (time1 - time0)) / iterations))

time0 = time.time()
for i in range(iterations):
black_scholes_numba(callResultNumba, putResultNumba, stockPrice,
optionStrike, optionYears, RISKFREE, VOLATILITY)
time1 = time.time()
print("Numba Time: %f msec" %
((1000 * (time1 - time0)) / iterations))

time0 = time.time()
blockdim = 1024, 1
griddim = int(math.ceil(float(OPT_N)/blockdim[0])), 1
stream = cuda.stream()
d_callResult = cuda.to_device(callResultNumbapro, stream)
d_putResult = cuda.to_device(putResultNumbapro, stream)
d_stockPrice = cuda.to_device(stockPrice, stream)
d_optionStrike = cuda.to_device(optionStrike, stream)
d_optionYears = cuda.to_device(optionYears, stream)
time1 = time.time()
for i in range(iterations):
black_scholes_cuda[griddim, blockdim, stream](
d_callResult, d_putResult, d_stockPrice, d_optionStrike,
d_optionYears, RISKFREE, VOLATILITY)
d_callResult.to_host(stream)
d_putResult.to_host(stream)
stream.synchronize()
time2 = time.time()
dt = (time1 - time0) * 10 + (time2 - time1)
print("numbapro.cuda time: %f msec" % ((1000 * dt) / iterations))

delta = np.abs(callResultNumpy - callResultNumba)
L1norm = delta.sum() / np.abs(callResultNumpy).sum()
print("L1 norm: %E" % L1norm)
print("Max absolute error: %E" % delta.max())

delta = np.abs(callResultNumpy - callResultNumbapro)
L1norm = delta.sum() / np.abs(callResultNumpy).sum()
print("L1 norm (Numbapro): %E" % L1norm)
print("Max absolute error (Numbapro): %E" % delta.max())

if __name__ == "__main__":
import sys
main(*sys.argv[1:])
　　運行結果是：

Numpy Time: 1178.500009 msec
Numba Time: 424.500012 msec
numbapro.cuda time: 138.099957 msec
　　能夠看出，該程序是經過引入cuda對象並使用即時編譯方法進行加速。 　　比較發現運行時間發現，運用numbapro方式加速效果明顯。

　　總結：

1.能夠經過GPU編寫數據並行處理的程序來加快處理速度。

2.可使用CUDA庫，例如cuRAND, cuBLAS, cuFFT 3.Python CUDA程序也能夠最大化的使用硬件資源。

其餘連接：https://dask.org/