cosmotool/python/_project.pyx

882 lines
25 KiB
Cython
Raw Normal View History

from cpython cimport bool
2014-06-05 16:24:05 +02:00
from cython cimport view
from cython.parallel import prange, parallel
2014-07-18 15:54:12 +02:00
from libc.math cimport sin, cos, abs, floor, round, sqrt
import numpy as np
cimport numpy as npx
cimport cython
2014-06-13 17:38:49 +02:00
from copy cimport *
2014-06-14 17:15:56 +02:00
from openmp cimport omp_get_max_threads, omp_get_thread_num
ctypedef npx.float64_t DTYPE_t
DTYPE=np.float64
2014-06-05 16:24:05 +02:00
FORMAT_DTYPE="d"
__all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"]
2014-06-04 18:06:48 +02:00
cdef extern from "project_tool.hpp" namespace "":
DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho) nogil
@cython.boundscheck(False)
@cython.cdivision(True)
2014-07-10 16:19:40 +02:00
cdef int interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
DTYPE_t z,
2014-07-10 16:19:40 +02:00
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
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
2014-07-10 16:19:40 +02:00
rx = (inv_delta*x)
ry = (inv_delta*y)
rz = (inv_delta*z)
ix = int(floor(rx))
iy = int(floor(ry))
iz = int(floor(rz))
rx -= ix
ry -= iy
rz -= iz
jx = (ix+1)%Ngrid
jy = (iy+1)%Ngrid
jz = (iz+1)%Ngrid
ix = ix%Ngrid
iy = iy%Ngrid
iz = iz%Ngrid
2014-07-09 13:58:48 +02:00
if (ix < 0) or (jx >= Ngrid):
2014-07-10 16:19:40 +02:00
return -1
2014-07-09 13:58:48 +02:00
if (iy < 0) or (jy >= Ngrid):
2014-07-10 16:19:40 +02:00
return -2
2014-07-09 13:58:48 +02:00
if (iz < 0) or (jz >= Ngrid):
2014-07-10 16:19:40 +02:00
return -3
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)
2014-07-10 16:19:40 +02:00
retval[0] = \
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]
2014-07-10 16:19:40 +02:00
return 0
2014-07-18 14:34:01 +02:00
@cython.boundscheck(False)
@cython.cdivision(True)
cdef int ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
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)
ry = (inv_delta*y)
rz = (inv_delta*z)
2014-07-18 15:54:12 +02:00
ix = int(round(rx))
iy = int(round(ry))
iz = int(round(rz))
2014-07-18 14:34:01 +02:00
ix = ix%Ngrid
iy = iy%Ngrid
iz = iz%Ngrid
retval[0] = d[ix ,iy ,iz ]
return 0
@cython.boundscheck(False)
@cython.cdivision(True)
cdef int ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
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)
ry = (inv_delta*y)
rz = (inv_delta*z)
2014-07-18 15:54:12 +02:00
ix = int(round(rx))
iy = int(round(ry))
iz = int(round(rz))
2014-07-18 14:34:01 +02:00
2014-07-18 14:36:55 +02:00
if (ix < 0 or ix >= Ngrid):
2014-07-18 14:34:01 +02:00
return -1
2014-07-18 14:36:55 +02:00
if (iy < 0 or iy >= Ngrid):
2014-07-18 14:34:01 +02:00
return -2
2014-07-18 14:36:55 +02:00
if (iz < 0 or iz >= Ngrid):
2014-07-18 14:34:01 +02:00
return -3
retval[0] = d[ix ,iy ,iz ]
return 0
@cython.boundscheck(False)
@cython.cdivision(True)
2014-07-10 16:19:40 +02:00
cdef int interp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
DTYPE_t z,
2014-07-10 16:19:40 +02:00
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
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
2014-07-10 16:19:40 +02:00
rx = (inv_delta*x)
ry = (inv_delta*y)
rz = (inv_delta*z)
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):
2014-07-10 16:19:40 +02:00
return -1
if ((iy < 0) or (iy+1) >= Ngrid):
2014-07-10 16:19:40 +02:00
return -2
if ((iz < 0) or (iz+1) >= Ngrid):
2014-07-10 16:19:40 +02:00
return -3
# 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)
2014-07-10 16:19:40 +02:00
retval[0] = \
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]
2014-07-10 16:19:40 +02:00
return 0
2014-07-09 13:58:48 +02:00
@cython.boundscheck(False)
def interp3d(x not None, y not None,
z not None,
npx.ndarray[DTYPE_t, ndim=3] d not None, DTYPE_t Lbox,
2014-07-18 14:34:01 +02:00
bool periodic=False, bool centered=True, bool ngp=False):
""" interp3d(x,y,z,d,Lbox,periodic=False,centered=True,ngp=False) -> interpolated values
2014-06-18 16:13:16 +02:00
Compute the tri-linear interpolation of the given field (d) at the given position (x,y,z). It assumes that they are box-centered coordinates. So (x,y,z) == (0,0,0) is equivalent to the pixel at (Nx/2,Ny/2,Nz/2) with Nx,Ny,Nz = d.shape. If periodic is set, it assumes the box is periodic
"""
cdef npx.ndarray[DTYPE_t] out
2014-07-09 13:58:48 +02:00
cdef DTYPE_t[:] out_slice
cdef DTYPE_t[:] ax, ay, az
2014-07-10 16:19:40 +02:00
cdef DTYPE_t[:,:,:] in_slice
cdef DTYPE_t retval
2014-07-09 13:58:48 +02:00
cdef long i
cdef long Nelt
2014-07-18 14:34:01 +02:00
cdef int myperiodic, myngp
2014-07-10 16:19:40 +02:00
cdef DTYPE_t shifter
2014-07-09 13:58:48 +02:00
myperiodic = periodic
2014-07-18 14:34:01 +02:00
myngp = ngp
2014-07-10 16:19:40 +02:00
if centered:
shifter = Lbox/2
else:
shifter = 0
if d.shape[0] != d.shape[1] or d.shape[0] != d.shape[2]:
raise ValueError("Grid must have a cubic shape")
2014-07-10 16:19:40 +02:00
ierror = IndexError("Interpolating outside range")
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)
2014-07-09 13:58:48 +02:00
out_slice = out
2014-07-10 16:19:40 +02:00
in_slice = d
2014-07-09 13:58:48 +02:00
Nelt = ax.size
with nogil:
2014-07-18 15:54:12 +02:00
if not myngp:
2014-07-18 14:34:01 +02:00
if myperiodic:
for i in prange(Nelt):
if interp3d_INTERNAL_periodic(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i]) < 0:
with gil:
raise ierror
else:
for i in prange(Nelt):
if interp3d_INTERNAL(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i]) < 0:
with gil:
raise ierror
2014-07-09 13:58:48 +02:00
else:
2014-07-18 14:34:01 +02:00
if myperiodic:
for i in prange(Nelt):
if ngp3d_INTERNAL_periodic(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i]) < 0:
with gil:
raise ierror
else:
for i in prange(Nelt):
if ngp3d_INTERNAL(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i]) < 0:
with gil:
raise ierror
return out
else:
2014-07-18 15:54:12 +02:00
if not myngp:
if periodic:
if interp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0:
raise ierror
else:
if interp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0:
raise ierror
else:
2014-07-18 15:54:12 +02:00
if periodic:
if ngp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0:
raise ierror
else:
if ngp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval) < 0:
raise ierror
2014-07-10 16:19:40 +02:00
return retval
2014-07-15 15:49:11 +02:00
@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
2014-06-14 17:15:56 +02:00
cdef double a[3], c[3]
cdef int b[3]
cdef int do_not_put
for i in range(x.shape[0]):
2014-07-10 16:19:40 +02:00
do_not_put = 0
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
2014-06-14 17:15:56 +02:00
cdef double a[3], c[3]
cdef int b[3], b1[3]
cdef int do_not_put
2014-07-15 15:49:11 +02:00
cdef DTYPE_t[:,:] ax
cdef DTYPE_t[:,:,:] ag
ax = x
ag = g
for i in range(x.shape[0]):
2014-07-10 16:19:40 +02:00
do_not_put = 0
for j in range(3):
2014-07-15 15:49:11 +02:00
a[j] = (ax[i,j]+half_Box)*delta_Box
b[j] = int(floor(a[j]))
2014-07-10 16:19:40 +02:00
b1[j] = (b[j]+1) % Ngrid
a[j] -= b[j]
c[j] = 1-a[j]
2014-07-10 16:19:40 +02:00
b[j] %= Ngrid
2014-07-15 15:49:11 +02:00
ag[b[0],b[1],b[2]] += c[0]*c[1]*c[2]
ag[b1[0],b[1],b[2]] += a[0]*c[1]*c[2]
ag[b[0],b1[1],b[2]] += c[0]*a[1]*c[2]
ag[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2]
2014-07-15 15:49:11 +02:00
ag[b[0],b[1],b1[2]] += c[0]*c[1]*a[2]
ag[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2]
ag[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2]
ag[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
2014-06-14 17:15:56 +02:00
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
2014-06-14 17:15:56 +02:00
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)
2014-06-04 18:06:48 +02:00
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 tmp_a
cdef DTYPE_t v[3]
cdef int i, j
alpha_max = 10.0 # A big number
2014-06-04 18:06:48 +02:00
for i in xrange(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
@cython.boundscheck(False)
@cython.cdivision(True)
2014-06-14 17:15:56 +02:00
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
2014-06-05 16:24:05 +02:00
cdef DTYPE_t I, tmp_a
2014-06-14 17:15:56 +02:00
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
2014-06-14 17:15:56 +02:00
2014-06-05 16:24:05 +02:00
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
2014-06-05 16:24:05 +02:00
return I
2014-06-04 18:06:48 +02:00
@cython.boundscheck(False)
cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
2014-06-14 17:15:56 +02:00
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil:
2014-06-04 18:06:48 +02:00
cdef DTYPE_t d
d = density[iu0[0], iu0[1], iu0[2]]
return cube_integral(u, u0, jumper)*d
@cython.boundscheck(False)
cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
2014-06-14 17:15:56 +02:00
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil:
2014-06-04 18:06:48 +02:00
cdef DTYPE_t vertex_value[8]
cdef DTYPE_t d
2014-06-14 17:15:56 +02:00
cdef int a[3][2], i
2014-06-04 18:06:48 +02:00
for i in xrange(3):
2014-06-05 16:24:05 +02:00
a[i][0] = iu0[i]
a[i][1] = iu0[i]+1
2014-06-14 17:15:56 +02:00
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]]
2014-06-14 17:15:56 +02:00
return cube_integral_trilin(u, u0, jumper, vertex_value)
2014-06-04 18:06:48 +02:00
@cython.boundscheck(False)
cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
DTYPE_t a_u[3],
DTYPE_t min_distance,
2014-06-14 17:15:56 +02:00
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 0:
2014-06-14 17:15:56 +02:00
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]
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]
2014-06-04 18:06:48 +02:00
cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
2014-06-14 17:15:56 +02:00
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1]) nogil
2014-06-04 18:06:48 +02:00
if integrator_id == 0:
integrator = integrator0
else:
integrator = integrator1
max_distance2 = max_distance**2
for i in range(3):
u[i] = a_u[i]
u0[i] = a_u[i]*min_distance
2014-06-05 16:24:05 +02:00
ifu0[i] = half_N+u0[i]+shifter[i]
if (ifu0[i] <= 0 or ifu0[i] >= N):
return 0
iu0[i] = int(floor(ifu0[i]))
u0[i] = ifu0[i]-iu0[i]
u_delta[i] = 1 if iu0[i] > 0 else -1
2014-06-14 17:15:56 +02:00
if (not ((iu0[i]>= 0) and (iu0[i] < N))):
with gil:
raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i]))
if (not (u0[i]>=0 and u0[i]<=1)):
with gil:
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)):
completed = 1
I0 = 0
jumper[0] = 0
while completed == 0:
2014-06-14 17:15:56 +02:00
I0 += integrator(density, u, u0, u_delta, iu0, jumper)
if u[jumper[0]] < 0:
iu0[jumper[0]] -= 1
u0[jumper[0]] = 1
else:
iu0[jumper[0]] += 1
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)):
completed = 1
else:
dist2 = 0
for i in range(3):
2014-06-13 17:38:49 +02:00
delta = iu0[i]+u0[i]-half_N-shifter[i]
dist2 += delta*delta
if (dist2 > max_distance2):
# Remove the last portion of the integral
2014-06-04 18:06:48 +02:00
#delta = sqrt(dist2) - max_distance
#I0 -= d*delta
completed = 1
return I0
def line_of_sight_projection(DTYPE_t[:,:,:] density not None,
DTYPE_t[:] a_u not None,
DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter not None, int integrator_id=0):
cdef DTYPE_t u[3]
u[0] = a_u[0]
u[1] = a_u[1]
u[2] = a_u[2]
2014-06-14 17:15:56 +02:00
C_line_of_sight_projection(density,
u,
min_distance,
2014-06-14 17:15:56 +02:00
max_distance, shifter, integrator_id)
cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density,
double min_distance, double max_distance,
2014-06-14 17:15:56 +02:00
DTYPE_t[:] shifter, int integrator_id) nogil:
cdef DTYPE_t u0[3]
stheta = sin(theta)
u0[0] = cos(phi)*stheta
u0[1] = sin(phi)*stheta
u0[2] = cos(theta)
2014-06-14 17:15:56 +02:00
return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id)
@cython.boundscheck(False)
def spherical_projection(int Nside,
npx.ndarray[DTYPE_t, ndim=3] density not None,
DTYPE_t min_distance,
2014-06-13 11:07:58 +02:00
DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=-1):
import healpy as hp
import progressbar as pb
cdef int i
cdef DTYPE_t[:] theta,phi
cdef DTYPE_t[:,:,:] density_view
cdef DTYPE_t[:] outm
cdef int[:] job_done
cdef npx.ndarray[DTYPE_t, ndim=1] outm_array
2014-06-14 17:15:56 +02:00
cdef long N, N0
cdef double stheta
2014-06-14 17:15:56 +02:00
cdef int tid
2014-06-05 16:24:05 +02:00
if shifter is None:
shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
shifter[:] = 0
outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE)
outm = outm_array
density_view = density
if progress != 0:
2014-06-12 15:35:54 +02:00
p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start()
N = omp_get_max_threads()
2014-06-14 17:15:56 +02:00
N0 = outm.size
if booster < 0:
booster = 1000
2014-06-12 15:35:54 +02:00
job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int))
job_done[:] = 0
2014-06-14 17:15:56 +02:00
theta,phi = hp.pix2ang(Nside, np.arange(N0))
with nogil, parallel():
2014-06-14 17:15:56 +02:00
tid = omp_get_thread_num()
for i in prange(N0):
if progress != 0 and (i%booster) == 0:
with gil:
p.update(_mysum(job_done))
2014-06-14 17:15:56 +02:00
outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id)
job_done[tid] += 1
if progress:
p.finish()
return outm_array