diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index c646f25..5bc6783 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -11,7 +11,7 @@ add_custom_command( add_custom_command( OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx) + DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx ${CMAKE_CURRENT_SOURCE_DIR}/project_tool.hpp ) ENDIF(CYTHON) diff --git a/python/_project.pyx b/python/_project.pyx index 89f0c28..c3d39a0 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -1,4 +1,5 @@ from cpython cimport bool +from cython cimport view from libc.math cimport sin, cos, abs, floor, sqrt import numpy as np cimport numpy as npx @@ -6,6 +7,7 @@ cimport cython ctypedef npx.float64_t DTYPE_t DTYPE=np.float64 +FORMAT_DTYPE="d" __all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"] @@ -541,7 +543,7 @@ cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil: @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 tmp_a + cdef DTYPE_t I, tmp_a cdef DTYPE_t v[3], term[4] cdef int i, j, q @@ -561,15 +563,21 @@ 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): + + with gil: + assert u0[i] >=0 + assert u0[i] <= 1 + + I = compute_projection(vertex_value, u, u0, alpha_max) + + for i in xrange(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) - + return I @cython.boundscheck(False) cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, @@ -592,14 +600,12 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, # a[i][0] = iu0[i]-1 # a[i][1] = iu0[i] # else: - a[i][0] = iu0[i]-1 - a[i][1] = iu0[i] + a[i][0] = iu0[i] + a[i][1] = iu0[i]+1 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]] @@ -619,7 +625,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density, def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, npx.ndarray[DTYPE_t] a_u, DTYPE_t min_distance, - DTYPE_t max_distance, int integrator_id=0): + DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id=0): cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3] cdef int u_delta[3] @@ -644,7 +650,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density, for i in range(3): u[i] = a_u[i] u0[i] = a_u[i]*min_distance - ifu0[i] = half_N+u0[i] + ifu0[i] = half_N+u0[i]+shifter[i] if (ifu0[i] <= 0 or ifu0[i] >= N): return 0 iu0[i] = int(floor(ifu0[i])) @@ -697,7 +703,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, int integrator_id=0): + DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None): import healpy as hp import progressbar as pb @@ -705,6 +711,10 @@ def spherical_projection(int Nside, cdef npx.ndarray[DTYPE_t, ndim=1] u cdef npx.ndarray[DTYPE_t, ndim=1] outm + if shifter is None: + shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t)) + shifter[:] = 0 + outm = np.empty(hp.nside2npix(Nside),dtype=DTYPE) if progress: @@ -714,7 +724,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, integrator_id=integrator_id) + outm[i] = line_of_sight_projection(density, u, min_distance, max_distance, shifter, integrator_id=integrator_id) if progress: p.finish() diff --git a/python/project_tool.hpp b/python/project_tool.hpp index c13f3f2..ad36c20 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); @@ -25,9 +25,9 @@ static T project_tool(T *vertex_value, T *u, T *u0) template -static T get_u0(T u0, int epsilon) +static T get_u0(const T& u0, int epsilon) { - return (epsilon < 0) ? u0 : (1-u0); + return (epsilon > 0) ? u0 : (1-u0); } template @@ -61,7 +61,7 @@ struct ProductTerm1 G[r] = get_u0(u0[r], epsilon[r]); } - double F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] }; + double F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] }; return F[q] * u[q] * epsilon[q]; } @@ -108,7 +108,7 @@ T compute_projection(T *vertex_value, T *u, T *u0, T rho) T ret; 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 / 2; // 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;