fixed merge

c_tools/mock/flash_loader.cpp

@@ -0,0 +1,111 @@
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadFlash.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class FlashLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  SimuData *gadget_header;
+  string snapshot_name;
+public:
+  FlashLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), gadget_header(header)
+  {
+  }
+  
+  ~FlashLoader()
+  {
+    delete gadget_header;
+  }
+  
+  SimuData *getHeader() {
+    return gadget_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+    if (onefile && id > 0)
+      return 0;
+    if (id >= _num_files)
+      return 0;
+
+    if (onefile)
+      d = loadFlashMulti(snapshot_name.c_str(), -1, load_flags);
+    else
+      d = loadFlashMulti(snapshot_name.c_str(), id, load_flags);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    applyTransformations(d);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *flashLoader(const std::string& snapshot, int flags)
+{
+  bool singleFile;
+  int num_files;
+  SimuData *d;
+
+  try
+    {
+      d = loadFlashMulti(snapshot.c_str(), -1, 0);
+      singleFile = true;
+      num_files = 1;
+    }
+  catch (const NoSuchFileException& e)
+    {
+      try
+        {
+          d = loadFlashMulti(snapshot.c_str(), 0, 0);
+	  num_files = 0;
+        }
+      catch(const NoSuchFileException& e)
+        {
+          return 0;
+        }
+    }
+    
+  assert(d != 0);
+  SimuData *header = d;
+    
+  if (!singleFile)
+    {
+      try
+	{
+	  while ((d = loadFlashMulti(snapshot.c_str(), num_files, 0)) != 0)
+	    {
+	      num_files++;
+	      delete d;
+	    }
+	}
+      catch(const NoSuchFileException& e)
+	{
+	}
+    }
+    
+  return new FlashLoader(snapshot, header, flags, singleFile, num_files);
+}

c_tools/mock/gadget_loader.cpp

@@ -0,0 +1,124 @@
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadGadget.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class GadgetLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  double unitMpc;
+  SimuData *gadget_header;
+  string snapshot_name;
+public:
+  GadgetLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num, double unit)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), unitMpc(1/unit), gadget_header(header)
+  {
+  }
+  
+  ~GadgetLoader()
+  {
+    delete gadget_header;
+  }
+  
+  SimuData *getHeader() {
+    return gadget_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+    if (onefile && id > 0)
+      return 0;
+    if (id >= _num_files)
+      return 0;
+
+    if (onefile)
+      d = loadGadgetMulti(snapshot_name.c_str(), -1, load_flags);
+    else
+      d = loadGadgetMulti(snapshot_name.c_str(), id, load_flags);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    for (int k = 0; k < 3; k++)
+      {
+        if (d->Pos[k] != 0)
+          {
+            for (long i = 0; i < d->NumPart; i++)
+              d->Pos[k][i] *= unitMpc;
+          }
+      }
+    d->BoxSize *= unitMpc;
+
+    applyTransformations(d);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags)
+{
+  bool singleFile;
+  int num_files;
+  SimuData *d;
+
+  try
+    {
+      d = loadGadgetMulti(snapshot.c_str(), -1, 0);
+      singleFile = true;
+      num_files = 1;
+    }
+  catch (const NoSuchFileException& e)
+    {
+      try
+        {
+          d = loadGadgetMulti(snapshot.c_str(), 0, 0);
+	  num_files = 0;
+        }
+      catch(const NoSuchFileException& e)
+        {
+          return 0;
+        }
+    }
+    
+  assert(d != 0);
+  SimuData *header = d;
+
+  header->BoxSize /= Mpc_unitLength;
+    
+  if (!singleFile)
+    {
+      try
+	{
+	  while ((d = loadGadgetMulti(snapshot.c_str(), num_files, 0)) != 0)
+	    {
+	      num_files++;
+	      delete d;
+	    }
+	}
+      catch(const NoSuchFileException& e)
+	{
+	}
+    }
+    
+  return new GadgetLoader(snapshot, header, flags, singleFile, num_files, Mpc_unitLength);
+}

c_tools/mock/generateMock.cpp

@@ -1,214 +1,28 @@
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/format.hpp>
 #include <cmath>
 #include <cassert>
 #include <iostream>
 #include <fstream>
+#include <sstream>
 #include <string>
+#include <algorithm>
 #include <CosmoTool/loadSimu.hpp>
-#include <CosmoTool/loadRamses.hpp>
-#include <CosmoTool/loadGadget.hpp>
-#include <CosmoTool/loadFlash.hpp>
 #include <CosmoTool/interpolate.hpp>
 #include <CosmoTool/fortran.hpp>
 #include "generateMock_conf.h"
 #include "gslIntegrate.hpp"
 #include <netcdfcpp.h>
+#include "simulation_loader.hpp"
 
 using namespace std;
 using namespace CosmoTool;
+using boost::format;
 
 #define LIGHT_SPEED 299792.458
 
-static double gadgetUnit=1e-3;
+typedef boost::function2<void, SimuData*, double*> MetricFunctor;
-
-SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRedshift)
-{
-  SimuData *d, *outd;
-
-  d = loadRamsesSimu(basename, baseid, -1, true, 0);
-  outd = new SimuData;
-
-  outd->NumPart = d->TotalNumPart;
-  outd->BoxSize = d->BoxSize;
-  outd->TotalNumPart = outd->NumPart;
-  outd->Hubble = d->Hubble;
-  outd->Omega_Lambda = d->Omega_Lambda;
-  outd->Omega_M = d->Omega_M;
-  outd->time = d->time;
-
-  for (int k = 0; k < 3; k++)
-    outd->Pos[k] = new float[outd->NumPart];
-  outd->Vel[2] = new float[outd->NumPart];
-  delete d;
-
-  int curCpu = 0;
-  cout << "loading cpu 0 " << endl;
-  while (d = loadRamsesSimu(basename, baseid, curCpu, true, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID))
-    {
-       for (int k = 0; k < 3; k++)
-         for (int i = 0; i < d->NumPart; i++)
-           {
-             assert(d->Id[i] >= 1);
-             assert(d->Id[i] <= outd->TotalNumPart);
-             outd->Pos[k][d->Id[i]-1] = d->Pos[k][i];
-             outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i];
-           }
-
-       if (goRedshift)
-         for (int i = 0; i < d->NumPart; i++)
-            outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.;
-
-       delete d;
-       curCpu++;
-       cout << "loading cpu " << curCpu  << endl;
-    }
-
-  return outd;
-}
-
-
-
-SimuData *myLoadGadget(const char *fname, int id, int flags)
-{
-  SimuData *sim = loadGadgetMulti(fname, id, flags);
-  sim->BoxSize *= gadgetUnit*1000;
-  for (int j = 0; j < 3; j++)
-  {
-    if (sim->Pos[j] != 0) {
-      for (long i = 0; i < sim->NumPart; i++)
-        sim->Pos[j][i] *= gadgetUnit*1000;
-    }
-  }
-  return sim;
-}
-
-SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, SimuData *(*loadFunction)(const char *fname, int id, int flags))
-{
-  SimuData *d, *outd;
-  bool singleFile = false;
-
-  try
-    {
-      d = loadFunction(gadgetname, -1, 0);
-      singleFile = true;
-    }
-  catch (const NoSuchFileException& e)
-    {
-      try
-	{
-	  d = loadFunction(gadgetname, 0, 0);
-	}
-      catch(const NoSuchFileException& e)
-	{
-	  return 0;
-	}
-    }
-  outd = new SimuData;
-
-  outd->NumPart = d->TotalNumPart;
-  outd->BoxSize = d->BoxSize/1000;
-  outd->TotalNumPart = outd->NumPart;
-  outd->Hubble = d->Hubble;
-  outd->Omega_Lambda = d->Omega_Lambda;
-  outd->Omega_M = d->Omega_M;
-  outd->time = d->time;
-
-  for (int k = 0; k < 3; k++)
-    outd->Pos[k] = new float[outd->NumPart];
-  outd->Vel[2] = new float[outd->NumPart];
-  delete d;
-
-  int curCpu = singleFile ? -1 : 0;
-  cout << "loading file 0 " << endl;
-  try
-    {
-      while (1)
-	{
-	  d = loadFunction(gadgetname, curCpu, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID);
-	  for (int k = 0; k < 3; k++)
-	    for (int i = 0; i < d->NumPart; i++)
-	      {
-		assert(d->Id[i] >= 1);
-		assert(d->Id[i] <= outd->TotalNumPart);
-		outd->Pos[k][d->Id[i]-1] = d->Pos[k][i]/1000;
-		outd->Vel[2][d->Id[i]-1] = d->Vel[velAxis][i];
-	      }
-	  
-	  if (goRedshift)
-	    for (int i = 0; i < d->NumPart; i++)
-	      outd->Pos[velAxis][d->Id[i]-1] += d->Vel[velAxis][i]/100.;
-
-	  delete d;
-          if (singleFile)
-            break;
-	  curCpu++;
-	  cout << "loading file " << curCpu  << endl;
-	}
-    }
-  catch (const NoSuchFileException& e)
-    {
-    }
-
-  return outd;
-}
-
-
-SimuData *doLoadMultidark(const char *multidarkname)
-{
-  SimuData *outd;
-  FILE *fp;
-  int actualNumPart;
-
-  outd = new SimuData;
-  cout << "opening multidark file " << multidarkname << endl;
-  fp = fopen(multidarkname, "r");
-  if (fp == NULL) {
-    cout << "could not open file!" << endl;
-    return 0;
-  }
-  fscanf(fp, "%f\n", &outd->BoxSize);
-  fscanf(fp, "%f\n", &outd->Omega_M);
-  fscanf(fp, "%f\n", &outd->Hubble);
-  fscanf(fp, "%f\n", &outd->time);
-  fscanf(fp, "%ld\n", &outd->NumPart);
-
-  outd->time = 1./(1.+outd->time); // convert to scale factor
-  outd->TotalNumPart = outd->NumPart;
-  outd->Omega_Lambda = 1.0 - outd->Omega_M;
-
-  for (int k = 0; k < 3; k++)
-    outd->Pos[k] = new float[outd->NumPart];
-  outd->Vel[2] = new float[outd->NumPart];
-  outd->Id = new int[outd->NumPart];
-  outd->uniqueID = new float[outd->NumPart];
-
-  cout << "loading multidark particles" << endl;
-  actualNumPart = 0;
-	for (int i = 0; i < outd->NumPart; i++) {
-    fscanf(fp, "%d %f %f %f %f\n", &outd->Id[i], 
-                                &outd->Pos[0][i], &outd->Pos[1][i], 
-                                &outd->Pos[2][i], &outd->Vel[2][i]);
-
-    outd->uniqueID[i] = 1.0 * outd->Id[i];
-
-    if (i < 10) printf("TEST %d %d\n", i, outd->Id[i]);
-    if (outd->Id[i] == -99 && 
-        outd->Pos[0][i] == -99 && outd->Pos[1][i] == -99 && 
-        outd->Pos[2][i] == -99 && outd->Vel[2][i] == -99) {
-   printf("FOUND END\n");
-      break;
-    } else {
-      actualNumPart++;
-    }
-  }
-  fclose(fp);
-
-  outd->NumPart = actualNumPart;
-  outd->TotalNumPart = actualNumPart;
-  return outd;
-}
-
-
-static double cubic(double a)
 {
   return a*a*a;
 }
@@ -247,7 +61,7 @@ Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z
   return buildFromVector(pairs);
 }
 
-void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double*& expfact, bool cosmo_flag)
+void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double* expfact, bool cosmo_flag)
 {
   int x0, x1, x2;
 
@@ -275,8 +89,6 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
   TotalExpansion e_computer;
   double baseComovingDistance;
 
-  expfact = new double[data->NumPart];
-
   cout << "Using base redshift z=" << z0 << " " << z0+8*data->BoxSize*100/LIGHT_SPEED << endl;
 
   e_computer.Omega_M = data->Omega_M;
@@ -304,7 +116,8 @@ void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double
           else
             z = reduced_red*LIGHT_SPEED/100.0;
 
-          expfact[i] = z / z_old; 
+          if (expfact)
+            expfact[i] = z / z_old; 
 	  // Add peculiar velocity
 	  if (pecvel)
 	    z += v/100;
@@ -393,39 +206,29 @@ void generateOutput(SimuData *data, int axis,
     }
   f.endCheckpoint();
 
-  cout << "Writing unique ID..." << endl;
+  long *uniqueID = data->as<long>("uniqueID");
-  f.beginCheckpoint();
+  if (uniqueID != 0)
-  for (uint32_t i = 0; i < data->NumPart; i++)
     {
-      f.writeReal32(data->uniqueID[i]);
+      cout << "Writing unique ID..." << endl;
+      f.beginCheckpoint();
+      for (uint32_t i = 0; i < data->NumPart; i++)
+	{
+	  f.writeInt64(uniqueID[i]);
+	}
+      f.endCheckpoint();
     }
-  f.endCheckpoint();
 }
 
-void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info)
+// This function prepares the list of targets for the specified snapshot. The target list is not
+// cleared. That way new particles can be appended if this is a multi-file snapshot.
+void selectBox(SimuData *simu, std::vector<long>& targets, generateMock_info& args_info)
 {
   float subsample = args_info.subsample_given ? args_info.subsample_arg : 1.0;
-  uint32_t goodParticles = 0;
   double ranges[3][2] = {
     { args_info.rangeX_min_arg, args_info.rangeX_max_arg },
     { args_info.rangeY_min_arg, args_info.rangeY_max_arg },
     { args_info.rangeZ_min_arg, args_info.rangeZ_max_arg }
   };
-  double mul[3];
-  float minmax[2][3];
-  int *particle_id;
-  bool *random_acceptance = 0;
-
-  boxed = new SimuData;
-  boxed->Hubble = simu->Hubble;
-  boxed->Omega_M = simu->Omega_M;
-  boxed->Omega_Lambda = simu->Omega_Lambda;
-  boxed->time = simu->time;
-  boxed->BoxSize = simu->BoxSize;
-
-  random_acceptance = new bool[simu->NumPart];
-
-  for (int j = 0; j < 3; j++) minmax[1][j] = minmax[0][j] = simu->Pos[j][0];
 
   for (uint32_t i = 0; i < simu->NumPart; i++)
     {
@@ -435,102 +238,134 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info&
 	acceptance = 
 	  acceptance &&
 	  (simu->Pos[j][i] > ranges[j][0]) && 
-	  (simu->Pos[j][i] < ranges[j][1]);
+	  (simu->Pos[j][i] < ranges[j][1]);	
-        minmax[0][j] = min(simu->Pos[j][i], minmax[0][j]);
-        minmax[1][j] = max(simu->Pos[j][i], minmax[1][j]);
       }
-      random_acceptance[i] = acceptance && (drand48() <= subsample);
+      
-      if (random_acceptance[i])
+      if (acceptance && (drand48() <= subsample))
-	goodParticles++;
+	targets.push_back(i);
     }
+}
 
-  cout << "Subsample fraction: " << subsample << endl;
+void createBox(SimuData *simu, vector<long>& targets, vector<long>& snapshot_split, SimuData *& boxed, generateMock_info& args_info)
-  cout << "Min range = " << ranges[0][0] << " " << ranges[1][0] << " " << ranges[2][0] << endl;
+{
-  cout << "Max range = " << ranges[0][1] << " " << ranges[1][1] << " " << ranges[2][1] << endl;
+  double *ranges = new double[6];
+  double *mul = new double[3];
+  long *simu_uniqueID = simu->as<long>("uniqueID");
 
-  cout << "Min position = " << minmax[0][0] << " " << minmax[0][1] << " " << minmax[0][2] << endl;
+  ranges[0] = args_info.rangeX_min_arg;
-  cout << "Max position = " << minmax[1][0] << " " << minmax[1][1] << " " << minmax[1][2] << endl;
+  ranges[1] = args_info.rangeX_max_arg;
- 
+  ranges[2] = args_info.rangeY_min_arg;
-  cout << "Number of accepted particles: " << goodParticles << endl; 
+  ranges[3] = args_info.rangeY_max_arg;
+  ranges[4] = args_info.rangeZ_min_arg;
+  ranges[5] = args_info.rangeZ_max_arg;
+
+  boxed = new SimuData;
+  boxed->Hubble = simu->Hubble;
+  boxed->Omega_M = simu->Omega_M;
+  boxed->Omega_Lambda = simu->Omega_Lambda;
+  boxed->time = simu->time;
+  boxed->BoxSize = simu->BoxSize;
+  boxed->NumPart = targets.size();
 
   for (int j = 0; j < 3; j++)
     {
-      boxed->Pos[j] = new float[goodParticles];
+      boxed->Pos[j] = new float[boxed->NumPart];
       boxed->Vel[j] = 0;
-      mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
+      mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]);
     }
-  boxed->uniqueID = new float[goodParticles];
+  cout << "Min range = " << ranges[0] << " " << ranges[2] << " " << ranges[4] << endl;
-
+  cout << "Max range = " << ranges[1] << " " << ranges[3] << " " << ranges[5] << endl;
+  cout << "Number of accepted particles: " << boxed->NumPart << endl; 
   cout << "Rescaling factors = " << mul[0] << " " << mul[1] << " " << mul[2] << endl;
-  boxed->NumPart = goodParticles;
 
-  particle_id = new int[goodParticles];
+  long *uniqueID = new long[boxed->NumPart];
-  double *expansion_fac = new double[goodParticles];
+  long *particle_id = new long[boxed->NumPart];
+  double *expansion_fac = new double[boxed->NumPart];
+  long *snap_split = new long[snapshot_split.size()];
+  int *numsnap_info = new int[1];
 
+  copy(targets.begin(), targets.end(), particle_id);
+  copy(snapshot_split.begin(), snapshot_split.end(), snap_split);
+  *numsnap_info = snapshot_split.size();
+
+  boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
+  boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
+  boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>); 
+  boxed->new_attribute("mul", mul, delete_adaptor<double>);
+  boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
+  boxed->new_attribute("snapshot_split", snap_split, delete_adaptor<long>);
+  boxed->new_attribute("num_snapshots", numsnap_info, delete_adaptor<int>);
+}
+
+void buildBox(SimuData *simu, long num_targets, long loaded, 
+	      SimuData *boxed, double *efac)
+{
   uint32_t k = 0;
-  for (uint32_t i = 0; i < simu->NumPart; i++)
+  long *uniqueID = boxed->as<long>("uniqueID");
+  long *simu_uniqueID = simu->as<long>("uniqueID");
+  double *expansion_fac = boxed->as<double>("expansion_fac");
+  double *mul = boxed->as<double>("mul");
+  double *ranges = boxed->as<double>("ranges");
+  long *particle_id = boxed->as<long>("particle_id");
+
+  for (uint32_t i = 0; i < num_targets; i++, loaded++)
     {
-      bool acceptance = random_acceptance[i]; 
+      long pid = particle_id[loaded];
       
-      if (acceptance)
+      for (int j = 0; j < 3; j++)
 	{
-	  for (int j = 0; j < 3; j++)
+	  boxed->Pos[j][loaded] = (simu->Pos[j][pid]-ranges[j*2])*mul[j];
-	    {
+	  assert(boxed->Pos[j][loaded] > 0);
-	      boxed->Pos[j][k] = (simu->Pos[j][i]-ranges[j][0])*mul[j];
+	  assert(boxed->Pos[j][loaded] < 1);
-	      assert(boxed->Pos[j][k] > 0);
-	      assert(boxed->Pos[j][k] < 1);
-	    }
-    boxed->uniqueID[k] = simu->uniqueID[i];
-
-	  particle_id[k] = i;
-          expansion_fac[k] = efac[i];
-	  k++;
 	}
+      uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[pid] : 0;
+      expansion_fac[loaded] = efac[pid];
     }
+}
 
-  delete[] random_acceptance;
+void saveBox(SimuData *&boxed, const std::string& outbox)
+{
+  double *ranges = boxed->as<double>("ranges");
+  NcFile f(outbox.c_str(), NcFile::Replace);
+  long *particle_id = boxed->as<long>("particle_id");
+  double *expansion_fac = boxed->as<double>("expansion_fac");
+  long *snapshot_split = boxed->as<long>("snapshot_split");
+  int num_snapshots = *boxed->as<int>("num_snapshots");
+  long *uniqueID = boxed->as<long>("uniqueID");
 
-  NcFile f(args_info.outputParameter_arg, NcFile::Replace);
+  f.add_att("range_x_min", ranges[0]);
-
+  f.add_att("range_x_max", ranges[1]);
-  f.add_att("range_x_min", ranges[0][0]);
+  f.add_att("range_y_min", ranges[2]);
-  f.add_att("range_x_max", ranges[0][1]);
+  f.add_att("range_y_max", ranges[3]);
-  f.add_att("range_y_min", ranges[1][0]);
+  f.add_att("range_z_min", ranges[4]);
-  f.add_att("range_y_max", ranges[1][1]);
+  f.add_att("range_z_max", ranges[5]);
-  f.add_att("range_z_min", ranges[2][0]);
-  f.add_att("range_z_max", ranges[2][1]);
   f.add_att("mask_index", -1);
 
   NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
-  NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim);
+  NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots);
+  NcVar *v = f.add_var("particle_ids", ncLong, NumPart_dim);
   NcVar *v2 = f.add_var("expansion", ncDouble, NumPart_dim);
+  NcVar *v3 = f.add_var("snapshot_split", ncLong, NumSnap_dim);
 
   v->put(particle_id, boxed->NumPart);
-  v2->put(expansion_fac, boxed->NumPart);
+  v2->put(expansion_fac, boxed->NumPart);  
-
+  v3->put(snapshot_split, num_snapshots);
-  delete[] particle_id;
+  if (uniqueID != 0)
-  delete[] expansion_fac;
+    {
-
+      NcVar *v4 = f.add_var("unique_ids", ncLong, NumPart_dim);
-
+      v4->put(uniqueID, boxed->NumPart);
-/*
+    }
-  FILE *fp = fopen("sample_info.txt", "w");
-  fprintf(fp, "x_min = %f\n", ranges[0][0]);
-  fprintf(fp, "x_max = %f\n", ranges[0][1]);
-  fprintf(fp, "y_min = %f\n", ranges[1][0]);
-  fprintf(fp, "y_max = %f\n", ranges[1][1]);
-  fprintf(fp, "z_min = %f\n", ranges[2][0]);
-  fprintf(fp, "z_max = %f\n", ranges[2][1]);
-  fprintf(fp, "mask_index = -1\n");
-  fprintf(fp, "total_particles = %d\n", boxed->NumPart);
-  fclose(fp);
-*/
 }
 
-void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info)
+void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& args_info)
 {
   NcFile f(args_info.inputParameter_arg);
   NcVar *v;
-  int *particle_id;
+  long *particle_id;
   double *expansion_fac;
+  long *uniqueID;
+  long *snapshot_split;
+  int *num_snapshots = new int[1];
 
   boxed = new SimuData;
   boxed->Hubble = simu->Hubble;
@@ -540,46 +375,79 @@ void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, genera
   boxed->BoxSize = simu->BoxSize;
 
   NcVar *v_id = f.get_var("particle_ids");
+  NcVar *v_snap = f.get_var("snapshot_split");
   long *edges1;
-  double ranges[3][2];
+  long *dim_snap;
-  double mul[3];
+  double *ranges;
+  double *mul;
 
   edges1 = v_id->edges();
+  dim_snap = v_snap->edges();
   assert(v_id->num_dims()==1);
+  assert(v_snap->num_dims()==1);
 
   boxed->NumPart = edges1[0];
+  *num_snapshots = dim_snap[0];
+  delete[] dim_snap;
   delete[] edges1;
 
-  particle_id = new int[boxed->NumPart];
+  particle_id = new long[boxed->NumPart];
+  uniqueID = new long[boxed->NumPart];
+  mul = new double[3];
+  ranges = new double[6];
+  snapshot_split = new long[*num_snapshots];
+
+
+  boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+  boxed->new_attribute("mul", mul, delete_adaptor<double>);
+  boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
+  boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
+  boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor<int>);
+  boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor<long>);
 
   v_id->get(particle_id, boxed->NumPart);
+  v_snap->get(snapshot_split, *num_snapshots);
 
