diff --git a/python/_project.pyx b/python/_project.pyx index 8489cc9..f2c31d9 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -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,8 +652,10 @@ 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 else: @@ -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: - - I0 += integrator(density, u, u0, u_delta, iu0, jumper) + err[0] = 0 + I0 += integrator(density, u, u0, u_delta, iu0, jumper, err) + + 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: