Added fast parallel interpolator. Some adjustment to BORG icgen

python/CMakeLists.txt

@@ -13,6 +13,11 @@ IF(CYTHON)
     COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx
     DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx) 
 
+  add_custom_command(
+    OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp
+    COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx
+    DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx) 
+
   add_custom_command(
     OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp
     COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_cic.pyx
@@ -29,6 +34,7 @@ ENDIF(CYTHON)
 add_library(_cosmotool MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp)
 add_library(_cosmo_power MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp)
 add_library(_cosmo_cic MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp)
+add_library(_fast_interp MODULE ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp)
 add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp)
 
 SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Bsymbolic-functions")
@@ -37,6 +43,7 @@ target_link_libraries(_cosmotool ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LI
 target_link_libraries(_cosmo_power ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LIBRARIES})
 target_link_libraries(_cosmo_cic ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LIBRARIES})
 target_link_libraries(_project ${PYTHON_LIBRARIES})
+target_link_libraries(_fast_interp ${CosmoTool_local} ${PYTHON_LIBRARIES})
 
 # Discover where to put packages
 if (NOT PYTHON_SITE_PACKAGES)
@@ -62,7 +69,7 @@ if (WIN32 AND NOT CYGWIN)
   SET_TARGET_PROPERTIES(_cosmotool PROPERTIES SUFFIX ".pyd")
 endif (WIN32 AND NOT CYGWIN)
 
-INSTALL(TARGETS _cosmotool _project _cosmo_power _cosmo_cic
+INSTALL(TARGETS _cosmotool _project _cosmo_power _cosmo_cic _fast_interp
   LIBRARY DESTINATION ${PYTHON_SITE_PACKAGES}/cosmotool
 )
 

python/_fast_interp.pyx

@@ -0,0 +1,32 @@
+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
+
+__all__=["fast_interp"]
+
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+def fast_interp(npx.float64_t xmin, npx.float64_t dx, npx.float64_t[:] A, npx.float64_t[:] y, npx.float64_t[:] out, npx.float64_t beyond_val=npx.nan):
+  cdef double rq, q
+  cdef int iq
+  cdef long i, Asize, ysize
+  cdef npx.float64_t[:] out0
+
+  ysize = y.size
+  out0 = out
+  Asize = A.size
+
+  with nogil:
+    for i in prange(ysize):
+      q = (y[i] - xmin) / dx
+      iq = int(floor(q))
+      rq = (q-iq)
+      if iq+1 >= Asize or iq < 0:
+        out0[i] = beyond_val
+      else:
+        out0[i] = rq * A[iq+1] + (1-rq)*A[iq]

python/cosmotool/__init__.py

@@ -2,6 +2,7 @@ from _cosmotool import *
 from _project import *
 from _cosmo_power import *
 from _cosmo_cic import *
+from _fast_interp import *
 from .grafic import writeGrafic, writeWhitePhase, readGrafic, readWhitePhase
 from .borg import read_borg_vol
 from .cic import cicParticles

python/cosmotool/fftw.py

@@ -8,7 +8,7 @@ class CubeFT(object):
   
       self.N = N
       self.align = pyfftw.simd_alignment
-      self.L = L
+      self.L = float(L)
       self.max_cpu = multiprocessing.cpu_count() if max_cpu < 0 else max_cpu
       self._dhat = pyfftw.n_byte_align_empty((self.N,self.N,self.N/2+1), self.align, dtype='complex64')
       self._density = pyfftw.n_byte_align_empty((self.N,self.N,self.N), self.align, dtype='float32')

python_sample/gen_2lpt_asmooth.py

@@ -7,7 +7,8 @@ import icgen.cosmogrowth as cg
 import sys
 import argparse
 
-ADAPT_SMOOTH="/home/bergeron1NS/lavaux/Software/cosmotool/build/sample/simple3DFilter"
+#ADAPT_SMOOTH="/home/bergeron1NS/lavaux/Software/cosmotool/build/sample/simple3DFilter"
+ADAPT_SMOOTH="/home/guilhem/PROJECTS/cosmotool/build/sample/simple3DFilter"
 cosmo={'omega_M_0':0.3175, 'h':0.6711}
 cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
 cosmo['omega_k_0'] = 0
@@ -32,13 +33,13 @@ args = parser.parse_args()
 
 
 
-for i in xrange(args.start, args.end, args.step):
+for i in [4629]:#xrange(args.start, args.end, args.step):
   print i
 
-  pos,vel,density,N,L,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
+  pos,vel,density,N,L,_,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
                                            cosmo, supersample=args.supersample, do_lpt2=True)
 
-  q = pos + vel
+  q = pos + vel + [np.ones(vel[0].shape[0])]
 
   with h5.File("particles.h5", mode="w") as f:
     f.create_dataset("particles", data=np.array(q).transpose())

python_sample/icgen/borgicgen.py

@@ -137,7 +137,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate
       p.normalize(cosmo['SIGMA8']*Dborg)
       density = do_supergenerate(density,mulfac=supergenerate,Pk=p,L=L,h=cosmo['h'])
 
-    lpt = LagrangianPerturbation(density, L, fourier=True, supersample=supersample)
+    lpt = LagrangianPerturbation(-density, L, fourier=True, supersample=supersample)
 
     # Generate grid
     posq = gen_posgrid(N*supergenerate, L)