Imported the projection library for python
This commit is contained in:
Guilhem Lavaux
parent
ae040d15ee
commit
1693c222e1
@ -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
|
||||
)
|
||||
|
||||
|
||||
664
python/_project.pyx
Normal file
664
python/_project.pyx
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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))
|
||||
|
||||
Loading…
Reference in New Issue
Block a user