general-purpose high-level scripting language
high-level technical computing language
The set of all points \(c \in \mathbb{C}\) for which the sequence \[ z_0:=c\] \[z_{n+1} = z_n^2 +c, \quad n \in \mathbb{N}\] is bounded.
def mandelbox (cent, xsize, ysize):
si = 2000
[X,Y] = np.mgrid[cent[0]-xsize/2.:cent[0]+xsize/2.:si*1j,
cent[1]-ysize/2.:cent[1]+ysize/2.:si*1j]
C = (X + 1j*Y)
Z = np.copy(C)
it_max = 500
Anz = np.zeros(Z.shape)
for k in range(1,it_max+1):
Z = Z**2+C
Anz += np.isfinite(Z)
function mandelbox(cent, xsize, ysize)
si = 2000
x = linspace(cent[1]-xsize/2., cent[1]+xsize/2., si)
y = linspace(cent[2]-ysize/2., cent[2]+ysize/2., si)
X, Y = [i for i in x, j in y], [j for i in x, j in y]
C = (X + Y*im)
Z = copy(C)
it_max = 500
Anz = zeros(si,si)
for k in range(1,it_max)
Z = Z.^2 + C
Anz += isfinite(Z)
end
end
for i in range(0,si):
for j in range(0,si):
Z[i,j] = Z[i,j]**2 + C[i,j]
@inbounds for i in eachindex(C)
Z[i] = Z[i]^2 + C[i]
end
def mandel(z):
c = z
maxiter = 50
for n in range(0,maxiter):
if abs(z) > 2:
return n-1
z = z**2 + c
return maxiter
for a in range(0,si):
for b in range(0,si):
Anz[a,b] = mandel(C[a,b])
@everywhere function mandel(z)
c = z
maxiter = 50
for n = range(1,maxiter)
if abs(z) > 2
return n-1
end
z = z^2 + c
end
return maxiter
end
for a in eachindex(C)
Anz[a] = mandel(C[a])
end
def mandelonep(C, imin, imax):
si = np.shape(C)
res = np.zeros((imax-imin,))
for idx in range(imin,imax):
res[idx-imin] = mandel(C[np.unravel_index(idx,si)])
p = Pool(2)
slices = distribute(int(si*si),2)
results = [p.apply_async(mandelonep, (C, imin, imax))
for (imin, imax) in slices]
for r, (imin, imax) in zip(results, slices):
Anz[np.unravel_index(range(imin,imax),(si,si))] = r.get()
function mandelonep(cent, xsize, ysize, si, C,Anz)
@sync @parallel for a in eachindex(C)
Anz[a] = mandel(C[a])
end
Anz = SharedArray(Float64,(si,si))
mandelonep(cent, xsize, ysize, si, C, Anz)
type | python | julia |
---|---|---|
matrix | 50s | 93s |
loop | 86m 20s | 38s |
pixelwise | 31m 59s | 11s |
parallel | 23m 22s | 6.15s |
Note: for matrix: internally parallized through BLAS
functions have specific methods for specific argument types
function mandel(z::Complex{Float64})
function mandel(z::Array{Float64,2})
type annotation + compilation
def cython_mandel(double x,double y):
cdef np.ndarray[double,ndim=1] x = linspace(-2.1,1.2,pointsx)
cdef np.ndarray[double,ndim=1] y = linspace(-1.1,1.1,pointsy)
::Int
see types:
@code_warntype somecode
can call python transparently (without knowing much of python)
using PyCall
@pyimport mayavi.mlab as m
a = m.surf(Anz)
m.colorbar()
m.show()