-  ranges[0][0] = f.get_att("range_x_min")->as_double(0);
+  ranges[0] = f.get_att("range_x_min")->as_double(0);
-  ranges[0][1] = f.get_att("range_x_max")->as_double(0);
+  ranges[1] = f.get_att("range_x_max")->as_double(0);
-  ranges[1][0] = f.get_att("range_y_min")->as_double(0);
+  ranges[2] = f.get_att("range_y_min")->as_double(0);
-  ranges[1][1] = f.get_att("range_y_max")->as_double(0);
+  ranges[3] = f.get_att("range_y_max")->as_double(0);
-  ranges[2][0] = f.get_att("range_z_min")->as_double(0);
+  ranges[4] = f.get_att("range_z_min")->as_double(0);
-  ranges[2][1] = f.get_att("range_z_max")->as_double(0);
+  ranges[5] = f.get_att("range_z_max")->as_double(0);
 
   for (int j = 0; j < 3; j++)
     {
       boxed->Pos[j] = new float[boxed->NumPart];
       boxed->Vel[j] = 0;
-      mul[j] = 1.0/(ranges[j][1] - ranges[j][0]);
+      mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]);
     }
   
   uint32_t k = 0;
-  for (uint32_t i = 0; i < boxed->NumPart; i++)
+  NcVar *v_uniq = f.get_var("unique_ids");
+  v_uniq->get(uniqueID, boxed->NumPart);
+}
+
+void makeBoxFromSimulation(SimulationLoader *loader, SimuData* &boxed, MetricFunctor metric, generateMock_info& args_info)
+{
+  vector<long> targets, split;
+  long previous_target_num = 0;
+
+  for (int nf = 0; nf < loader->num_files(); nf++)
     {
-      int id = particle_id[i];
+      SimuData *simu;
+      double *expfact;
 
-      for (int j = 0; j < 3; j++)
+      cout << format("Analyzing and selecting targets in file number %d / %d") % (nf+1) % loader->num_files() << endl;
-	{
+      simu = loader->loadFile(nf);      
-	  boxed->Pos[j][i] = (simu->Pos[j][id]-ranges[j][0])*mul[j];
-	}
-    } 
 
-  delete[] particle_id;
+      metric(simu, 0);
+
+      selectBox(simu, targets, args_info);
+      split.push_back(targets.size() - previous_target_num);
+      previous_target_num = targets.size();
+    }
+
+  createBox(loader->getHeader(), targets, split, boxed, args_info);
 }
 
 int main(int argc, char **argv)
@@ -587,6 +455,7 @@ int main(int argc, char **argv)
   generateMock_info args_info;
   generateMock_conf_params args_params;
   SimuData *simu, *simuOut;
+  SimulationLoader *loader;
  
   generateMock_conf_init(&args_info);
   generateMock_conf_params_init(&args_params);
@@ -612,78 +481,92 @@ int main(int argc, char **argv)
   
   generateMock_conf_print_version();
 
-  gadgetUnit=args_info.gadgetUnit_arg;
-
   if (args_info.ramsesBase_given || args_info.ramsesId_given)
     {
       if (args_info.ramsesBase_given && args_info.ramsesId_given)   {
-	simu = doLoadRamses(args_info.ramsesBase_arg, 
+	loader = ramsesLoader(args_info.ramsesBase_arg, 
-			    args_info.ramsesId_arg, 
+			      args_info.ramsesId_arg, 
-			    args_info.axis_arg, false);
+			      false,
+			      NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID);
       }
       else
 	{
 	  cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl;
 	  return 1;
 	}
-
-      if (simu == 0)
-	{
-	  cerr << "Error while loading" << endl;
-	  return 1;
-	}
     }
-  else if (args_info.gadget_given || args_info.flash_given || args_info.multidark_given)
+  else if (args_info.gadget_given)
     {
-      if (args_info.gadget_given && args_info.flash_given)
+      loader = gadgetLoader(args_info.gadget_arg, 1/args_info.gadgetUnit_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID);
-	{
+    }
-	  cerr << "Do not know which file to use: Gadget or Flash ?" << endl;
+  else if (args_info.flash_given)
-	  return 1;
+    {
-	}
+      loader = flashLoader(args_info.flash_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID); 
-
+  }
-      if (args_info.multidark_given) {
+  else if (args_info.multidark_given)
-        simu = doLoadMultidark(args_info.multidark_arg);      
+    {
-      }
+      loader = multidarkLoader(args_info.multidark_arg);      
-
-      if (args_info.gadget_given) {
-	simu = doLoadSimulation(args_info.gadget_arg, args_info.axis_arg, false, myLoadGadget);      
-      }
-      if (args_info.flash_given) {
-	simu = doLoadSimulation(args_info.flash_arg, args_info.axis_arg, false, loadFlashMulti); 
-      }
-
-      if (simu == 0)
-	{
-	  cerr << "Error while loading " << endl;
-	  return 1;
-	}
     }
   else
     {
-      cerr << "Either a ramses snapshot or a gadget snapshot is required." << endl;
+      cerr << "A simulation snapshot is required to generate a mock catalog." << endl;
       return 1;
     }
-  cout << "Hubble = " << simu->Hubble << endl;
+     
-  cout << "Boxsize = " << simu->BoxSize << endl;
+  if (loader == 0)
-  cout << "Omega_M = " << simu->Omega_M << endl;
+    {
-  cout << "Omega_Lambda = " << simu->Omega_Lambda << endl;
+      cerr << "Error while loading " << endl;
+      return 1;
+    }
+  simu = loader->getHeader();
 
+  {
+    SimuData *header = loader->getHeader();
+    cout << "Hubble = " << header->Hubble << endl;
+    cout << "Boxsize = " << header->BoxSize << endl;
+    cout << "Omega_M = " << header->Omega_M << endl;
+    cout << "Omega_Lambda = " << header->Omega_Lambda << endl;
+    cout << "Subsample fraction: " << (args_info.subsample_given ? 1 : args_info.subsample_arg) << endl;
+  }
   double *expfact;
 
-  metricTransform(simu, args_info.axis_arg, args_info.preReShift_flag,
+  boost::function2<void, SimuData*, double*> metricOperation=
-		  args_info.peculiarVelocities_flag, expfact, 
+    boost::bind(metricTransform, _1, args_info.axis_arg, args_info.preReShift_flag,
-		  args_info.cosmo_flag);
+                      args_info.peculiarVelocities_flag, _2, 
+                      args_info.cosmo_flag);
 
   if (args_info.inputParameter_given)
-    makeBoxFromParameter(simu, expfact, simuOut, args_info);
+    makeBoxFromParameter(loader->getHeader(), simuOut, args_info);
   else
