From 1693c222e113c4b9827c836519a0a27ffbbf0c2f Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 4 Jun 2014 14:00:54 +0200 Subject: [PATCH] Imported the projection library for python --- python/CMakeLists.txt | 9 +- python/_project.pyx | 664 +++++++++++++++++++++++++++++++++++ python/cosmotool/__init__.py | 1 + python/cosmotool/grafic.py | 5 +- 4 files changed, 676 insertions(+), 3 deletions(-) create mode 100644 python/_project.pyx diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index e7d4653..cb0f7eb 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -8,14 +8,21 @@ add_custom_command( COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx) +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) + ENDIF(CYTHON) add_library(_cosmotool MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp) +add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp) SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Bsymbolic-functions") target_link_libraries(_cosmotool ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LIBRARIES}) +target_link_libraries(_project ${PYTHON_LIBRARIES}) # Discover where to put packages if (NOT PYTHON_SITE_PACKAGES) @@ -41,7 +48,7 @@ if (WIN32 AND NOT CYGWIN) SET_TARGET_PROPERTIES(_cosmotool PROPERTIES SUFFIX ".pyd") endif (WIN32 AND NOT CYGWIN) -INSTALL(TARGETS _cosmotool +INSTALL(TARGETS _cosmotool _project LIBRARY DESTINATION ${PYTHON_SITE_PACKAGES}/cosmotool ) diff --git a/python/_project.pyx b/python/_project.pyx new file mode 100644 index 0000000..fe4f3a9 --- /dev/null +++ b/python/_project.pyx @@ -0,0 +1,664 @@ +from cpython cimport bool +from libc.math cimport sin, cos, abs, floor, sqrt +import numpy as np +cimport numpy as npx +cimport cython + +ctypedef npx.float64_t DTYPE_t +DTYPE=np.float64 + +__all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"] + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef DTYPE_t interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, + DTYPE_t z, + npx.ndarray[DTYPE_t, ndim=3] d, DTYPE_t Lbox) except? 0: + + cdef int Ngrid = d.shape[0] + cdef DTYPE_t inv_delta = Ngrid/Lbox + cdef int ix, iy, iz + cdef DTYPE_t f[2][2][2] + cdef DTYPE_t rx, ry, rz + cdef int jx, jy, jz + + rx = (inv_delta*x + Ngrid/2) + ry = (inv_delta*y + Ngrid/2) + rz = (inv_delta*z + Ngrid/2) + + ix = int(floor(rx)) + iy = int(floor(ry)) + iz = int(floor(rz)) + + + rx -= ix + ry -= iy + rz -= iz + + while ix < 0: + ix += Ngrid + while iy < 0: + iy += Ngrid + while iz < 0: + iz += Ngrid + + jx = (ix+1)%Ngrid + jy = (iy+1)%Ngrid + jz = (iz+1)%Ngrid + + ix = ix%Ngrid + iy = iy%Ngrid + iz = iz%Ngrid + + assert ((ix >= 0) and ((jx) < Ngrid)) + assert ((iy >= 0) and ((jy) < Ngrid)) + assert ((iz >= 0) and ((jz) < Ngrid)) + + f[0][0][0] = (1-rx)*(1-ry)*(1-rz) + f[1][0][0] = ( rx)*(1-ry)*(1-rz) + f[0][1][0] = (1-rx)*( ry)*(1-rz) + f[1][1][0] = ( rx)*( ry)*(1-rz) + + f[0][0][1] = (1-rx)*(1-ry)*( rz) + f[1][0][1] = ( rx)*(1-ry)*( rz) + f[0][1][1] = (1-rx)*( ry)*( rz) + f[1][1][1] = ( rx)*( ry)*( rz) + + return \ + d[ix ,iy ,iz ] * f[0][0][0] + \ + d[jx ,iy ,iz ] * f[1][0][0] + \ + d[ix ,jy ,iz ] * f[0][1][0] + \ + d[jx ,jy ,iz ] * f[1][1][0] + \ + d[ix ,iy ,jz ] * f[0][0][1] + \ + d[jx ,iy ,jz ] * f[1][0][1] + \ + d[ix ,jy ,jz ] * f[0][1][1] + \ + d[jx ,jy ,jz ] * f[1][1][1] + + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef DTYPE_t interp3d_INTERNAL(DTYPE_t x, DTYPE_t y, + DTYPE_t z, + npx.ndarray[DTYPE_t, ndim=3] d, DTYPE_t Lbox) except? 0: + + cdef int Ngrid = d.shape[0] + cdef DTYPE_t inv_delta = Ngrid/Lbox + cdef int ix, iy, iz + cdef DTYPE_t f[2][2][2] + cdef DTYPE_t rx, ry, rz + + rx = (inv_delta*x + Ngrid/2) + ry = (inv_delta*y + Ngrid/2) + rz = (inv_delta*z + Ngrid/2) + + ix = int(floor(rx)) + iy = int(floor(ry)) + iz = int(floor(rz)) + + rx -= ix + ry -= iy + rz -= iz + + if ((ix < 0) or (ix+1) >= Ngrid): + raise IndexError("X coord out of bound (ix=%d, x=%g)" % (ix,x)) + if ((iy < 0) or (iy+1) >= Ngrid): + raise IndexError("Y coord out of bound (iy=%d, y=%g)" % (iy,y)) + if ((iz < 0) or (iz+1) >= Ngrid): + raise IndexError("Z coord out of bound (iz=%d, z=%g)" % (iz,z)) + # assert ((ix >= 0) and ((ix+1) < Ngrid)) +# assert ((iy >= 0) and ((iy+1) < Ngrid)) +# assert ((iz >= 0) and ((iz+1) < Ngrid)) + + f[0][0][0] = (1-rx)*(1-ry)*(1-rz) + f[1][0][0] = ( rx)*(1-ry)*(1-rz) + f[0][1][0] = (1-rx)*( ry)*(1-rz) + f[1][1][0] = ( rx)*( ry)*(1-rz) + + f[0][0][1] = (1-rx)*(1-ry)*( rz) + f[1][0][1] = ( rx)*(1-ry)*( rz) + f[0][1][1] = (1-rx)*( ry)*( rz) + f[1][1][1] = ( rx)*( ry)*( rz) + + return \ + d[ix ,iy ,iz ] * f[0][0][0] + \ + d[ix+1,iy ,iz ] * f[1][0][0] + \ + d[ix ,iy+1,iz ] * f[0][1][0] + \ + d[ix+1,iy+1,iz ] * f[1][1][0] + \ + d[ix ,iy ,iz+1] * f[0][0][1] + \ + d[ix+1,iy ,iz+1] * f[1][0][1] + \ + d[ix ,iy+1,iz+1] * f[0][1][1] + \ + d[ix+1,iy+1,iz+1] * f[1][1][1] + +def interp3d(x not None, y not None, + z not None, + npx.ndarray[DTYPE_t, ndim=3] d not None, DTYPE_t Lbox, + bool periodic=False): + cdef npx.ndarray[DTYPE_t] out + cdef npx.ndarray[DTYPE_t] ax, ay, az + cdef int i + + if d.shape[0] != d.shape[1] or d.shape[0] != d.shape[2]: + raise ValueError("Grid must have a cubic shape") + + if type(x) == np.ndarray or type(y) == np.ndarray or type(z) == np.ndarray: + if type(x) != np.ndarray or type(y) != np.ndarray or type(z) != np.ndarray: + raise ValueError("All or no array. No partial arguments") + + ax = x + ay = y + az = z + assert ax.size == ay.size and ax.size == az.size + + out = np.empty(x.shape, dtype=DTYPE) + if periodic: + for i in range(ax.size): + out[i] = interp3d_INTERNAL_periodic(ax[i], ay[i], az[i], d, Lbox) + else: + for i in range(ax.size): + out[i] = interp3d_INTERNAL(ax[i], ay[i], az[i], d, Lbox) + + return out + else: + if periodic: + return interp3d_INTERNAL_periodic(x, y, z, d, Lbox) + else: + return interp3d_INTERNAL(x, y, z, d, Lbox) + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef DTYPE_t interp2d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, + npx.ndarray[DTYPE_t, ndim=2] d, DTYPE_t Lbox) except? 0: + + cdef int Ngrid = d.shape[0] + cdef DTYPE_t inv_delta = Ngrid/Lbox + cdef int ix, iy + cdef DTYPE_t f[2][2] + cdef DTYPE_t rx, ry + cdef int jx, jy + + rx = (inv_delta*x + Ngrid/2) + ry = (inv_delta*y + Ngrid/2) + + ix = int(floor(rx)) + iy = int(floor(ry)) + + rx -= ix + ry -= iy + + while ix < 0: + ix += Ngrid + while iy < 0: + iy += Ngrid + + jx = (ix+1)%Ngrid + jy = (iy+1)%Ngrid + + assert ((ix >= 0) and ((jx) < Ngrid)) + assert ((iy >= 0) and ((jy) < Ngrid)) + + f[0][0] = (1-rx)*(1-ry) + f[1][0] = ( rx)*(1-ry) + f[0][1] = (1-rx)*( ry) + f[1][1] = ( rx)*( ry) + + return \ + d[ix ,iy ] * f[0][0] + \ + d[jx ,iy ] * f[1][0] + \ + d[ix ,jy ] * f[0][1] + \ + d[jx ,jy ] * f[1][1] + + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef DTYPE_t interp2d_INTERNAL(DTYPE_t x, DTYPE_t y, + npx.ndarray[DTYPE_t, ndim=2] d, DTYPE_t Lbox) except? 0: + + cdef int Ngrid = d.shape[0] + cdef DTYPE_t inv_delta = Ngrid/Lbox + cdef int ix, iy + cdef DTYPE_t f[2][2] + cdef DTYPE_t rx, ry + + rx = (inv_delta*x + Ngrid/2) + ry = (inv_delta*y + Ngrid/2) + + ix = int(floor(rx)) + iy = int(floor(ry)) + + rx -= ix + ry -= iy + + if ((ix < 0) or (ix+1) >= Ngrid): + raise IndexError("X coord out of bound (ix=%d, x=%g)" % (ix,x)) + if ((iy < 0) or (iy+1) >= Ngrid): + raise IndexError("Y coord out of bound (iy=%d, y=%g)" % (iy,y)) +# assert ((ix >= 0) and ((ix+1) < Ngrid)) +# assert ((iy >= 0) and ((iy+1) < Ngrid)) +# assert ((iz >= 0) and ((iz+1) < Ngrid)) + + f[0][0] = (1-rx)*(1-ry) + f[1][0] = ( rx)*(1-ry) + f[0][1] = (1-rx)*( ry) + f[1][1] = ( rx)*( ry) + + return \ + d[ix ,iy ] * f[0][0] + \ + d[ix+1,iy ] * f[1][0] + \ + d[ix ,iy+1] * f[0][1] + \ + d[ix+1,iy+1] * f[1][1] + +def interp2d(x not None, y not None, + npx.ndarray[DTYPE_t, ndim=2] d not None, DTYPE_t Lbox, + bool periodic=False): + cdef npx.ndarray[DTYPE_t] out + cdef npx.ndarray[DTYPE_t] ax, ay + cdef int i + + if d.shape[0] != d.shape[1]: + raise ValueError("Grid must have a square shape") + + if type(x) == np.ndarray or type(y) == np.ndarray: + if type(x) != np.ndarray or type(y) != np.ndarray: + raise ValueError("All or no array. No partial arguments") + + ax = x + ay = y + assert ax.size == ay.size + + out = np.empty(x.shape, dtype=DTYPE) + if periodic: + for i in range(ax.size): + out[i] = interp2d_INTERNAL_periodic(ax[i], ay[i], d, Lbox) + else: + for i in range(ax.size): + out[i] = interp2d_INTERNAL(ax[i], ay[i], d, Lbox) + + return out + else: + if periodic: + return interp2d_INTERNAL_periodic(x, y, d, Lbox) + else: + return interp2d_INTERNAL(x, y, d, Lbox) + + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef void INTERNAL_project_cic_no_mass(npx.ndarray[DTYPE_t, ndim=3] g, + npx.ndarray[DTYPE_t, ndim=2] x, int Ngrid, float Lbox): + cdef double half_Box = 0.5*Lbox + cdef double delta_Box = Ngrid/Lbox + cdef int i + cdef double a[3], c[3] + cdef int b[3] + cdef int do_not_put + + for i in range(x.shape[0]): + + do_not_put = False + for j in range(3): + a[j] = (x[i,j]+half_Box)*delta_Box + b[j] = int(floor(a[j])) + a[j] -= b[j] + c[j] = 1-a[j] + if b[j] < 0 or b[j]+1 >= Ngrid: + do_not_put = True + + print("b = %g %g %g" % (b[0], b[1], b[2])) + print("a = %g %g %g" % (a[0], a[1], a[2])) + + if not do_not_put: + g[b[0],b[1],b[2]] += c[0]*c[1]*c[2] + g[b[0]+1,b[1],b[2]] += a[0]*c[1]*c[2] + g[b[0],b[1]+1,b[2]] += c[0]*a[1]*c[2] + g[b[0]+1,b[1]+1,b[2]] += a[0]*a[1]*c[2] + + g[b[0],b[1],b[2]+1] += c[0]*c[1]*a[2] + g[b[0]+1,b[1],b[2]+1] += a[0]*c[1]*a[2] + g[b[0],b[1]+1,b[2]+1] += c[0]*a[1]*a[2] + g[b[0]+1,b[1]+1,b[2]+1] += a[0]*a[1]*a[2] + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef void INTERNAL_project_cic_no_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g, + npx.ndarray[DTYPE_t, ndim=2] x, int Ngrid, double Lbox): + cdef double half_Box = 0.5*Lbox + cdef double delta_Box = Ngrid/Lbox + cdef int i + cdef double a[3], c[3] + cdef int b[3], b1[3] + cdef int do_not_put + + for i in range(x.shape[0]): + + do_not_put = False + for j in range(3): + a[j] = (x[i,j]+half_Box)*delta_Box + b[j] = int(floor(a[j])) + b1[j] = b[j]+1 + while b1[j] < 0: + b1[j] += Ngrid + while b1[j] >= Ngrid: + b1[j] -= Ngrid + + a[j] -= b[j] + c[j] = 1-a[j] + + g[b[0],b[1],b[2]] += c[0]*c[1]*c[2] + g[b1[0],b[1],b[2]] += a[0]*c[1]*c[2] + g[b[0],b1[1],b[2]] += c[0]*a[1]*c[2] + g[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2] + + g[b[0],b[1],b1[2]] += c[0]*c[1]*a[2] + g[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2] + g[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2] + g[b1[0],b1[1],b1[2]] += a[0]*a[1]*a[2] + + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef void INTERNAL_project_cic_with_mass(npx.ndarray[DTYPE_t, ndim=3] g, + npx.ndarray[DTYPE_t, ndim=2] x, + npx.ndarray[DTYPE_t, ndim=1] mass, + int Ngrid, double Lbox): + cdef double half_Box = 0.5*Lbox + cdef double delta_Box = Ngrid/Lbox + cdef int i + cdef double a[3], c[3] + cdef DTYPE_t m0 + cdef int b[3] + + for i in range(x.shape[0]): + + do_not_put = False + for j in range(3): + a[j] = (x[i,j]+half_Box)*delta_Box + b[j] = int(a[j]) + a[j] -= b[j] + c[j] = 1-a[j] + if b[j] < 0 or b[j]+1 >= Ngrid: + do_not_put = True + + if not do_not_put: + m0 = mass[i] + + g[b[0],b[1],b[2]] += c[0]*c[1]*c[2]*m0 + g[b[0]+1,b[1],b[2]] += a[0]*c[1]*c[2]*m0 + g[b[0],b[1]+1,b[2]] += c[0]*a[1]*c[2]*m0 + g[b[0]+1,b[1]+1,b[2]] += a[0]*a[1]*c[2]*m0 + + g[b[0],b[1],b[2]+1] += c[0]*c[1]*a[2]*m0 + g[b[0]+1,b[1],b[2]+1] += a[0]*c[1]*a[2]*m0 + g[b[0],b[1]+1,b[2]+1] += c[0]*a[1]*a[2]*m0 + g[b[0]+1,b[1]+1,b[2]+1] += a[0]*a[1]*a[2]*m0 + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef void INTERNAL_project_cic_with_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g, + npx.ndarray[DTYPE_t, ndim=2] x, + npx.ndarray[DTYPE_t, ndim=1] mass, + int Ngrid, double Lbox): + cdef double half_Box = 0.5*Lbox, m0 + cdef double delta_Box = Ngrid/Lbox + cdef int i + cdef double a[3], c[3] + cdef int b[3], b1[3] + + for i in range(x.shape[0]): + + for j in range(3): + a[j] = (x[i,j]+half_Box)*delta_Box + b[j] = int(floor(a[j])) + b1[j] = b[j]+1 + while b1[j] < 0: + b1[j] += Ngrid + while b1[j] >= Ngrid: + b1[j] -= Ngrid + + a[j] -= b[j] + c[j] = 1-a[j] + + m0 = mass[i] + g[b[0],b[1],b[2]] += c[0]*c[1]*c[2]*m0 + g[b1[0],b[1],b[2]] += a[0]*c[1]*c[2]*m0 + g[b[0],b1[1],b[2]] += c[0]*a[1]*c[2]*m0 + g[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2]*m0 + + g[b[0],b[1],b1[2]] += c[0]*c[1]*a[2]*m0 + g[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2]*m0 + g[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2]*m0 + g[b1[0],b1[1],b1[2]] += a[0]*a[1]*a[2]*m0 + + +def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, ndim=1] mass, int Ngrid, + double Lbox, bool periodic = False): + """ + project_cic(x array (N,3), mass (may be None), Ngrid, Lbox, periodict) + + This function does a Cloud-In-Cell projection of a 3d unstructured dataset. First argument is a Nx3 array of coordinates. + Second argument is an optinal mass. Ngrid is the size output grid and Lbox is the physical size of the grid. + """ + cdef npx.ndarray[DTYPE_t, ndim=3] g + + if x.shape[1] != 3: + raise ValueError("Invalid shape for x array") + + if mass != None and mass.shape[0] != x.shape[0]: + raise ValueError("Mass array and coordinate array must have the same number of elements") + + g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE) + + if not periodic: + if mass == None: + INTERNAL_project_cic_no_mass(g, x, Ngrid, Lbox) + else: + INTERNAL_project_cic_with_mass(g, x, mass, Ngrid, Lbox) + else: + if mass == None: + INTERNAL_project_cic_no_mass_periodic(g, x, Ngrid, Lbox) + else: + INTERNAL_project_cic_with_mass_periodic(g, x, mass, Ngrid, Lbox) + + return g + +def tophat_fourier_internal(npx.ndarray[DTYPE_t, ndim=1] x not None): + cdef int i + cdef npx.ndarray[DTYPE_t] y + cdef DTYPE_t x0 + + y = np.empty(x.size, dtype=DTYPE) + + for i in range(x.size): + x0 = x[i] + if abs(x0)<1e-5: + y[i] = 1 + else: + y[i] = (3*(sin(x0) - x0 * cos(x0))/(x0**3)) + + return y + +def tophat_fourier(x not None): + cdef npx.ndarray[DTYPE_t, ndim=1] b + + if type(x) != np.ndarray: + raise ValueError("x must be a Numpy array") + + b = np.array(x, dtype=DTYPE).ravel() + + b = tophat_fourier_internal(b) + + return b.reshape(x.shape) + + + +@cython.boundscheck(False) +@cython.cdivision(True) +cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]): + cdef DTYPE_t alpha_max + cdef DTYPE_t tmp_a + cdef DTYPE_t v[3] + cdef int i, j + + alpha_max = 10.0 # A big number + + for i in range(3): + if u[i] == 0.: + continue + + if u[i] < 0: + tmp_a = -u0[i]/u[i] + else: + tmp_a = (1-u0[i])/u[i] + + if tmp_a < alpha_max: + alpha_max = tmp_a + j = i + + for i in range(3): + u0[i] += u[i]*alpha_max + + r[0] = j + + return alpha_max + +#cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]): +# cdef DTYPE_t alpha_max +# cdef DTYPE_t 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.: +# continue + +# if u[i] < 0: +# tmp_a = -u0[i]/u[i] +# else: +# tmp_a = (1-u0[i])/u[i] + +# if tmp_a < alpha_max: +# alpha_max = tmp_a +# j = i + + # alpha_max is the integration length + # now we compute the integration of a trilinearly interpolated field + # There are four terms. + + + # First term +# term[0]= (u0[0]*u0[1]*u0[2])*sum(vertex_value) + + # Second term +# term[1] = 0 +# +# for q in range(3): +# for r in range(8): +# pass + +# for i in range(3): +# u0[i] += u[i]*alpha_max + +# r[0] = j + +# return 0#alpha_max + +@cython.boundscheck(False) +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): + + cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3] + cdef int iu0[3] + cdef int i + cdef int N = density.shape[0] + cdef int half_N = density.shape[0]/2 + cdef int completed + cdef DTYPE_t I0, d, dist2, delta, s, max_distance2 + cdef int jumper[1] + + max_distance2 = max_distance**2 + + for i in range(3): + u[i] = a_u[i] + u0[i] = a_u[i]*min_distance + ifu0[i] = half_N+u0[i] + if (ifu0[i] <= 0 or ifu0[i] >= N): + return 0 + iu0[i] = int(floor(ifu0[i])) + u0[i] = ifu0[i]-iu0[i] + 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)): + raise RuntimeError("u0[%d] = %g !" % (i,u0[i])) + + completed = 0 + if ((int(iu0[0]) >= N) or (int(iu0[0]) <= 0) or + (int(iu0[1]) >= N) or (int(iu0[1]) <= 0) or + (int(iu0[2]) >= N) or (int(iu0[2]) <= 0)): + completed = 1 + + I0 = 0 + jumper[0] = 0 + while completed == 0: + + d = density[iu0[0], iu0[1], iu0[2]] + s = cube_integral(u, u0, jumper) + I0 += s*d + + if u[jumper[0]] < 0: + iu0[jumper[0]] -= 1 + u0[jumper[0]] = 1 + else: + iu0[jumper[0]] += 1 + u0[jumper[0]] = 0 + + + if ((int(iu0[0]) >= N) or (int(iu0[0]) <= 0) or + (int(iu0[1]) >= N) or (int(iu0[1]) <= 0) or + (int(iu0[2]) >= N) or (int(iu0[2]) <= 0)): + completed = 1 + else: + dist2 = 0 + for i in range(3): + delta = iu0[i]+u0[i]-half_N + dist2 += delta*delta + + if (dist2 > max_distance2): + # Remove the last portion of the integral + delta = sqrt(dist2) - max_distance + I0 -= d*delta + completed = 1 + + return I0 + + +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): + + import healpy as hp + import progressbar as pb + cdef int i + cdef npx.ndarray[DTYPE_t, ndim=1] u + cdef npx.ndarray[DTYPE_t, ndim=1] outm + + outm = np.empty(hp.nside2npix(Nside),dtype=DTYPE) + + if progress: + p = pb.ProgressBar(maxval=outm.size).start() + + for i in range(outm.size): + 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) + + if progress: + p.finish() + + + return outm diff --git a/python/cosmotool/__init__.py b/python/cosmotool/__init__.py index de8244c..7a1fad3 100644 --- a/python/cosmotool/__init__.py +++ b/python/cosmotool/__init__.py @@ -1,4 +1,5 @@ from _cosmotool import * +from _project import * from grafic import writeGrafic, writeWhitePhase, readGrafic from borg import read_borg_vol from cic import cicParticles diff --git a/python/cosmotool/grafic.py b/python/cosmotool/grafic.py index 13d93ac..7eb5f41 100644 --- a/python/cosmotool/grafic.py +++ b/python/cosmotool/grafic.py @@ -42,8 +42,9 @@ def writeGrafic(filename, field, BoxSize, scalefac, **cosmo): bad, bad, bad, scalefac, cosmo['omega_M_0'], cosmo['omega_lambda_0'], 100*cosmo['h'], checkPoint)) - checkPoint = 4*Ny*Nz - for i in xrange(Nx): + checkPoint = 4*Ny*Nx + field = field.reshape(field.shape, order='F') + for i in xrange(Nz): f.write(struct.pack("I", checkPoint)) f.write(field[i].astype(np.float32).tostring()) f.write(struct.pack("I", checkPoint))