more work on spherical projection. testlet for this code

python/_project.pyx

@@ -11,7 +11,7 @@ __all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE"
 
 cdef extern from "project_tool.hpp" namespace "":
 
-  DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho)
+  DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho) nogil
 
 
 @cython.boundscheck(False)
@@ -561,15 +561,19 @@ 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):
+        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)
 
 
 @cython.boundscheck(False)
 cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
-                         DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil:
+                         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]]
@@ -578,13 +582,37 @@ cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
 
 @cython.boundscheck(False)
 cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
-                         DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil:
+                         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
     
-    d = density[iu0[0], iu0[1], iu0[2]]
+    for i in xrange(3):
-    return cube_integral(u, u0, jumper)*d
+#      if u[i] < 0:
-    cube_integral_trilin(u, u0, jumper, vertex_value)
+#        a[i][0] = iu0[i]-1
+#        a[i][1] = iu0[i]
+#      else:
+      a[i][0] = iu0[i]-1
+      a[i][1] = iu0[i]
+      
+      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]]
+    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)
@@ -594,6 +622,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
                              DTYPE_t max_distance, int integrator_id=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]
@@ -603,7 +632,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
     cdef int jumper[1]
     
     cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
-                               DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil
+                               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
@@ -620,6 +649,7 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
           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))):
                raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i]))
         if (not (u0[i]>=0 and u0[i]<=1)):
@@ -635,15 +665,13 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
     jumper[0] = 0
     while completed == 0:
         
-        I0 += integrator(density, u, u0, iu0, jumper)
+        I0 += integrator(density, u, u0, u_delta, iu0, jumper)
 
         if u[jumper[0]] < 0:
             iu0[jumper[0]] -= 1
-            direction = -1
             u0[jumper[0]] = 1
         else:
             iu0[jumper[0]] += 1
-            direction = 1
             u0[jumper[0]] = 0
 
             
@@ -669,7 +697,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):
+                         DTYPE_t max_distance, int progress=1, int integrator_id=0):
                          
     import healpy as hp
     import progressbar as pb
@@ -686,7 +714,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)
+        outm[i] = line_of_sight_projection(density, u, min_distance, max_distance, 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);
@@ -24,6 +24,12 @@ static T project_tool(T *vertex_value, T *u, T *u0)
 }
 
 
+template<typename T>
+static T get_u0(T u0, int epsilon)
+{
+  return (epsilon < 0) ? u0 : (1-u0);
+}
+
 template<typename T>
 struct ProductTerm0
 {
@@ -34,11 +40,12 @@ struct ProductTerm0
     T a = 1;
     
     for (int r = 0; r < 3; r++)
-      a *= (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+      a *= get_u0(u0[r], epsilon[r]);
     return a;
   }
 };
 
+
 template<typename T>
 struct ProductTerm1
 {
@@ -51,12 +58,12 @@ struct ProductTerm1
     
     for (int r = 0; r < 3; r++)
       {
-        G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+        G[r] = get_u0(u0[r], epsilon[r]);
       }
 
-    double F[3] = { G[0]*u[1]*u[2], u[0]*G[1]*u[2], u[0]*u[1]*G[2] };    
+    double F[3] = { G[1]*G[2], G[0]*G[2], G[0]*G[1] };    
 
-    return F[q] * epsilon[q];
+    return F[q] * u[q] * epsilon[q];
   }
 };
 
@@ -72,16 +79,17 @@ struct ProductTerm2
     
     for (int r = 0; r < 3; r++)
       {
-        G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
+        G[r] = get_u0(u0[r], epsilon[r]);
       }
 
-    double F[3] = { u[0]*G[1]*G[2], G[0]*u[1]*G[2], G[0]*G[1]*u[2] };    
+    double F[3] = { epsilon[1]*epsilon[2]*u[1]*u[2], epsilon[0]*epsilon[2]*u[0]*u[2], epsilon[0]*epsilon[1]*u[0]*u[1] };    
 
-    return F[q] * epsilon[q];
+    return F[q] * G[q];
   }
 };
 
 
+
 template<typename T>
 struct ProductTerm3
 {
@@ -89,9 +97,7 @@ struct ProductTerm3
   
   static T product(T *u, T *u0, int *epsilon, int q)
   {
-    T a = 1;
+    return epsilon[0]*epsilon[1]*epsilon[2]*u[0]*u[1]*u[2];
-    
-    return epsilon[0]*epsilon[1]*epsilon[2];
   }
 };
 
@@ -103,8 +109,8 @@ T compute_projection(T *vertex_value, T *u, T *u0, T rho)
   
   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, ProductTerm2<T> >(vertex_value, u, u0) * rho * rho * rho / 3;
+//  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;
+//  ret += project_tool<T, ProductTerm3<T> >(vertex_value, u, u0) * rho * rho * rho * rho / 4;
   return ret;
 }
 

python_sample/test_spheric_proj.py

@@ -0,0 +1,10 @@
+import cosmotool as ct
+import numpy as np
+import healpy as hp
+
+d = np.zeros((64,64,64), ct.DTYPE)
+
+d[33,33,33] = 1
+
+proj0 = ct.spherical_projection(32, d, 0, 20, integrator_id=0)
+proj1 = ct.spherical_projection(32, d, 0, 20, integrator_id=1)