Fixed CIC errors. Added support for CosmoPower in _cosmotool.pyx

python/_cosmotool.pyx

@@ -6,6 +6,47 @@ from cpython cimport PyObject, Py_INCREF
 
 np.import_array()
 
+cdef extern from "cosmopower.hpp" namespace "CosmoTool":
+
+  cdef enum CosmoFunction "CosmoTool::CosmoPower::CosmoFunction":
+      POWER_EFSTATHIOU "CosmoTool::CosmoPower::POWER_EFSTATHIOU",
+      HU_WIGGLES "CosmoTool::CosmoPower::HU_WIGGLES",
+      HU_BARYON "CosmoTool::CosmoPower::HU_BARYON",
+      OLD_POWERSPECTRUM,
+      POWER_BARDEEN "CosmoTool::CosmoPower::POWER_BARDEEN",
+      POWER_SUGIYAMA "CosmoTool::CosmoPower::POWER_SUGIYAMA",
+      POWER_BDM,
+      POWER_TEST
+
+  cdef cppclass CosmoPower:
+    double n
+    double K0
+    double V_LG_CMB
+    
+    double CMB_VECTOR[3]
+    double h
+    double SIGMA8
+    double OMEGA_B
+    double OMEGA_C
+    double omega_B
+    double omega_C
+    double Theta_27
+    double OMEGA_0
+    double Omega
+    double beta
+    double OmegaEff
+    double Gamma0
+    double normPower
+    
+        
+    CosmoPower()
+    void setFunction(CosmoFunction)
+    void updateCosmology()
+    void updatePhysicalCosmology()
+    void normalize()
+    void setNormalization(double)
+    double power(double)
+
 
 cdef extern from "loadSimu.hpp" namespace "CosmoTool":
 
@@ -302,4 +343,60 @@ def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision =
 
   return _PySimulationAdaptor(wrap_simudata(data, flags))
   