-    makeBox(simu, expfact, simuOut, args_info);
+    makeBoxFromSimulation(loader, simuOut, metricOperation, args_info);
 
-  delete simu;
+  long loaded = 0;
+  for (int nf = 0; nf < loader->num_files(); nf++)
+    {
+      long num_targets = simuOut->as<long>("snapshot_split")[nf];
+      cout << format("Building box from particles in %d / %d") % (nf+1) % loader->num_files() << endl;
+     
+      if (num_targets == 0)
+        {
+          cout << "No particles selected there. Skipping." << endl;
+          continue;
+        }
+      
+      SimuData *simu = loader->loadFile(nf);
+      double *efac = new double[simu->NumPart];
+      metricOperation(simu, efac);
+      buildBox(simu, num_targets, loaded, simuOut, efac);
+      
+      loaded += num_targets;
+      assert(loaded <= simuOut->NumPart);
 
-  generateOutput(simuOut, args_info.axis_arg, args_info.output_arg);
+      delete[] efac;
+    }
 
-  delete simuOut;
+  saveBox(simuOut, args_info.outputParameter_arg);
+  generateOutput(simuOut, args_info.axis_arg, 
+                 args_info.output_arg);
+  
  
   printf("Done!\n"); 
   return 0;

c_tools/mock/multidark_loader.cpp

@@ -0,0 +1,108 @@
+#include <cassert>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class MultiDarkLoader: public SimulationLoader
+{
+protected:
+  SimuData *header;
+  string darkname;
+public:
+  MultiDarkLoader(const std::string& name, SimuData *h)
+    : darkname(name), header(h)
+  {
+  }
+
+  ~MultiDarkLoader()
+  {
+    delete header;
+  }
+
+  int num_files()
+  {
+    return 1;
+  }
+
+  SimuData *getHeader()
+  {
+    return header;
+  }
+
+  SimuData *loadFile(int id)
+  {
+    if (id != 0)
+      return 0;
+
+    ifstream fp(darkname.c_str());
+    SimuData *simu;
+
+    fp >> simu->BoxSize >> simu->Omega_M >> simu->Hubble >> simu->time >> simu->NumPart;
+    simu->time = 1./(1.+simu->time); // convert to scale factor
+    simu->TotalNumPart = simu->NumPart;
+    simu->Omega_Lambda = 1.0 - simu->Omega_M;
+
+    for (int k = 0; k < 3; k++)
+      simu->Pos[k] = new float[simu->NumPart];
+    simu->Vel[2] = new float[simu->NumPart];
+    simu->Id = new int[simu->NumPart];
+    long *uniqueID = new long[simu->NumPart];
+
+    simu->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+
+    cout << "loading multidark particles" << endl;
+    long actualNumPart = 0;
+    for (int i = 0; i < simu->NumPart; i++) {
+
+      fp >> simu->Id[i] >> simu->Pos[0][i] >> simu->Pos[1][i]
+	 >> simu->Pos[2][i] >> simu->Vel[2][i];
+
+      uniqueID[i] = 1;
+      //uniqueID[i] = 1 * simu->Id[i];
+
+      if (simu->Id[i] == -99 && 
+          simu->Pos[0][i] == -99 && simu->Pos[1][i] == -99 && 
+          simu->Pos[2][i] == -99 && simu->Vel[2][i] == -99) {
+        break;
+      } else {
+        actualNumPart++;
+      }
+    }
+
+    simu->NumPart = actualNumPart;
+    simu->TotalNumPart = actualNumPart;
+    return simu;
+  }
+};
+
+SimulationLoader *multidarkLoader(const string& multidarkname)
+{
+  SimuData *header;
+  int actualNumPart;
+  ifstream fp(multidarkname.c_str());
+
+  cout << "opening multidark file " << multidarkname << endl;
+  if (!fp)
+    {
+      cout << "could not open file!" << endl;
+      return 0;
+    }
+
+  header = new SimuData();
+  fp >> header->BoxSize >> header->Omega_M >> header->Hubble >> header->time >> header->NumPart;
+  
+  header->time = 1./(1.+header->time); // convert to scale factor
+  header->TotalNumPart = header->NumPart;
+  header->Omega_Lambda = 1.0 - header->Omega_M;
+
+  return new MultiDarkLoader(multidarkname, header);
+}
+
+  
+

c_tools/mock/ramses_loader.cpp

@@ -0,0 +1,81 @@
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class RamsesLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  int _num_files;
+  int baseid;
+  bool double_precision;
+  SimuData *ramses_header;
+  string snapshot_name;
+public:
+  RamsesLoader(const string& basename, int baseid, bool dp, SimuData *header, int flags, int _num)
+    : snapshot_name(basename), load_flags(flags), _num_files(_num), double_precision(dp),
+      ramses_header(header)
+  {
+  }
+  
+  ~RamsesLoader()
+  {
+    delete ramses_header;
+  }
+  
+  SimuData *getHeader() {
+    return ramses_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+    if (id >= _num_files)
+      return 0;
+
+    d = loadRamsesSimu(snapshot_name.c_str(), baseid, id, double_precision, load_flags);
+    assert(d != 0);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    applyTransformations(d);
+
+    return d;
+  }
+};
+
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags)
+{
+  SimuData *d, *header;
+  int num_files = 0;
+
+  header = loadRamsesSimu(snapshot.c_str(), baseid, 0, double_precision, 0);
+  if (header == 0)
+    return 0;
+
+  while ((d = loadRamsesSimu(snapshot.c_str(), baseid, num_files, double_precision, 0)) != 0)
+    {
+      num_files++;
+      delete d;
+    }
+
+  return new RamsesLoader(snapshot, baseid, double_precision, header, flags, num_files);
+}
+

c_tools/mock/simulation_loader.cpp

@@ -0,0 +1,15 @@
+#include <CosmoTool/loadSimu.hpp>
+#include "simulation_loader.hpp"
+
+using namespace CosmoTool;
+
+void SimulationLoader::applyTransformations(SimuData *s)
+{
+  float redshift_gravity = do_redshift ? 1.0 : 0.0;
+
+  for (int i = 0; i < s->NumPart; i++)
+    {
+      s->Pos[redshift_axis][i] += 
+	redshift_gravity*s->Vel[redshift_axis][i]/100.;
+    }
+}

c_tools/mock/simulation_loader.hpp

