From cef0fb3786985d859c7ea0197e795c2d8b896c52 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Thu, 12 Jun 2014 17:16:59 +0200
Subject: [PATCH] Added a lean CIC filter

---
 python/_cosmotool.pyx                    | 46 ++++++++++++++++++++++++
 python/_project.pyx                      |  2 +-
 python_sample/icgen/test_ic_from_borg.py |  2 +-
 src/cic.cpp                              |  2 ++
 src/cic.hpp                              |  4 ++-
 5 files changed, 53 insertions(+), 3 deletions(-)

diff --git a/python/_cosmotool.pyx b/python/_cosmotool.pyx
index 0181b47..bd89d4e 100644
--- a/python/_cosmotool.pyx
+++ b/python/_cosmotool.pyx
@@ -3,6 +3,7 @@ from libcpp cimport string as cppstring
 import numpy as np
 cimport numpy as np
 from cpython cimport PyObject, Py_INCREF
+cimport cython
 
 np.import_array()
 
@@ -47,6 +48,22 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
     void setNormalization(double)
     double power(double)
 
+cdef extern from "cic.hpp" namespace "CosmoTool":
+
+  ctypedef float CICType
+  ctypedef float Coordinates[3]
+  
+  cdef cppclass CICParticles:
+    float mass
+    Coordinates coords
+
+  cdef cppclass CICFilter:
+  
+    CICFilter(np.uint32_t resolution, double L) nogil
+    
+    void resetMesh() nogil
+    void putParticles(CICParticles* particles, np.uint32_t N) nogil
+    void getDensityField(CICType *& field, np.uint32_t& res) nogil
 
 cdef extern from "loadSimu.hpp" namespace "CosmoTool":
 
@@ -400,3 +417,32 @@ cdef class CosmologyPower:
       return out
     else:
       return self._compute(k)
+      
+@cython.boundscheck(False)
+def leanCic(float[:,:] particles, float L, int Resolution):
+  cdef CICParticles p
+  cdef CICFilter *cic
+  cdef np.uint64_t i
+  cdef CICType *field
+  cdef np.uint32_t dummyRes
+
+  cic = new CICFilter(Resolution, L)
+  cic.resetMesh()
+  
+  if particles.shape[1] != 3:
+    raise ValueError("Particles must be Nx3 array")
+  
+  p.mass = 1
+  for i in xrange(particles.shape[0]):
+    p.coords[0] = particles[i,0]
+    p.coords[1] = particles[i,1]
+    p.coords[2] = particles[i,2]
+    cic.putParticles(&p, 1)
+    
+  field = <CICType*>0
+  dummyRes = 0
+  cic.getDensityField(field, dummyRes)
+  
+  
+  
+  del cic
diff --git a/python/_project.pyx b/python/_project.pyx
index 769c89e..29df4c8 100644
--- a/python/_project.pyx
+++ b/python/_project.pyx
@@ -748,7 +748,7 @@ cdef npx.uint64_t _mysum(int[:] jobs) nogil:
 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):
+                         DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=1000):
                          
     import healpy as hp
     import progressbar as pb
diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py
index 0ce037e..4da87a6 100644
--- a/python_sample/icgen/test_ic_from_borg.py
+++ b/python_sample/icgen/test_ic_from_borg.py
@@ -19,7 +19,7 @@ snap_id=int(sys.argv[1])
 astart=1/100.
 
 if doSimulation:
-  s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim2/", snap_id, doublePrecision=True)
+  s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim3/", snap_id, doublePrecision=True)
   astart=s.getTime()
   L = s.getBoxsize()
 
diff --git a/src/cic.cpp b/src/cic.cpp
index bfec994..8be6016 100644
--- a/src/cic.cpp
+++ b/src/cic.cpp
@@ -38,6 +38,8 @@ knowledge of the CeCILL license and that you accept its terms.
 #include <inttypes.h>
 #include "cic.hpp"
 
+using namespace CosmoTool;
+
 CICFilter::CICFilter(uint32_t N, double len)
 {
   spatialLen = len;
diff --git a/src/cic.hpp b/src/cic.hpp
index f76ec81..2792426 100644
--- a/src/cic.hpp
+++ b/src/cic.hpp
@@ -39,7 +39,7 @@ knowledge of the CeCILL license and that you accept its terms.
 #include "config.hpp"
 #include <inttypes.h>
 
-using namespace CosmoTool;
+namespace CosmoTool {
 
 typedef float CICType;
 
@@ -67,4 +67,6 @@ typedef float CICType;
     uint32_t szGrid;
   };
 
+};
+
 #endif