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

Warning fixes, Performance improvements

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-06-13 09:45:48 +02:00
parent 0ad97c6e66
commit ae16e77de7
1 changed files with 25 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
python/_project.pyx
View File

@ -1,6 +1,7 @@
from cpython cimport bool from cpython cimport bool
from cython cimport view from cython cimport view
from cython.parallel import prange, parallel from cython.parallel import prange, parallel
from openmp cimport omp_get_max_threads, omp_get_thread_num
from libc.math cimport sin, cos, abs, floor, sqrt from libc.math cimport sin, cos, abs, floor, sqrt
import numpy as np import numpy as np
cimport numpy as npx cimport numpy as npx
@ -296,7 +297,8 @@ cdef void INTERNAL_project_cic_no_mass(npx.ndarray[DTYPE_t, ndim=3] g,
cdef double half_Box = 0.5*Lbox cdef double half_Box = 0.5*Lbox
cdef double delta_Box = Ngrid/Lbox cdef double delta_Box = Ngrid/Lbox
cdef int i cdef int i
cdef double a[3], c[3] cdef double a[3]
cdef double c[3]
cdef int b[3] cdef int b[3]
cdef int do_not_put cdef int do_not_put
@ -332,8 +334,10 @@ cdef void INTERNAL_project_cic_no_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g,
cdef double half_Box = 0.5*Lbox cdef double half_Box = 0.5*Lbox
cdef double delta_Box = Ngrid/Lbox cdef double delta_Box = Ngrid/Lbox
cdef int i cdef int i
cdef double a[3], c[3] cdef double a[3]
cdef int b[3], b1[3] cdef double c[3]
cdef int b[3]
cdef int b1[3]
cdef int do_not_put cdef int do_not_put
for i in range(x.shape[0]): for i in range(x.shape[0]):
@ -371,7 +375,8 @@ cdef void INTERNAL_project_cic_with_mass(npx.ndarray[DTYPE_t, ndim=3] g,
cdef double half_Box = 0.5*Lbox cdef double half_Box = 0.5*Lbox
cdef double delta_Box = Ngrid/Lbox cdef double delta_Box = Ngrid/Lbox
cdef int i cdef int i
cdef double a[3], c[3] cdef double a[3]
cdef double c[3]
cdef DTYPE_t m0 cdef DTYPE_t m0
cdef int b[3] cdef int b[3]
@ -408,8 +413,10 @@ cdef void INTERNAL_project_cic_with_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g
cdef double half_Box = 0.5*Lbox, m0 cdef double half_Box = 0.5*Lbox, m0
cdef double delta_Box = Ngrid/Lbox cdef double delta_Box = Ngrid/Lbox
cdef int i cdef int i
cdef double a[3], c[3] cdef double a[3]
cdef int b[3], b1[3] cdef double c[3]
cdef int b[3]
cdef int b1[3]
for i in range(x.shape[0]): for i in range(x.shape[0]):
@ -545,7 +552,8 @@ cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
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]) nogil:
cdef DTYPE_t alpha_max cdef DTYPE_t alpha_max
cdef DTYPE_t I, tmp_a cdef DTYPE_t I, tmp_a
cdef DTYPE_t v[3], term[4] cdef DTYPE_t v[3]
cdef DTYPE_t term[4]
cdef int i, j, q cdef int i, j, q
alpha_max = 10.0 # A big number alpha_max = 10.0 # A big number
@ -594,7 +602,8 @@ 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]) nogil:
cdef DTYPE_t vertex_value[8] cdef DTYPE_t vertex_value[8]
cdef DTYPE_t d cdef DTYPE_t d
cdef int a[3][2], i cdef int a[3][2]
cdef int i
for i in xrange(3): for i in xrange(3):
# if u[i] < 0: # if u[i] < 0:
@ -629,7 +638,10 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
DTYPE_t min_distance, 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) nogil except? 0:
cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3] cdef DTYPE_t u[3]
cdef DTYPE_t ifu0[3]
cdef DTYPE_t u0[3]
cdef DTYPE_t utot[3]
cdef int u_delta[3] cdef int u_delta[3]
cdef int iu0[3] cdef int iu0[3]
cdef int i cdef int i
@ -772,17 +784,17 @@ def spherical_projection(int Nside,
if progress != 0: if progress != 0:
p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start() p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start()
N = outm.size N = omp_get_max_threads()
job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int)) job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int))
job_done[:] = 0 job_done[:] = 0
theta,phi = hp.pix2ang(Nside, np.arange(N)) theta,phi = hp.pix2ang(Nside, np.arange(outm.size))
with nogil, parallel(): with nogil, parallel():
for i in prange(N): for i in prange(N):
if progress != 0 and (i%booster) == 0: if omp_get_thread_num() == 0 and progress != 0 and (i%booster) == 0:
with gil: with gil:
p.update(_mysum(job_done)) 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)
job_done[i] = 1 job_done[omp_get_thread_num()] = i
if progress: if progress:
p.finish() p.finish()