-def loadAllRamses(str basepath, int 
+ 
+cdef class CosmologyPower:
+
+  cdef CosmoPower power
+
+  def __init__(self,**cosmo):
+  
+    self.power = CosmoPower()
+    self.power.OMEGA_B = cosmo['omega_B_0']
+    self.power.OMEGA_C = cosmo['omega_M_0']-cosmo['omega_B_0']
+    self.power.h = cosmo['h']
+    if 'ns' in cosmo:
+      self.power.n = cosmo['ns'] 
+    
+    assert self.power.OMEGA_C > 0
+    
+    self.power.updateCosmology()
+    
+  def normalize(self,s8):
+    self.power.SIGMA8 = s8
+    self.power.normalize()
+    
+    
+  def setFunction(self,funcname):
+    cdef CosmoFunction f
+    
+    f = POWER_EFSTATHIOU
+  
+    if funcname=='EFSTATHIOU':
+      f = POWER_EFSTATHIOU
+    elif funcname=='HU_WIGGLES':
+      f = HU_WIGGLES
+    elif funcname=='HU_BARYON':
+      f = HU_BARYON
+    elif funcname=='BARDEEN':
+      f = POWER_BARDEEN
+    elif funcname=='SUGIYAMA':
+      f = POWER_SUGIYAMA
+
+    self.power.setFunction(f)
+    
+  cdef double _compute(self, double k):
+    k *= self.power.h
+    return self.power.power(k)
+    
+  def compute(self, k):
+    cdef np.ndarray out
+    cdef double kval
+    cdef tuple i
+    
+    if isinstance(k, np.ndarray):
+      out = np.empty(k.shape, dtype=np.float64)
+      for i,kval in np.ndenumerate(k):
+        out[i] = self._compute(kval)
+      return out
+    else:
+      return self._compute(k)

python/cosmotool/__init__.py

@@ -1,5 +1,5 @@
 from _cosmotool import *
-from grafic import writeGrafic, writeWhitePhase
+from grafic import writeGrafic, writeWhitePhase, readGrafic
 from borg import read_borg_vol
 
 
@@ -13,9 +13,9 @@ class SimulationBare(PySimulationBase):
       raise TypeError("Simulation object to mirror must be a PySimulationBase")
     s = args[0]
 
-    self.positions = s.getPositions()
+    self.positions = [q.copy() for q in s.getPositions()] if s.getPositions() is not None else None
-    self.velocities = s.getVelocities()
+    self.velocities = [q.copy() for q in s.getVelocities()] if s.getVelocities() is not None else None
-    self.identifiers = s.getIdentifiers()
+    self.identifiers = s.getIdentifiers().copy() if s.getIdentifiers() is not None else None
     self.boxsize = s.getBoxsize()
     self.time = s.getTime()
     self.Hubble = s.getHubble()
@@ -28,7 +28,7 @@ class SimulationBare(PySimulationBase):
     def _safe_merge(a, b):
       if b:
         if a:
-          a = [np.append(q, r) for q,r in zip(a,other.positions)]
+          a = [np.append(q, r) for q,r in zip(a,b)]
         else:
           a = b
       return a
@@ -42,15 +42,15 @@ class SimulationBare(PySimulationBase):
       return a
   
   
-    assert self.time == other.time
+    assert self.time == other.getTime()
-    assert self.Hubble == other.Hubble
+    assert self.Hubble == other.getHubble()
-    assert self.boxsize == other.boxsize
+    assert self.boxsize == other.getBoxsize()
-    assert self.Omega_M == other.Omega_M
+    assert self.Omega_M == other.getOmega_M()
-    assert self.Omega_Lambda == other.Omega_Lambda
+    assert self.Omega_Lambda == other.getOmega_Lambda()
     
-    self.positions = _safe_merge(self.positions, other.positions)
+    self.positions = _safe_merge(self.positions, other.getPositions())
-    self.velocities = _safe_merge(self.velocities, other.velocities)
+    self.velocities = _safe_merge(self.velocities, other.getVelocities())
-    self.identifiers = _safe_merge0(self.idenfiers, other.identifiers)
+    self.identifiers = _safe_merge0(self.identifiers, other.getIdentifiers())
     
   def getPositions(self):
     return self.positions
@@ -106,7 +106,7 @@ def loadRamsesAll(basepath, snapshot_id, **kwargs):
   cpu_id = 0
   output = None
   while True:
-    s = loadRamses("%s/output_%05d" % (basepath,snapshot_id), cpu_id, **kwargs)
+    s = loadRamses("%s/output_%05d" % (basepath,snapshot_id), snapshot_id, cpu_id, **kwargs)
     if s == None:
       break
     if output == None:
@@ -115,3 +115,5 @@ def loadRamsesAll(basepath, snapshot_id, **kwargs):
       output.merge(s)
       
     cpu_id += 1
+
+  return output 

python/cosmotool/grafic.py

@@ -1,6 +1,32 @@
 import struct
 import numpy as np
 
+def readGrafic(filename):
+
+  with file(filename, mode="rb") as f:
+        p = struct.unpack("IIIIffffffffI", f.read(4*11 + 2*4))
+        checkPoint0, Nx, Ny, Nz, delta, _, _, _, scalefac, omega0, omegalambda0, h, checkPoint1 = p
+        if checkPoint0 != checkPoint1 or checkPoint0 != 4*11:
+          raise ValueError("Invalid unformatted access")
+          
+        a = np.empty((Nx,Ny,Nz), dtype=np.float32)
+        
+        BoxSize = delta * Nx * h
+        
+        checkPoint = 4*Ny*Nz
+        for i in xrange(Nx):
+            checkPoint = struct.unpack("I", f.read(4))[0]
+            if checkPoint != 4*Ny*Nz:
+              raise ValueError("Invalid unformatted access")
+              
+            a[i, :, :] = np.fromfile(f, dtype=np.float32, count=Ny*Nz).reshape((Ny, Nz))
+
+            checkPoint = struct.unpack("I", f.read(4))[0]
+            if checkPoint != 4*Ny*Nz:
+              raise ValueError("Invalid unformatted access")
+
+  return a, BoxSize, scalefac, omega0, omegalambda0, h      
+
 def writeGrafic(filename, field, BoxSize, scalefac, **cosmo):
 
     with file(filename, mode="wb") as f:

python_sample/icgen/borgicgen.py

@@ -15,6 +15,40 @@ def gen_posgrid(N, L):
 
     return x.flatten(), y.flatten(), z.flatten()
 
+def bin_power(P, L, bins=20, range=(0,1.)):
+
+  N = P.shape[0]
+  ik = np.fft.fftfreq(N, d=L/N)*2*np.pi
+
+  k = np.sqrt(ik[:,None,None]**2 + ik[None,:,None]**2 + ik[None,None,:(N/2+1)]**2)
+
+  H,b = np.histogram(k, bins=bins, range=range)
+  Hw,b = np.histogram(k, bins=bins, weights=P, range=range)
+  
+  return Hw/H, 0.5*(b[1:]+b[0:bins])
+
+def compute_power_from_borg(input_borg, a_borg, cosmo, bins=10, range=(0,1)):
+  borg_vol = ct.read_borg_vol(input_borg)
+  N = borg_vol.density.shape[0]
+
+  cgrowth = CosmoGrowth(**cosmo)
+  D1 = cgrowth.D(1)
+  D1_0 = D1/cgrowth.D(a_borg)
+
+  density_hat, L = ba.half_pixel_shift(borg_vol)
+
+  return bin_power(D1_0**2*np.abs(density_hat)**2/L**3, L, bins=bins, range=range)
+  
+def compute_ref_power(L, N, cosmo, bins=10, range=(0,1), func='HU_WIGGLES'):
+  ik = np.fft.fftfreq(N, d=L/N)*2*np.pi
+
+  k = np.sqrt(ik[:,None,None]**2 + ik[None,:,None]**2 + ik[None,None,:(N/2+1)]**2)
+  p = ct.CosmologyPower(**cosmo)
+  p.setFunction(func)
+  p.normalize(cosmo['SIGMA8'])
+
+  return bin_power(p.compute(k)*cosmo['h']**3, L, bins=bins, range=range)
+
 def run_generation(input_borg, a_borg, a_ic, **cosmo):
     """ Generate particles and velocities from a BORG snapshot. Returns a tuple of 
     (positions,velocities,N,BoxSize,scale_factor)."""
@@ -26,7 +60,6 @@ def run_generation(input_borg, a_borg, a_ic, **cosmo):
 
     density_hat, L = ba.half_pixel_shift(borg_vol)
 
-
     lpt = LagrangianPerturbation(density_hat, L, fourier=True)
 
     # Generate grid

python_sample/icgen/gen_ic_from_borg.py

@@ -1,11 +1,14 @@
+import cosmotool as ct
 import borgicgen as bic
 
 cosmo={'omega_M_0':0.3175, 'h':0.6711}
 cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
 cosmo['omega_k_0'] = 0
+cosmo['omega_B_0']=0.049
+cosmo['SIGMA8']=0.8344
 
-zstart=50
+zstart=0
 astart=1/(1.+zstart)
 
-if __name__=="__main__":
+#if __name__=="__main__":
-    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo), **cosmo)
+#    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo), **cosmo)

