Fixes

python/CMakeLists.txt

@@ -11,7 +11,7 @@ add_custom_command(
 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) 
+  DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx ${CMAKE_CURRENT_SOURCE_DIR}/project_tool.hpp ) 
 
 ENDIF(CYTHON)
 

python/_project.pyx

@@ -1,4 +1,5 @@
 from cpython cimport bool
+from cython cimport view
 from libc.math cimport sin, cos, abs, floor, sqrt
 import numpy as np
 cimport numpy as npx
@@ -6,6 +7,7 @@ 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"]
 
@@ -541,7 +543,7 @@ cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
 @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 tmp_a
+    cdef DTYPE_t I, tmp_a
     cdef DTYPE_t v[3], term[4]
     cdef int i, j, q
 
@@ -561,15 +563,21 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
             alpha_max = tmp_a
             j = i
 
-    for i in range(3):
+
+    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 compute_projection(vertex_value, u, u0, alpha_max)
+    return I
-
 
 @cython.boundscheck(False)
 cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
@@ -592,14 +600,12 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
 #        a[i][0] = iu0[i]-1
 #        a[i][1] = iu0[i]
 #      else:
-      a[i][0] = iu0[i]-1
+      a[i][0] = iu0[i]
-      a[i][1] = iu0[i]
+      a[i][1] = iu0[i]+1
       
       with gil:    
         assert a[i][0] >= 0
         assert a[i][1] < density.shape[i]
-        assert u0[i] >=0
-        assert u0[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]]
@@ -619,7 +625,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
 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, int integrator_id=0):
+                             DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id=0):
 
     cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3]
     cdef int u_delta[3]
@@ -644,7 +650,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
     for i in range(3):
         u[i] = a_u[i]
         u0[i] = a_u[i]*min_distance
-        ifu0[i] = half_N+u0[i]
+        ifu0[i] = half_N+u0[i]+shifter[i]
         if (ifu0[i] <= 0 or ifu0[i] >= N):
           return 0
         iu0[i] = int(floor(ifu0[i]))
@@ -697,7 +703,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
 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 max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None):
                          
     import healpy as hp
     import progressbar as pb
@@ -705,6 +711,10 @@ def spherical_projection(int Nside,
     cdef npx.ndarray[DTYPE_t, ndim=1] u
     cdef npx.ndarray[DTYPE_t, ndim=1] outm
     
+    if shifter is None:
+        shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
+        shifter[:] = 0
+        
     outm = np.empty(hp.nside2npix(Nside),dtype=DTYPE)
 
     if progress:
@@ -714,7 +724,7 @@ def spherical_projection(int Nside,
         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, integrator_id=integrator_id)
+        outm[i] = line_of_sight_projection(density, u, min_distance, max_distance, shifter, integrator_id=integrator_id)
 
     if progress:
       p.finish()

python/project_tool.hpp

@@ -12,7 +12,7 @@ static T project_tool(T *vertex_value, T *u, T *u0)
       T ret = 0;
       
       for (int q = 0; q < 3; q++)
-        epsilon[q] = -(2*c[q]-1);
+        epsilon[q] = (2*c[q]-1);
 
       for (int q = 0; q < ProdType::numProducts; q++)
         ret += ProdType::product(u, u0, epsilon, q);
@@ -25,9 +25,9 @@ static T project_tool(T *vertex_value, T *u, T *u0)
 
 
 template<typename T>
-static T get_u0(T u0, int epsilon)
+static T get_u0(const T& u0, int epsilon)
 {
-  return (epsilon < 0) ? u0 : (1-u0);
+  return (epsilon > 0) ? u0 : (1-u0);
 }
 
 template<typename T>
@@ -108,7 +108,7 @@ T compute_projection(T *vertex_value, T *u, T *u0, T rho)
   T ret;
   
   ret = project_tool<T, ProductTerm0<T> >(vertex_value, u, u0) * rho;
-  ret += project_tool<T, ProductTerm1<T> >(vertex_value, u, u0) * rho * rho / 2;
+//  ret += project_tool<T, ProductTerm1<T> >(vertex_value, u, u0) * rho * rho / 2;
 //  ret += project_tool<T, ProductTerm2<T> >(vertex_value, u, u0) * rho * rho * rho / 3;
 //  ret += project_tool<T, ProductTerm3<T> >(vertex_value, u, u0) * rho * rho * rho * rho / 4;
   return ret;