pysphereproj/ext_src/_project.pyx

from cpython cimport bool
from cython cimport view
from cython.parallel import prange, parallel
from libc.math cimport sin, cos, abs, floor, round, 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

cdef extern from "openmp.hpp" namespace "sphereproj":
  int smp_get_max_threads() nogil
  int smp_get_thread_id() nogil


@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef void interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
                                        DTYPE_t z,
                                        DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) noexcept 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)

    ix = int(floor(rx))
    iy = int(floor(ry))
    iz = int(floor(rz))


    rx -= ix
    ry -= iy
    rz -= iz

    ix = ix % Ngrid
    iy = iy % Ngrid
    iz = iz % Ngrid

    jx = (ix+1)%Ngrid
    jy = (iy+1)%Ngrid
    jz = (iz+1)%Ngrid

    ix = ix%Ngrid
    iy = iy%Ngrid
    iz = iz%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)

    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]

@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef void ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
                                 DTYPE_t z,
                                 DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) noexcept 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)

    ix = int(round(rx))
    iy = int(round(ry))
    iz = int(round(rz))


    ix = ix%Ngrid
    iy = iy%Ngrid
    iz = iz%Ngrid

    retval[0] = d[ix ,iy ,iz ]


@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef void ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
                         DTYPE_t z,
                         DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) noexcept 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)

    ix = int(round(rx))
    iy = int(round(ry))
    iz = int(round(rz))

    if ((ix < 0)  or (ix+1) >= Ngrid or (iy < 0)  or (iy+1) >= Ngrid or (iz < 0)  or (iz+1) >= Ngrid):
      retval[0] = inval
      return

    retval[0] = d[ix ,iy ,iz ]


@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef void interp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
                               DTYPE_t z,
                               DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) noexcept 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

    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 or (iy < 0)  or (iy+1) >= Ngrid or (iz < 0)  or (iz+1) >= Ngrid):
      retval[0] = inval
      return
    #   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)

    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]

