From 21507e301393b2b09e850c2e65914832d1c0572b Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 14 Jul 2015 12:21:20 +0200 Subject: [PATCH] Fixed integrator to stop at exactly the right distance --- python/_project.pyx | 46 ++++++++++++++++++++++----------------------- src/mykdtree.tcc | 2 +- 2 files changed, 23 insertions(+), 25 deletions(-) diff --git a/python/_project.pyx b/python/_project.pyx index e432495..47c7091 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -600,14 +600,11 @@ def tophat_fourier(x not None): @cython.boundscheck(False) @cython.cdivision(True) -cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) nogil: - cdef DTYPE_t alpha_max +cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t alpha_max) nogil: cdef DTYPE_t tmp_a cdef DTYPE_t v[3] cdef int i, j - alpha_max = 10.0 # A big number - for i in xrange(3): if u[i] == 0.: continue @@ -630,14 +627,11 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) 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 alpha_max +cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], DTYPE_t alpha_max) nogil: cdef DTYPE_t I, tmp_a cdef DTYPE_t v[3], term[4] cdef int i, j, q - alpha_max = 10.0 # A big number - j = 0 for i in range(3): if u[i] == 0.: @@ -651,7 +645,7 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t if tmp_a < alpha_max: alpha_max = tmp_a j = i - + I = compute_projection(vertex_value, u, u0, alpha_max) for i in xrange(3): @@ -665,16 +659,16 @@ 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], DTYPE_t alpha_max) nogil: cdef DTYPE_t d d = density[iu0[0], iu0[1], iu0[2]] - return cube_integral(u, u0, jumper)*d + return cube_integral(u, u0, jumper, alpha_max)*d @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], DTYPE_t alpha_max) nogil: cdef DTYPE_t vertex_value[8] cdef DTYPE_t d cdef int a[3][2], i @@ -682,7 +676,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, for i in xrange(3): a[i][0] = iu0[i] a[i][1] = iu0[i]+1 - + vertex_value[0 + 2*0 + 4*0] = density[a[0][0], a[1][0], a[2][0]] vertex_value[1 + 2*0 + 4*0] = density[a[0][1], a[1][0], a[2][0]] vertex_value[0 + 2*1 + 4*0] = density[a[0][0], a[1][1], a[2][0]] @@ -693,7 +687,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, vertex_value[0 + 2*1 + 4*1] = density[a[0][0], a[1][1], a[2][1]] vertex_value[1 + 2*1 + 4*1] = density[a[0][1], a[1][1], a[2][1]] - return cube_integral_trilin(u, u0, jumper, vertex_value) + return cube_integral_trilin(u, u0, jumper, vertex_value, alpha_max) @@ -714,7 +708,7 @@ 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], DTYPE_t alpha_max) nogil if integrator_id == 0: integrator = integrator0 @@ -741,16 +735,17 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, raise RuntimeError("u0[%d] = %g !" % (i,u0[i])) completed = 0 - if ((iu0[0] >= N) or (iu0[0] <= 0) or - (iu0[1] >= N) or (iu0[1] <= 0) or - (iu0[2] >= N) or (iu0[2] <= 0)): + if ((iu0[0] >= N-1) or (iu0[0] <= 0) or + (iu0[1] >= N-1) or (iu0[1] <= 0) or + (iu0[2] >= N-1) or (iu0[2] <= 0)): completed = 1 I0 = 0 jumper[0] = 0 + dist2 = 0 while completed == 0: - I0 += integrator(density, u, u0, u_delta, iu0, jumper) + I0 += integrator(density, u, u0, u_delta, iu0, jumper, max_distance-sqrt(dist2)) if u[jumper[0]] < 0: iu0[jumper[0]] -= 1 @@ -760,9 +755,9 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, u0[jumper[0]] = 0 - if ((iu0[0] >= N) or (iu0[0] <= 0) or - (iu0[1] >= N) or (iu0[1] <= 0) or - (iu0[2] >= N) or (iu0[2] <= 0)): + if ((iu0[0] >= N-1) or (iu0[0] <= 0) or + (iu0[1] >= N-1) or (iu0[1] <= 0) or + (iu0[2] >= N-1) or (iu0[2] <= 0)): completed = 1 else: dist2 = 0 @@ -788,7 +783,7 @@ def line_of_sight_projection(DTYPE_t[:,:,:] density not None, u[1] = a_u[1] u[2] = a_u[2] - C_line_of_sight_projection(density, + return C_line_of_sight_projection(density, u, min_distance, max_distance, shifter, integrator_id) @@ -846,10 +841,13 @@ def spherical_projection(int Nside, shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t)) shifter[:] = 0 + print("allocating map") outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE) + print("initializing views") outm = outm_array density_view = density + print("progress?") if progress != 0: p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start() @@ -864,7 +862,7 @@ def spherical_projection(int Nside, theta,phi = hp.pix2ang(Nside, np.arange(N0)) with nogil, parallel(): tid = smp_get_thread_id() - for i in prange(N0): + for i in prange(N0,schedule='dynamic',chunksize=256): if progress != 0 and (i%booster) == 0: with gil: p.update(_mysum(job_done)) diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc index e9ce02c..e93575d 100644 --- a/src/mykdtree.tcc +++ b/src/mykdtree.tcc @@ -280,7 +280,7 @@ namespace CosmoTool { coords tmpBound; NodeIntType nodeId; -#pragma omp atomic +#pragma omp atomic capture nodeId = (this->lastNode)++; node = &nodes[nodeId];