from cpython cimport bool from cython cimport view from cython.parallel import prange, parallel from openmp cimport omp_get_max_threads, omp_get_thread_num from libc.math cimport sin, cos, abs, floor, sqrt import numpy as np cimport numpy as npx cimport cython from copy cimport * 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] cdef double 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] cdef double c[3] cdef int b[3] cdef int 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] cdef double 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] cdef double c[3] cdef int b[3] cdef int 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 cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], int err[1]) nogil: cdef DTYPE_t alpha_max cdef DTYPE_t I, tmp_a cdef DTYPE_t v[3] cdef DTYPE_t 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 # if (u0[i] < 0 or u0[i] > 1): # err[0] = 1 # with gil: # print("Bouh! u0[%d] = %lg" % (i, u0[i])) # return 0 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], int err[1]) nogil except? 0: 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], int err[1]) nogil except? 0: cdef DTYPE_t vertex_value[8] cdef DTYPE_t d cdef int a[3][2] cdef int i for i in xrange(3): 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, err) @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, int out_err[1]) nogil except? 0: cdef DTYPE_t u[3] cdef DTYPE_t ifu0[3] cdef DTYPE_t u0[3] cdef DTYPE_t utot[3] cdef int u_delta[3] cdef int iu0[3] cdef int err[1] 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], int err[1]) nogil except? 0 out_err[0] = 0 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): with gil: print "BLABLA" 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))): # out_err[0] = 1 # return 0 # if (not (u0[i]>=0 and u0[i]<=1)): # out_err[0] = 1 # return 0 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: err[0] = 0 I0 += integrator(density, u, u0, u_delta, iu0, jumper, err) # if err[0] == 1: # with gil: # print("Bah! error!") # out_err[0] = 1 # return 0 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-shifter[i] 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] cdef int out_err[1] cdef DTYPE_t v u[0] = a_u[0] u[1] = a_u[1] u[2] = a_u[2] out_err[0] = 0 v = C_line_of_sight_projection(density, u, min_distance, max_distance, shifter, integrator_id, out_err) # if out_err[0] == 1: # raise RuntimeError("Error occured during integration") return v cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density, double min_distance, double max_distance, DTYPE_t[:] shifter, int integrator_id, int out_err[1]) 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, out_err) @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=-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 cdef long N cdef double stheta cdef int out_err[1] 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.Bar(), pb.ETA()]).start() N = omp_get_max_threads() job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int)) job_done[:] = 0 theta,phi = hp.pix2ang(Nside, np.arange(outm.size)) if booster <= 0: booster = density.size / 100 / N with nogil, parallel(): for i in prange(N): if omp_get_thread_num() == 0 and progress != 0 and (i%booster) == 0: with gil: p.update(_mysum(job_done)) # out_err[0] = 0 outm[i] = _spherical_projloop(theta[i], phi[i], density_view, min_distance, max_distance, shifter, integrator_id, out_err) job_done[omp_get_thread_num()] = i if progress: p.finish() return outm_array