@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,
             bool periodic=False, bool centered=True, bool ngp=False, DTYPE_t inval = 0):
    """ interp3d(x,y,z,d,Lbox,periodic=False,centered=True,ngp=False) -> interpolated values

    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
    cdef DTYPE_t[:] out_slice
    cdef DTYPE_t[:] ax, ay, az
    cdef DTYPE_t[:,:,:] in_slice
    cdef DTYPE_t retval
    cdef long i
    cdef long Nelt
    cdef int myperiodic, myngp
    cdef DTYPE_t shifter

    myperiodic = periodic
    myngp = ngp

    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")


    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)
        out_slice = out
        in_slice = d
        Nelt = ax.size
        with nogil:
          if not myngp:
            if myperiodic:
              for i in prange(Nelt):
                interp3d_INTERNAL_periodic(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i])
            else:
              for i in prange(Nelt):
                interp3d_INTERNAL(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i], inval)
          else:
            if myperiodic:
              for i in prange(Nelt):
                ngp3d_INTERNAL_periodic(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i])
            else:
              for i in prange(Nelt):
                ngp3d_INTERNAL(shifter+ax[i], shifter+ay[i], shifter+az[i], in_slice, Lbox, &out_slice[i], inval)
        return out
    else:
      if not myngp:
        if periodic:
          interp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval)
        else:
          interp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval, inval)
      else:
        if periodic:
          ngp3d_INTERNAL_periodic(shifter+x, shifter+y, shifter+z, d, Lbox, &retval)
        else:
          ngp3d_INTERNAL(shifter+x, shifter+y, shifter+z, d, Lbox, &retval, inval)
      return retval

@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(DTYPE_t[:,:,:] g,
                                       DTYPE_t[:,:] x, int Ngrid, double Lbox, double shifter) noexcept nogil:
    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 = 0
        for j in range(3):
            a[j] = (x[i,j]+shifter)*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

        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(DTYPE_t[:,:,:] g,
                                                DTYPE_t[:,:] x, int Ngrid, double Lbox, double shifter) noexcept nogil:
    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
    cdef DTYPE_t[:,:] ax
    cdef DTYPE_t[:,:,:] ag

    ax = x
    ag = g

    for i in range(x.shape[0]):

        do_not_put = 0
        for j in range(3):
            a[j] = (ax[i,j]+shifter)*delta_Box
            b[j] = int(floor(a[j]))
            b1[j] = (b[j]+1) % Ngrid

            a[j] -= b[j]
            c[j] = 1-a[j]

            b[j] %= Ngrid

        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]

        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(DTYPE_t[:,:,:] g,
                                         DTYPE_t[:,:] x,
                                         DTYPE_t[:] mass,
                                         int Ngrid, double Lbox, double shifter) noexcept nogil:
    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]+shifter)*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(DTYPE_t[:,:,:] g,
                                                  DTYPE_t[:,:] x,
                                                  DTYPE_t[:] mass,
                                                  int Ngrid, double Lbox, double shifter) noexcept nogil:
    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]+shifter)*delta_Box
            b[j] = int(floor(a[j]))
            b1[j] = (b[j]+1) % Ngrid

            a[j] -= b[j]
            c[j] = 1-a[j]

            b[j] %= Ngrid


        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, centered=True, output=None):
    """
    project_cic(x array (N,3), mass (may be None), Ngrid, Lbox, periodict, centered=True, output=None)

    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.
    if output is not None, it must be a numpy array with dimension NgridxNgridxNgrid. The result will be accumulated in that array.
    """
    cdef npx.ndarray[DTYPE_t, ndim=3] g
    cdef double shifter
    cdef bool local_periodic
    cdef DTYPE_t[:] d_mass

    local_periodic = periodic

    if centered:
      shifter = 0.5*Lbox
    else:
      shifter = 0

    if x.shape[1] != 3:
        raise ValueError("Invalid shape for x array")

    if mass is not None and mass.shape[0] != x.shape[0]:
        raise ValueError("Mass array and coordinate array must have the same number of elements")

    if output is None:
        g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE)
    else:
        if type(output) != np.ndarray:
            raise ValueError("Invalid array type")
        g = output

    cdef DTYPE_t[:,:,:] d_g = g
    cdef DTYPE_t[:,:] d_x = x

    if not local_periodic:
       if mass is None:
         with nogil:
           INTERNAL_project_cic_no_mass(d_g, d_x, Ngrid, Lbox, shifter)
       else:
         d_mass = mass
         with nogil:
           INTERNAL_project_cic_with_mass(d_g, d_x, d_mass, Ngrid, Lbox, shifter)
    else:
       if mass is None:
         with nogil:
          INTERNAL_project_cic_no_mass_periodic(d_g, d_x, Ngrid, Lbox, shifter)
       else:
         d_mass = mass
         with nogil:
          INTERNAL_project_cic_with_mass_periodic(d_g, d_x, d_mass, Ngrid, Lbox, shifter)

    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], DTYPE_t alpha_max) nogil:
    cdef DTYPE_t tmp_a
    cdef DTYPE_t v[3]
    cdef int i, j

    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], DTYPE_t alpha_max) nogil:
    cdef DTYPE_t I, tmp_a
    cdef DTYPE_t v[3]
    cdef DTYPE_t term[4]
    cdef int i, j, q

    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

    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], DTYPE_t alpha_max) noexcept nogil:
    cdef DTYPE_t d

    d = density[iu0[0], iu0[1], iu0[2]]

    return cube_integral(u, u0, jumper, alpha_max)*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], DTYPE_t alpha_max) noexcept nogil:
    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

    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_trilin(u, u0, jumper, vertex_value, alpha_max)



@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) except? 0 nogil:

    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 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], DTYPE_t alpha_max) noexcept 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-1) or (iu0[0] <= 0) or
        (iu0[1] >= N-1) or (iu0[1] <= 0) or
        (iu0[2] >= N-1) or (iu0[2] <= 0)):
      completed = 1

    I0 = 0
    jumper[0] = 0
    dist2 = 0
    while completed == 0:
        I0 += integrator(density, u, u0, u_delta, iu0, jumper, max_distance-sqrt(dist2))

        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-1) or (iu0[0] <= 0) or
            (iu0[1] >= N-1) or (iu0[1] <= 0) or
            (iu0[2] >= N-1) 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]

  u[0] = a_u[0]
  u[1] = a_u[1]
  u[2] = a_u[2]

  return 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)
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):
    """
    spherical_projection(Nside, density, min_distance, max_distance, progress=1, integrator_id=0, shifter=None, booster=-1)

    Keyword arguments:
      progress (int): show progress if it is equal to 1
      integrator_id (int): specify the order of integration along the line of shift
      shifter (DTYPE_t array): this is an array of size 3. It specifies the amount of shift to apply to the center, in unit of voxel
      booster (int): what is the frequency of refreshment of the progress bar. Small number decreases performance by locking the GIL.


    Arguments:
      Nside (int): Nside of the returned map
      density (NxNxN array): this is the density field, expressed as a cubic array
      min_distance (float): lower bound of the integration
      max_distance (float): upper bound of the integration

    Returns:
      an healpix map, as a 1-dimensional array.
    """
    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, N0
    cdef double stheta
    cdef int tid

    if shifter is None:
        shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
        shifter[:] = 0

    print("allocating map")
    outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE)
    print("initializing views")
    outm = outm_array
    density_view = density

    print("progress?")
    if progress != 0:
      p = pb.ProgressBar(maxval=outm.size,widgets=[pb.Bar(), pb.ETA()]).start()

    N = smp_get_max_threads()
    N0 = outm.size

    if booster < 0:
      booster = 1#000

    job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int))
    job_done[:] = 0
    theta,phi = hp.pix2ang(Nside, np.arange(N0))
    with nogil, parallel():
      tid = smp_get_thread_id()
      for i in prange(N0,schedule='dynamic',chunksize=256):
          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[tid] += 1

    if progress:
      p.finish()


    return outm_array