sample/CMakeLists.txt

@@ -88,6 +88,8 @@ if (Boost_FOUND)
   target_link_libraries(simple3DFilter ${tolink})
 endif (Boost_FOUND)
 
-add_executable(gadgetToDensity gadgetToDensity.cpp)
+IF (ENABLE_OPENMP)
-target_link_libraries(gadgetToDensity ${tolink})
+  add_executable(gadgetToDensity gadgetToDensity.cpp)
+  target_link_libraries(gadgetToDensity ${tolink})
+ENDIF (ENABLE_OPENMP)
 

src/CMakeLists.txt

@@ -11,9 +11,12 @@ SET(CosmoTool_SRCS
 	miniargs.cpp
 	growthFactor.cpp
 	cosmopower.cpp
-  cic.cpp
 )
 
+IF (ENABLE_OPENMP)
+  set(CosmoTool_SRCS ${CosmoTool_SRCS} cic.cpp)
+ENDIF (ENABLE_OPENMP)
+
 IF(FOUND_NETCDF3)
 	SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc3.cpp)
 ELSE(FOUND_NETCDF3)

src/cic.cpp

@@ -32,7 +32,7 @@ same conditions as regards security.
 The fact that you are presently reading this means that you have had
 knowledge of the CeCILL license and that you accept its terms.
 +*/
-
+#include <omp.h>
 #include <assert.h>
 #include <math.h>
 #include <inttypes.h>