@@ -0,0 +1,49 @@
+#ifndef _MOCK_SIMULATION_LOADER_HPP
+#define _MOCK_SIMULATION_LOADER_HPP
+
+#include <string>
+#include <CosmoTool/loadSimu.hpp>
+
+class SimulationLoader
+{
+protected:
+  bool do_redshift;
+  int redshift_axis;
+
+  SimulationLoader() 
+  {
+    do_redshift = false;
+    redshift_axis = 2;
+  }
+ 
+  void applyTransformations(CosmoTool::SimuData *s);
+
+public:
+  virtual ~SimulationLoader() {}
+  
+  void doRedshift(bool set = true) { do_redshift = set; }
+  void setVelAxis(int axis) { redshift_axis = axis; }
+
+  virtual CosmoTool::SimuData *getHeader() = 0;
+  virtual int num_files() = 0;
+  virtual CosmoTool::SimuData* loadFile(int id) = 0;
+
+};
+
+template<typename T>
+void delete_adaptor(void *ptr)
+{
+  T *ptr_T = reinterpret_cast<T *>(ptr);
+
+  delete[] ptr_T;
+}
+
+
+// Unit length is the size of one Mpc in the simulation units
+SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags);
+SimulationLoader *flashLoader(const std::string& snapshot, int flags);
+SimulationLoader *multidarkLoader(const std::string& snapshot);
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags);
+
+
+#endif

external/cosmotool/src/loadGadget.cpp

@@ -10,61 +10,41 @@
 using namespace CosmoTool;
 using namespace std;
 
-PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
-{
-  PurePositionData *data;
-  int p, n;
-  UnformattedRead f(fname);
-  GadgetHeader h;
-  
-  data = new PurePositionData;
-  f.beginCheckpoint();
-  for (int i = 0; i < 6; i++)
-    h.npart[i] = f.readInt32();
-  for (int i = 0; i < 6; i++)    
-    h.mass[i] = f.readReal64();
-  h.time = f.readReal64();
-  h.redshift = f.readReal64();
-  h.flag_sfr = f.readInt32();
-  h.flag_feedback = f.readInt32();
-  for (int i = 0; i < 6; i++)
-    h.npartTotal[i] = f.readInt32();
-  h.flag_cooling = f.readInt32();
-  h.num_files = f.readInt32();
-  data->BoxSize = h.BoxSize = f.readReal64();
-  h.Omega0 = f.readReal64();
-  h.OmegaLambda = f.readReal64();
-  h.HubbleParam = f.readReal64();
-  f.endCheckpoint(true);
-  
-  data->NumPart = 0;
-  for(int k=0; k<5; k++)
-    data->NumPart += h.npart[k];
-  
-  data->pos = new FCoordinates[data->NumPart];
-  
-  f.beginCheckpoint();
-  for(int k = 0, p = 0; k < 5; k++) {
-    for(int n = 0; n < h.npart[k]; n++) {
-      data->pos[p][0] = f.readReal32();
-      data->pos[p][1] = f.readReal32();
-      data->pos[p][2] = f.readReal32();
-      p++;
-    }
-  }
-  f.endCheckpoint();
-  
-  // Skip velocities
-  f.skip((long)data->NumPart*3+2*4);
-  // Skip ids
 
-  return data;
+void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int id)
+{
+  f->beginCheckpoint();
+  for (int i = 0; i < 6; i++)
+    h.npart[i] = f->readInt32();
+  for (int i = 0; i < 6; i++)    
+    h.mass[i] = f->readReal64();
+  data->time = h.time = f->readReal64();
+  h.redshift = f->readReal64();
+  h.flag_sfr = f->readInt32();
+  h.flag_feedback = f->readInt32();
+  for (int i = 0; i < 6; i++)
+    h.npartTotal[i] = f->readInt32();
+  h.flag_cooling = f->readInt32();
+  h.num_files = f->readInt32();
+  data->BoxSize = h.BoxSize = f->readReal64();
+  data->Omega_M = h.Omega0 = f->readReal64();
+  data->Omega_Lambda = h.OmegaLambda = f->readReal64();
+  data->Hubble = h.HubbleParam = f->readReal64();
+  f->endCheckpoint(true);
+  
+  long NumPart = 0, NumPartTotal = 0;
+  for(int k=0; k<6; k++)
+    {
+      NumPart += h.npart[k];
+      NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
+    }
+  data->NumPart = NumPart;
+  data->TotalNumPart = NumPartTotal;
 }
 
-
+SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
-
+				     int loadflags, int GadgetFormat, 
-
+				     SimuFilter filter)
-SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat)
 {
   SimuData *data;
   int p, n;
@@ -98,35 +78,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
   }
 
   long NumPart = 0, NumPartTotal = 0;
-
+  
   try
     {
-      f->beginCheckpoint();
+      loadGadgetHeader(f, h, data, id);
-      for (int i = 0; i < 6; i++)
+
-	h.npart[i] = f->readInt32();
-      for (int i = 0; i < 6; i++)    
-	h.mass[i] = f->readReal64();
-      data->time = h.time = f->readReal64();
-      h.redshift = f->readReal64();
-      h.flag_sfr = f->readInt32();
-      h.flag_feedback = f->readInt32();
-      for (int i = 0; i < 6; i++)
-	h.npartTotal[i] = f->readInt32();
-      h.flag_cooling = f->readInt32();
-      h.num_files = f->readInt32();
-      data->BoxSize = h.BoxSize = f->readReal64();
-      data->Omega_M = h.Omega0 = f->readReal64();
-      data->Omega_Lambda = h.OmegaLambda = f->readReal64();
-      data->Hubble = h.HubbleParam = f->readReal64();
-      f->endCheckpoint(true);
-      
-      for(int k=0; k<6; k++)
-	{
-	  NumPart += h.npart[k];
-	  NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
-	}
-      data->NumPart = NumPart;
-      data->TotalNumPart = NumPartTotal;
       if (GadgetFormat == 1)
         velmul = sqrt(h.time);
       else if (GadgetFormat == 2)
@@ -135,6 +91,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
         cerr << "unknown gadget format" << endl;
 	abort();
       }
