guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
6967a1f55c
cosmotool/python/_project.pyx

792 lines
23 KiB
Cython
Raw Blame History

from cpython cimport bool
from cython cimport view
from cython.parallel import prange, parallel
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
FORMAT_DTYPE="d"
__all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"]
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)
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]) 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
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)
cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
cdef int i
cdef DTYPE_t s
s = 0
for i in xrange(q):
s += v[i]
return s
@cython.boundscheck(False)
@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 I, 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
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 I
@cython.boundscheck(False)
cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
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]]
return cube_integral(u, u0, jumper)*d
@cython.boundscheck(False)
cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
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
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]
a[i][1] = iu0[i]+1
with gil:
assert a[i][0] >= 0
assert a[i][1] < density.shape[i]
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)
cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
DTYPE_t a_u[3],
DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 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]
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]
cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
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
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
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
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:
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):
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 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]
C_line_of_sight_projection(density,
u,
min_distance,
max_distance, shifter, integrator_id)
cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density,
double min_distance, double max_distance,
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)
return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id)
@cython.boundscheck(False)
cdef npx.uint64_t _mysum(int[:] jobs) nogil:
cdef npx.uint64_t s
cdef npx.uint64_t N
cdef int i
s = 0
N = jobs.shape[0]
for i in xrange(N):
s += jobs[i]
return s
@cython.boundscheck(False)
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[:] shifter = None, int booster=10):
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
cdef long N
cdef double stheta
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:
p = pb.ProgressBar(maxval=outm.size,widgets=[pb.BouncingBar(), pb.ETA()]).start()
N = outm.size
job_done = view.array(shape=(1,), format="i", itemsize=sizeof(int))
job_done[:] = 0
theta,phi = hp.pix2ang(Nside, np.arange(N))
with nogil, parallel():
for i in prange(N):
if progress != 0 and (i%booster) == 0:
with gil:
p.update(_mysum(job_done))
outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id)
job_done[i] = 1
if progress:
p.finish()
return outm_array
Reference in New Issue View Git Blame Copy Permalink