@@ -60,89 +60,29 @@ void CICFilter::resetMesh()
 
 void CICFilter::putParticles(CICParticles *particles, uint32_t N)
 {
-#if 0
+  int threadUsed = omp_get_max_threads();
-  uint32_t numCorners = 1 << NUMDIMS;
+  double *threadedDensity[threadUsed];
+  bool threadActivated[threadUsed];
+  uint32_t tUsedMin[threadUsed], tUsedMax[threadUsed];
 
-  for (uint32_t i = 0; i < N; i++)
+  for (int t = 0; t < threadUsed; t++)
-    {
+  {
-      Coordinates xyz;
+    threadedDensity[t] = new double[totalSize];
-      int32_t ixyz[NUMDIMS];
+    std::fill(threadedDensity[t], threadedDensity[t]+totalSize, 0);
-      int32_t rxyz[NUMDIMS];
+  }
-      CICType alpha[NUMDIMS];
-      CICType beta[NUMDIMS];
-      for (int j = 0; j < NUMDIMS; j++)
-	{
-	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
-	  ixyz[j] = (int32_t)floor(xyz[j] - 0.5);
-	  beta[j] = xyz[j] - ixyz[j] - 0.5;
-	  alpha[j] = 1 - beta[j];
-	  if (ixyz[j] < 0)
-	    ixyz[j] = szGrid-1;
-	}
 
-      CICType tot_mass = 0;
+  std::fill(threadActivated, threadActivated+threadUsed, false);
-      for (int j = 0; j < numCorners; j++)
+  std::fill(tUsedMin, tUsedMin+threadUsed, totalSize);
-	{
+  std::fill(tUsedMax, tUsedMax+threadUsed, 0);
-	  CICType rel_mass = 1;
-	  uint32_t idx = 0;
-	  uint32_t mul = 1;
-	  uint32_t mul2 = 1;
 
-	  for (int k = 0; k < NUMDIMS; k++)
+#pragma omp parallel
-	    {
+  {
-	      uint32_t ipos = ((j & mul2) != 0);
+    int thisThread = omp_get_thread_num();
-
+    double *dg = threadedDensity[thisThread];
-	      if (ipos == 1)
-		{
-		  rel_mass *= beta[k];
-		}
-	      else
-		{
-		  rel_mass *= alpha[k];
-		}
-
-	      rxyz[k] = ixyz[k] + ipos;
-	      
-	      if (rxyz[k] >= szGrid)
-		idx += (rxyz[k] - szGrid) * mul;
-	      else
-		idx += rxyz[k] * mul;
-
-	      mul2 *= 2;
-	      mul *= szGrid;
-	    }
-
-	  assert(rel_mass > 0);
-	  assert(rel_mass < 1);
-	  assert(idx < totalSize);
-	  densityGrid[idx] += rel_mass * particles[i].mass;
-	  tot_mass += rel_mass;
-	}
-      assert(tot_mass < 1.1);
-      assert(tot_mass > 0.9);
-    }
-#endif
-#if 0  
-  for (uint32_t i = 0; i < N; i++)
-    {
-      Coordinates xyz;
-      int32_t ixyz[NUMDIMS];
-      for (int j = 0; j < NUMDIMS; j++)
-	{
-	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
-	  ixyz[j] = (int32_t)round(xyz[j] - 0.5);
-	  if (ixyz[j] < 0)
-	    ixyz[j] = szGrid-1;
-	  else if (ixyz[j] >= szGrid)
-	    ixyz[j] = 0;
-	}
-
-      uint32_t idx = ixyz[0] + ixyz[1] * szGrid + ixyz[2] * szGrid * szGrid;
-      densityGrid[idx] += particles[i].mass;      
-    }
-
-#endif
 
+    threadActivated[thisThread] = true;
+  
+#pragma omp for schedule(static)
   for (uint32_t i = 0; i < N; i++)
     {
       CICType x, y, z;
@@ -193,44 +133,57 @@ void CICFilter::putParticles(CICParticles *particles, uint32_t N)
 
       // 000
       idx = ix + (iy + iz * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * beta_x * beta_y * beta_z;
       
       // 100
       idx = ix2 + (iy + iz * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * alpha_x * beta_y * beta_z;
       
       // 010
       idx = ix + (iy2 + iz * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * beta_x * alpha_y * beta_z;
       
       // 110
       idx = ix2 + (iy2 + iz * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * alpha_x * alpha_y * beta_z;
 
       // 001
       idx = ix + (iy + iz2 * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * beta_x * beta_y * alpha_z;
 
       // 101
       idx = ix2 + (iy + iz2 * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * alpha_x * beta_y * alpha_z;
 
       // 011
       idx = ix + (iy2 + iz2 * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * beta_x * alpha_y * alpha_z;
 
       // 111
       idx = ix2 + (iy2 + iz2 * szGrid) * szGrid;
-      densityGrid[idx] +=
+      dg[idx] +=
 	mass * alpha_x * alpha_y * alpha_z;
+
+      tUsedMin[thisThread] = std::min(tUsedMin[thisThread], idx);
+      tUsedMax[thisThread] = std::max(tUsedMax[thisThread], idx);
     }
+  }
+
+  for (int t = 0; t < threadUsed; t++)
+  {
+    if (!threadActivated[t])
+      continue;
+#pragma omp parallel for schedule(static)
+    for (long p = tUsedMin[t]; p < tUsedMax[t]; p++)
+      densityGrid[p] += threadedDensity[t][p];
+  }
 }
   
 void CICFilter::getDensityField(CICType*& field, uint32_t& res)
@@ -238,3 +191,4 @@ void CICFilter::getDensityField(CICType*& field, uint32_t& res)
   field = densityGrid;
   res = totalSize;
 }
+