from numba import autojit
from numba.decorators import jit
import numexpr as ne
import numpy as np
import math
from numba import int32, float32, int16
import time
@jit(argtypes=[int16, float32[:,:]], restype=float32[:,:])
def ra_numba(doy, lat):
M, N = lat.shape
ra = np.zeros_like(lat)
Gsc = 0.0820
# math.pi doesnt work?
# NumbaError: 11:31: Binary operations mul on values typed object_ and object_ not (yet) supported)
pi = 3.1415926535897932384626433832795
dr = 1 + 0.033 * math.cos( 2 * pi / 365 * doy)
decl = 0.409 * math.sin( 2 * pi / 365 * doy - 1.39 )
for i in range(M):
for j in range(N):
# it crashes without the float() wrapped around the array slicing?!
ws = math.acos(-1 * math.tan(float(lat[i,j])) * math.tan(decl))
ra[i,j] = 24 * 60 / pi * Gsc * dr * ( ws * math.sin(float(lat[i,j])) * math.sin(decl) + math.cos(float(lat[i,j])) * math.cos(decl) * math.sin(ws)) * 11.6
return ra
def ra_python(doy, lat):
M, N = lat.shape
ra = np.zeros_like(lat)
Gsc = 0.0820
pi = math.pi
dr = 1 + 0.033 * math.cos( 2 * pi / 365 * doy)
decl = 0.409 * math.sin( 2 * pi / 365 * doy - 1.39 )
for i in range(M):
for j in range(N):
ws = math.acos(-1 * math.tan(lat[i,j]) * math.tan(decl))
ra[i,j] = 24 * 60 / pi * Gsc * dr * ( ws * math.sin(lat[i,j]) * math.sin(decl) + math.cos(lat[i,j]) * math.cos(decl) * math.sin(ws)) * 11.6
return ra
def ra_numpy(doy, lat):
Gsc = 0.0820
pi = math.pi
dr = 1 + 0.033 * np.cos( 2 * pi / 365 * doy)
decl = 0.409 * np.sin( 2 * pi / 365 * doy - 1.39 )
ws = np.arccos(-np.tan(lat) * np.tan(decl))
ra = 24 * 60 / pi * Gsc * dr * ( ws * np.sin(lat) * np.sin(decl) + np.cos(lat) * np.cos(decl) * np.sin(ws)) * 11.6
return ra
def ra_numxpr(doy, lat):
Gsc = 0.0820
pi = math.pi
dr = ne.evaluate("1 + 0.033 * cos( 2 * pi / 365 * doy)")
decl = ne.evaluate("0.409 * sin( 2 * pi / 365 * doy - 1.39 )")
ws = ne.evaluate("arccos(-tan(lat) * tan(decl))")
ra = ne.evaluate("24 * 60 / pi * Gsc * dr * ( ws * sin(lat) * sin(decl) + cos(lat) * cos(decl) * sin(ws)) * 11.6")
return ra
doy = 120 # day of year
py = []
nump = []
numb = []
nume = []
dims = []
for dim in [25,50,100,200,400,800,1600,3200]:
dims.append(dim)
lat = np.deg2rad(np.ones((dim,dim), dtype=np.float32) * 45.) # array of 45 degrees latitude converted to rad
tic = time.clock()
ra_nb = ra_numba(doy, lat)
numb.append(time.clock() - tic)
tic = time.clock()
ra_np = ra_numpy(doy, lat)
nump.append(time.clock() - tic)
tic = time.clock()
ra_py = ra_python(doy, lat)
py.append(time.clock() - tic)
tic = time.clock()
ra_ne = ra_numxpr(doy, lat)
nume.append(time.clock() - tic)
dims = np.array(dims)**2
py = np.array(py)
numb = np.array(numb)
nump = np.array(nump)
nume = np.array(nume)
fig, axs = plt.subplots(1,2, sharex=True)
fig.set_figwidth(20)
fig.set_figheight(5)
axs[0].plot(dims, numb, 'o-', label='Numba', ms=5, mew=0)
axs[0].plot(dims, nump, 'o-', label='Numpy', ms=5, mew=0)
axs[0].plot(dims, py, 'o-', label='Python', ms=5, mew=0)
axs[0].plot(dims, nume, 'o-', label='NumExpr', ms=5, mew=0)
axs[0].set_title('Absolute duration')
axs[0].set_xlabel('Array size [n]')
axs[0].set_ylabel('Duration [seconds]')
axs[0].legend(frameon=False, loc=2)
axs[0].set_xscale('log')
axs[0].set_yscale('log')
axs[0].grid(True)
axs[1].plot(dims, py / numb, 'o-', label='Python vs Numba', ms=5, mew=0)
axs[1].plot(dims, py / nump, 'o-', label='Python vs Numpy', ms=5, mew=0)
axs[1].plot(dims, py / nume, 'o-', label='Python vs NumExpr', ms=5, mew=0)
axs[1].set_title('Relative to python (py / xx)')
axs[1].set_xlabel('Array size [n]')
axs[1].set_ylabel('Relative duration [-]')
axs[1].legend(frameon=False, loc=2)
axs[1].set_xscale('log')
axs[1].set_ylim(0,120)
axs[1].grid(True)
#ra_nb = ra_numba(doy, lat)
#ra_np = ra_numpy(doy, lat)
#ra_py = ra_python(doy, lat)
print 'Allclose() test:'
print 'Numpy / Numba:\t\t', np.allclose(ra_np,ra_nb)
print 'Numpy / Python:\t\t', np.allclose(ra_np,ra_py)
print 'Numpy / NumExpr:\t', np.allclose(ra_np,ra_ne)
Allclose() test: Numpy / Numba: True Numpy / Python: True Numpy / NumExpr: True