diff --git a/python/_project.pyx b/python/_project.pyx index af65590..89f0c28 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -11,7 +11,7 @@ __all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE" cdef extern from "project_tool.hpp" namespace "": - DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho) + DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho) nogil @cython.boundscheck(False) @@ -561,15 +561,19 @@ 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 + for i in range(3): + u0[i] += u[i]*alpha_max + # alpha_max is the integration length # we integrate between 0 and alpha_max (curvilinear coordinates) + r[0] = j return compute_projection(vertex_value, u, u0, alpha_max) @cython.boundscheck(False) cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, - DTYPE_t u[3], DTYPE_t u0[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 d d = density[iu0[0], iu0[1], iu0[2]] @@ -578,13 +582,37 @@ 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 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 d + cdef int a[3][2], i - d = density[iu0[0], iu0[1], iu0[2]] - return cube_integral(u, u0, jumper)*d - cube_integral_trilin(u, u0, jumper, vertex_value) + for i in xrange(3): +# if u[i] < 0: +# a[i][0] = iu0[i]-1 +# a[i][1] = iu0[i] +# else: + a[i][0] = iu0[i]-1 + a[i][1] = iu0[i] + + with gil: + assert a[i][0] >= 0 + assert a[i][1] < density.shape[i] + assert u0[i] >=0 + assert u0[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]] + vertex_value[1 + 2*1 + 4*0] = density[a[0][1], a[1][1], a[2][0]] + + vertex_value[0 + 2*0 + 4*1] = density[a[0][0], a[1][0], a[2][1]] + vertex_value[1 + 2*0 + 4*1] = density[a[0][1], a[1][0], a[2][1]] + 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(u, u0, jumper)*d + return cube_integral_trilin(u, u0, jumper, vertex_value) @cython.boundscheck(False) @@ -594,6 +622,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, DTYPE_t max_distance, int integrator_id=0): cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3] + cdef int u_delta[3] cdef int iu0[3] cdef int i cdef int N = density.shape[0] @@ -603,7 +632,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, cdef int jumper[1] cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:], - DTYPE_t u[3], DTYPE_t u0[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 if integrator_id == 0: integrator = integrator0 @@ -620,6 +649,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, return 0 iu0[i] = int(floor(ifu0[i])) 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))): raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i])) if (not (u0[i]>=0 and u0[i]<=1)): @@ -635,15 +665,13 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, jumper[0] = 0 while completed == 0: - I0 += integrator(density, u, u0, iu0, jumper) + I0 += integrator(density, u, u0, u_delta, iu0, jumper) if u[jumper[0]] < 0: iu0[jumper[0]] -= 1 - direction = -1 u0[jumper[0]] = 1 else: iu0[jumper[0]] += 1 - direction = 1 u0[jumper[0]] = 0 @@ -669,7 +697,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, 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): + DTYPE_t max_distance, int progress=1, int integrator_id=0): import healpy as hp import progressbar as pb @@ -686,7 +714,7 @@ def spherical_projection(int Nside, if progress: p.update(i) u = np.array(hp.pix2vec(Nside, i), dtype=DTYPE) - outm[i] = line_of_sight_projection(density, u, min_distance, max_distance) + outm[i] = line_of_sight_projection(density, u, min_distance, max_distance, integrator_id=integrator_id) if progress: p.finish() diff --git a/python/project_tool.hpp b/python/project_tool.hpp index 9322b5f..c13f3f2 100644 --- a/python/project_tool.hpp +++ b/python/project_tool.hpp @@ -12,7 +12,7 @@ static T project_tool(T *vertex_value, T *u, T *u0) T ret = 0; for (int q = 0; q < 3; q++) - epsilon[q] = 2*c[q] - 1; + epsilon[q] = -(2*c[q]-1); for (int q = 0; q < ProdType::numProducts; q++) ret += ProdType::product(u, u0, epsilon, q); @@ -24,6 +24,12 @@ static T project_tool(T *vertex_value, T *u, T *u0) } +template +static T get_u0(T u0, int epsilon) +{ + return (epsilon < 0) ? u0 : (1-u0); +} + template struct ProductTerm0 { @@ -34,11 +40,12 @@ struct ProductTerm0 T a = 1; for (int r = 0; r < 3; r++) - a *= (epsilon[r] < 0) ? u0[r] : (1-u0[r]); + a *= get_u0(u0[r], epsilon[r]); return a; } }; + template struct ProductTerm1 { @@ -51,12 +58,12 @@ struct ProductTerm1 for (int r = 0; r < 3; r++) { - G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]); + G[r] = get_u0(u0[r], epsilon[r]); } - double F[3] = { G[0]*u[1]*u[2], u[0]*G[1]*u[2], u[0]*u[1]*G[2] }; + double F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] }; - return F[q] * epsilon[q]; + return F[q] * u[q] * epsilon[q]; } }; @@ -72,16 +79,17 @@ struct ProductTerm2 for (int r = 0; r < 3; r++) { - G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]); + G[r] = get_u0(u0[r], epsilon[r]); } - double F[3] = { u[0]*G[1]*G[2], G[0]*u[1]*G[2], G[0]*G[1]*u[2] }; + double F[3] = { epsilon[1]*epsilon[2]*u[1]*u[2], epsilon[0]*epsilon[2]*u[0]*u[2], epsilon[0]*epsilon[1]*u[0]*u[1] }; - return F[q] * epsilon[q]; + return F[q] * G[q]; } }; + template struct ProductTerm3 { @@ -89,9 +97,7 @@ struct ProductTerm3 static T product(T *u, T *u0, int *epsilon, int q) { - T a = 1; - - return epsilon[0]*epsilon[1]*epsilon[2]; + return epsilon[0]*epsilon[1]*epsilon[2]*u[0]*u[1]*u[2]; } }; @@ -103,8 +109,8 @@ T compute_projection(T *vertex_value, T *u, T *u0, T rho) ret = project_tool >(vertex_value, u, u0) * rho; ret += project_tool >(vertex_value, u, u0) * rho * rho / 2; - ret += project_tool >(vertex_value, u, u0) * rho * rho * rho / 3; - ret += project_tool >(vertex_value, u, u0) * rho * rho * rho * rho / 4; +// ret += project_tool >(vertex_value, u, u0) * rho * rho * rho / 3; +// ret += project_tool >(vertex_value, u, u0) * rho * rho * rho * rho / 4; return ret; } diff --git a/python_sample/test_spheric_proj.py b/python_sample/test_spheric_proj.py new file mode 100644 index 0000000..312a387 --- /dev/null +++ b/python_sample/test_spheric_proj.py @@ -0,0 +1,10 @@ +import cosmotool as ct +import numpy as np +import healpy as hp + +d = np.zeros((64,64,64), ct.DTYPE) + +d[33,33,33] = 1 + +proj0 = ct.spherical_projection(32, d, 0, 20, integrator_id=0) +proj1 = ct.spherical_projection(32, d, 0, 20, integrator_id=1)