guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

attempt to boost python performance

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-06-13 11:07:58 +02:00
parent ae16e77de7
commit 7a61ad7519
1 changed files with 38 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
python/_project.pyx
View File

@ -549,7 +549,7 @@ cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
@cython.boundscheck(False)
@cython.cdivision(True)
cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]) nogil:
cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], int err[1]) nogil:
cdef DTYPE_t alpha_max
cdef DTYPE_t I, tmp_a
cdef DTYPE_t v[3]
@ -572,10 +572,9 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
alpha_max = tmp_a
j = i
with gil:
assert u0[i] >=0
assert u0[i] <= 1
if (u0[i] >= 0 || u0[i] <= 1)
err[0] = 1
return 0
I = compute_projection(vertex_value, u, u0, alpha_max)
@ -590,7 +589,7 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
@cython.boundscheck(False)
cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil:
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil:
cdef DTYPE_t d
d = density[iu0[0], iu0[1], iu0[2]]
@ -599,7 +598,7 @@ cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
@cython.boundscheck(False)
cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil:
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil:
cdef DTYPE_t vertex_value[8]
cdef DTYPE_t d
cdef int a[3][2]
@ -628,7 +627,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
vertex_value[1 + 2*1 + 4*1] = density[a[0][1], a[1][1], a[2][1]]
# return cube_integral(u, u0, jumper)*d
return cube_integral_trilin(u, u0, jumper, vertex_value)
return cube_integral_trilin(u, u0, jumper, vertex_value, err)
@ -636,7 +635,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
DTYPE_t a_u[3],
DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 0:
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id, int out_err[1]) nogil:
cdef DTYPE_t u[3]
cdef DTYPE_t ifu0[3]
@ -644,6 +643,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
cdef DTYPE_t utot[3]
cdef int u_delta[3]
cdef int iu0[3]
cdef int err[1]
cdef int i
cdef int N = density.shape[0]
cdef int half_N = density.shape[0]/2
@ -652,7 +652,9 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
cdef int jumper[1]
cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], int err[1]) nogil
out_err[0] = 0
if integrator_id == 0:
integrator = integrator0
@ -671,12 +673,12 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
u0[i] = ifu0[i]-iu0[i]
u_delta[i] = 1 if iu0[i] > 0 else -1
if (not ((iu0[i]>= 0) and (iu0[i] < N))):
with gil:
raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i]))
out_err[0] = 1
return 0
if (not (u0[i]>=0 and u0[i]<=1)):
with gil:
raise RuntimeError("u0[%d] = %g !" % (i,u0[i]))
out_err[0] = 1
return 0
completed = 0
if ((iu0[0] >= N) or (iu0[0] <= 0) or
@ -687,8 +689,12 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
I0 = 0
jumper[0] = 0
while completed == 0:
err[0] = 0
I0 += integrator(density, u, u0, u_delta, iu0, jumper, err)
I0 += integrator(density, u, u0, u_delta, iu0, jumper)
if err[0] == 1:
out_err[0] = 1
return 0
if u[jumper[0]] < 0:
iu0[jumper[0]] -= 1
@ -721,19 +727,24 @@ def line_of_sight_projection(DTYPE_t[:,:,:] density not None,
DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter not None, int integrator_id=0):
cdef DTYPE_t u[3]
cdef int out_err[1]
cdef DTYPE_t v
u[0] = a_u[0]
u[1] = a_u[1]
u[2] = a_u[2]
C_line_of_sight_projection(density,
v = C_line_of_sight_projection(density,
u,
min_distance,
max_distance, shifter, integrator_id)
max_distance, shifter, integrator_id, out_err)
if out_err[0] == 1:
raise RuntimeError("Error occured during integration")
return v
cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density,
double min_distance, double max_distance,
DTYPE_t[:] shifter, int integrator_id) nogil:
DTYPE_t[:] shifter, int integrator_id, int out_err[1]) nogil:
cdef DTYPE_t u0[3]
stheta = sin(theta)
@ -741,7 +752,7 @@ cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density
u0[1] = sin(phi)*stheta
u0[2] = cos(theta)
return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id)
return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id, out_err)
@cython.boundscheck(False)
cdef npx.uint64_t _mysum(int[:] jobs) nogil:
@ -760,7 +771,7 @@ cdef npx.uint64_t _mysum(int[:] jobs) nogil:
def spherical_projection(int Nside,
npx.ndarray[DTYPE_t, ndim=3] density not None,
DTYPE_t min_distance,
DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=1000):
DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=-1):
import healpy as hp
import progressbar as pb
@ -772,6 +783,7 @@ def spherical_projection(int Nside,
cdef npx.ndarray[DTYPE_t, ndim=1] outm_array
cdef long N
cdef double stheta
cdef int out_err[1]
if shifter is None:
shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
@ -788,12 +800,16 @@ def spherical_projection(int Nside,
job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int))
job_done[:] = 0
theta,phi = hp.pix2ang(Nside, np.arange(outm.size))
if booster <= 0:
booster = density.size / 100 / N
with nogil, parallel():
for i in prange(N):
if omp_get_thread_num() == 0 and progress != 0 and (i%booster) == 0:
with gil:
p.update(_mysum(job_done))
outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id)
outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id, out_err)
job_done[omp_get_thread_num()] = i
if progress: