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/
https://anaconda.org/talley/repo?type=conda&label=main
---------------------
作者：WeisongZhao
來源：CSDN
原文：https://blog.csdn.net/weixin_41923961/article/details/83687809
版權聲明：本文為博主原創文章，轉載請附上博文鏈接！