+
+      NumPart = data->NumPart;
+      NumPartTotal = data->TotalNumPart;
     }
   catch (const InvalidUnformattedAccess& e)
     {
@@ -256,7 +215,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
       }
     catch (const InvalidUnformattedAccess& e)
       {
-	cerr << "Invalid formatted while reading ID" << endl;
+	cerr << "Invalid unformatted access while reading ID" << endl;
 	delete f;
 	delete data;
 	return 0;

external/cosmotool/src/loadGadget.hpp

@@ -5,10 +5,9 @@
 #include "loadSimu.hpp"
 
 namespace CosmoTool {
-  
-  PurePositionData *loadGadgetPosition(const char *fname);
 
-  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1);
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, 
+			    int GadgetFormat = 1, SimuFilter filter = 0);
  
   // Only single snapshot supported
   void writeGadget(const char *fname, SimuData *data, int GadgetFormat = 1);

external/cosmotool/src/loadSimu.hpp

@@ -1,6 +1,8 @@
 #ifndef __COSMOTOOLBOX_HPP
 #define __COSMOTOOLBOX_HPP
 
+#include <map>
+#include <string>
 
 namespace CosmoTool
 {
@@ -9,9 +11,24 @@ namespace CosmoTool
   static const int NEED_VELOCITY = 4;
   static const int NEED_TYPE = 8;
 
+  struct SimuParticle
+  {
+    float Pos[3];
+    float Vel[3];
+    int type;
+    int id;
+
+    bool flag_vel, flag_type, flag_id;
+  };
+
+  typedef bool (*SimuFilter)(const SimuParticle& p);
+
   class SimuData 
   {
   public:
+    typedef void (*FreeFunction)(void *);
+    typedef std::map<std::string, std::pair<void *, FreeFunction> > AttributeMap;
+
     float BoxSize;
     float time;
     float Hubble;
@@ -24,10 +41,12 @@ namespace CosmoTool
     int *Id;
     float *Pos[3];
     float *Vel[3];
-    float *uniqueID;
     int *type;
+
+    AttributeMap attributes;
+    
   public:
-    SimuData() : Id(0),NumPart(0),type(0), uniqueID(0)  { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; uniqueID=0;}
+    SimuData() : Id(0),NumPart(0),type(0)  { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
     ~SimuData() 
     {
       for (int j = 0; j < 3; j++)
@@ -41,9 +60,37 @@ namespace CosmoTool
 	delete[] type;
       if (Id)
 	delete[] Id;
-      if (uniqueID)
+
-  delete[] uniqueID;
+      for (AttributeMap::iterator i = attributes.begin();
+	   i != attributes.end();
+	   ++i)
+	{
+	  if (i->second.second)
+	    i->second.second(i->second.first);
+	}
     }    
+    
+    template<typename T>
+    T *as(const std::string& n) 
+    {
+      AttributeMap::iterator i = attributes.find(n);
+      if (i == attributes.end())
+	return 0;
+      
+      return reinterpret_cast<T *>(i->second.first);
+    }
+
+    void new_attribute(const std::string& n, void *p, FreeFunction free_func)
+    {
+      AttributeMap::iterator i = attributes.find(n);
+      if (i != attributes.end())
+	{
+	  if (i->second.second)
+	    i->second.second(i->second.first);
+	}
+      attributes[n] = std::make_pair(p, free_func);
+    }
+    
   };
 
 };

external/external_build.cmake

@@ -26,7 +26,7 @@ IF(INTERNAL_GENGETOPT)
 ENDIF(INTERNAL_GENGETOPT)
 
 IF(INTERNAL_HDF5)
-  SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/current/src/hdf5-1.8.9.tar.gz" CACHE STRING "URL to download HDF5 from")
+  SET(HDF5_URL "http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.9/src/hdf5-1.8.9.tar.gz" CACHE STRING "URL to download HDF5 from")
   mark_as_advanced(HDF5_URL)
 ENDIF(INTERNAL_HDF5)
 

pipeline/generateCatalog.py

@@ -9,8 +9,8 @@ import pickle
 
 # ------------------------------------------------------------------------------
 
-if (len(sys.argv) == 0):
+if (len(sys.argv) == 1):
-  print "Usage: ./anylyzeVoids.py parameter_file.py"
+  print "Usage: ./analyzeVoids.py parameter_file.py"
   exit(-1)
 
 if (len(sys.argv) > 1):

pipeline/prepareCatalogs.in.py

@@ -81,8 +81,8 @@ endCatalogStage   = 4
 startAPStage = 1
 endAPStage = 7
 
-ZOBOV_PATH = os.getenv("PWD")+"/../zobov/"
+ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/"
-CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/"
+CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/"
 freshStack = True
 errorBars = "CALCULATED"
 numIncoherentRuns = 100

python_tools/void_python_tools/backend/launchers.py

@@ -60,7 +60,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
     file(parmFile, mode="w").write(conf)
 
     if not (continueRun and jobSuccessful(logFile, "Done!\n")):
-      cmd = "%s --configFile=%s >& %s" % (binPath,parmFile,logFile)
+      cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
         print "done"
@@ -148,7 +148,7 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
     file(parmFile, mode="w").write(conf)
 
     if not (continueRun and jobSuccessful(logFile, "Done!\n")):
-      cmd = "%s --configFile=%s >& %s" % (binPath,parmFile,logFile)
+      cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile)
       os.system(cmd)
       if jobSuccessful(logFile, "Done!\n"):
         print "done"
@@ -199,7 +199,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
     if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK):
       os.unlink(zobovDir+"/voidDesc_"+sampleName+".out")
 
-    cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s >& %s" % \
+    cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s &> %s" % \
           (binPath, datafile,  numZobovDivisions, \
                       "_"+sampleName, numZobovThreads, \
                       binPath, zobovDir, maskIndex, logFile)
@@ -295,7 +295,7 @@ def launchPrune(sample, binPath, thisDataPortion=None,
     cmd += " --outDistances=" + zobovDir+"/boundaryDistances_"+\
                           str(thisDataPortion)+"_"+\
                           str(sampleName)+".out"
-    cmd += " >& " + logFile
+    cmd += " &> " + logFile
     os.system(cmd)
 
     if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \
@@ -423,7 +423,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
   fp.close()
 
   if not (continueRun and jobSuccessful(logFile, "Done!\n")):
-    cmd = "%s --configFile=%s >& %s" % \
+    cmd = "%s --configFile=%s &> %s" % \
           (binPath, parmFile, logFile)
     os.system(cmd)