renamed c_tools to c_src and moved zobov to that directory

c_src/prep/loaders/CMakeLists.txt

@@ -0,0 +1,11 @@
+SET(SIMU_LOADERS_SRCS 
+   simulation_loader.cpp ramses_loader.cpp 
+   gadget_loader.cpp flash_loader.cpp)
+
+IF (SDF_SUPPORT)
+ add_definitions(-DSDF_SUPPORT)
+ SET(SIMU_LOADERS_SRCS ${SIMU_LOADERS_SRCS} sdf_loader.cpp multidark_loader.cpp)
+ENDIF (SDF_SUPPORT)
+
+
+add_library(simu_loaders ${SIMU_LOADERS_SRCS})

c_src/prep/loaders/basic_loader.cpp

@@ -0,0 +1,149 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/basic_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <iostream>
+#include <boost/format.hpp>
+#include <fstream>
+#include <CosmoTool/yorick.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+using boost::format;
+using boost::str;
+
+class BasicGroupLoader: public SimulationLoader
+{
+private:
+  string storage;
+  SimuData *header;
+  int flags, numFiles;
+public:
+  BasicGroupLoader(const string& storage_path, SimuData *header, int flags, int numfiles)
+  {
+    this->header = header;
+    this->storage = storage_path;
+    this->flags = flags;
+    this->numFiles = numfiles;
+  }
+  
+  SimuData *getHeader() {
+    return header;
+  }
+  
+  int num_files() {
+    return numFiles;
+  }
+  
+  SimuData *loadFile(int id) {    
+    if (id < 0 || id >= numFiles)
+      return 0;
+
+    SimuData *simu = new SimuData;
+    uint32_t *dimlist, dimrank;
+    string fname;
+
+    simu->time = header->time;
+    simu->TotalNumPart = header->NumPart;
+    simu->Omega_M = header->Omega_M;
+    simu->Omega_Lambda = header->Omega_Lambda;
+    simu->NumPart = -1;
+
+    if (flags & NEED_POSITION)
+      {
+	loadArray(str(format("%s/x_%d.nc") % storage % id), 
+		  simu->Pos[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	simu->NumPart = dimlist[0];
+	
+	loadArray(str(format("%s/y_%d.nc") % storage % id), 
+		  simu->Pos[1], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/z_%d.nc") % storage % id), 
+		  simu->Pos[2], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+      }
+
+    if (flags & NEED_VELOCITY)
+      {
+	loadArray(str(format("%s/vx_%d.nc") % storage % id), 
+		  simu->Vel[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	if (simu->NumPart < 0)
+	  simu->NumPart = dimlist[0];
+
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/vy_%d.nc") % storage % id), 
+		  simu->Vel[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/vz_%d.nc") % storage % id), 
+		  simu->Vel[2], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+      }
+
+    if (flags & NEED_GADGET_ID)
+      {
+	loadArray(str(format("%s/id_%d.nc") % storage % id), 
+		  simu->Id, dimlist, dimrank);
+	assert(dimrank == 1);
+	if (simu->NumPart < 0)
+	  simu->NumPart = dimlist[0];
+
+	assert(simu->NumPart == dimlist[0]);
+      }
+    
+    return simu;
+  }
+
+  ~BasicGroupLoader()
+  {
+    delete header;
+  }
+};
+
+SimulationLoader *basicGroupLoader(const std::string& simupath, 
+				   int flags)
+{
+  SimuData *header;
+  ifstream f;
+  string header_path = simupath + "/header.txt";
+  int numFiles;
+
+  header = new SimuData;
+  f.open(header_path.c_str());
+
+  f >> header->time
+    >> header->Omega_M
+    >> header->Omega_Lambda 
+    >> header->NumPart
+    >> numFiles;
+
+  return new BasicGroupLoader(simupath, header, flags, numFiles);
+}

c_src/prep/loaders/flash_loader.cpp

@@ -0,0 +1,134 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/flash_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#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;
+  SimulationPreprocessor *preproc;
+public:
+  FlashLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), gadget_header(header), preproc(p)
+  {
+  }
+  
+  ~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);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *flashLoader(const std::string& snapshot, int flags, SimulationPreprocessor *p)
+{
+  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, p);
+}

c_src/prep/loaders/gadget_loader.cpp

@@ -0,0 +1,151 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/gadget_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <vector>
+#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;
+  int gadgetFormat;
+  SimuData *gadget_header;
+  string snapshot_name;
+  SimulationPreprocessor *preproc;
+public:
+  GadgetLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num, double unit, int gadgetFormat, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), unitMpc(1/unit), gadget_header(header), gadgetFormat(gadgetFormat), preproc(p)
+  {
+  }
+  
+  ~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, gadgetFormat);
+    else
+      d = loadGadgetMulti(snapshot_name.c_str(), id, load_flags, gadgetFormat);
+
+    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);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags, int gadgetFormat, SimulationPreprocessor *p)
+{
+  bool singleFile = false;
+  int num_files;
+  SimuData *d;
+
+  cout << " File to load is:" << snapshot.c_str() << endl;
+
+  try
+    {
+      d = loadGadgetMulti(snapshot.c_str(), -1, 0, gadgetFormat);
+      singleFile = true;
+      num_files = 1;
+    }
+  catch (const NoSuchFileException& e)
+    {
+      try
+        {
+          d = loadGadgetMulti(snapshot.c_str(), 0, 0, gadgetFormat);
+	  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, gadgetFormat)) != 0)
+	    {
+	      num_files++;
+	      delete d;
+	    }
+	}
+      catch(const NoSuchFileException& e)
+	{
+	}
+    }
+    
+  return new GadgetLoader(snapshot, header, flags, singleFile, num_files, Mpc_unitLength, gadgetFormat, p);
+}

c_src/prep/loaders/multidark_loader.cpp

@@ -0,0 +1,177 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/multidark_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <boost/format.hpp>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using boost::format;
+using namespace std;
+using namespace CosmoTool;
+
+static const double one_kpc = 3.08567802e16; /* km */
+static const double one_Gyr = 3.1558149984e16; /* sec */
+
+class MultiDarkLoader: public SimulationLoader
+{
+protected:
+  SimuData *header;
+  string darkname;
+  SimulationPreprocessor *preproc;
+public:
+  MultiDarkLoader(const std::string& name, SimuData *h, SimulationPreprocessor *p)
+    : preproc(p), 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 = new SimuData;
+
+    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;
+
+    long estimated = (preproc == 0) ? simu->NumPart : preproc->getEstimatedPostprocessed(simu->NumPart);
+    long allocated = estimated; 
+
+    for (int k = 0; k < 3; k++) {
+      simu->Pos[k] = new float[allocated];
+      simu->Vel[k] = new float[allocated];
+    }
+    simu->Id = new int64_t[allocated];
+    long *uniqueID = new long[allocated];
+    long *index = new long[allocated];
+
+    double tempData;
+
+    double shift = 0.5*simu->BoxSize;
+    double rescale_position = simu->Hubble*1.e-5*100./simu->time;
+    double rescale_velocity = one_kpc/one_Gyr;
+#define INFINITY std::numeric_limits<float>::max()
+    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
+
+    cout << "loading multidark particles" << endl;
+    long actualNumPart = 0;
+
+    for (long i = 0; i < simu->NumPart; i++) {
+      SingleParticle p;
+
+      fp >> p.ID >> p.Pos[0] >> p.Pos[1]
+         >> p.Pos[2] >> p.Vel[2] >> p.Vel[1] >> p.Vel[0] >> tempData;
+
+      if (p.ID == -99 && 
+          p.Pos[0] == -99 && p.Pos[1] == -99 && 
+          p.Pos[2] == -99 && p.Vel[2] == -99)
+        break;
+
+      //p.Pos[0] = p.Pos[0]*rescale_position + shift;
+      //p.Pos[1] = p.Pos[1]*rescale_position + shift;
+      //p.Pos[2] = p.Pos[2]*rescale_position + shift;
+      //p.Vel[2] = p.Vel[2]*rescale_velocity;
+
+      // enforce box size in case of roundoff error
+      for (int k = 0; k < 3; k++) {
+        if (p.Pos[k] < 0) p.Pos[k] += simu->BoxSize;
+        if (p.Pos[k] >= simu->BoxSize) p.Pos[k] -= simu->BoxSize;
+      }
+
+      if (preproc != 0 && !preproc->accept(p))
+        continue;
+
+      copyParticleToSimu(p, simu, actualNumPart);
+      uniqueID[actualNumPart]= p.ID;
+      index[actualNumPart] = i;
+      actualNumPart++;
+      if (actualNumPart == allocated)
+        {
+          allocated += (estimated+9)/10;
+          reallocSimu(simu, allocated); 
+          reallocArray(uniqueID, allocated, actualNumPart);
+          reallocArray(index, allocated, actualNumPart);
+        }
+
+       for (int k = 0; k < 3; k++) {
+         min_pos[k] = std::min(min_pos[k], p.Pos[k]);
+         max_pos[k] = std::max(max_pos[k], p.Pos[k]);
+       }
+    }
+    for (int k = 0; k < 3; k++) cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
+
+    applyTransformations(simu);
+    simu->NumPart = actualNumPart;
+    simu->TotalNumPart = actualNumPart;
+    simu->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+    simu->new_attribute("index", index, delete_adaptor<long>);
+    return simu;
+  }
+};
+
+SimulationLoader *multidarkLoader(const string& multidarkname, SimulationPreprocessor *p)
+{
+  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, p);
+}
+
+  
+

c_src/prep/loaders/ramses_loader.cpp

@@ -0,0 +1,183 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/ramses_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+// ben edit
+
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <sys/types.h>
+#include <regex.h>
+#include <cstring>
+#include <cstdlib>
+#include <cassert>
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <map>
+
+
+// end ben edit
+
+
+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;
+  SimulationPreprocessor *preproc;
+public:
+  RamsesLoader(const string& basename, int baseid, bool dp, SimuData *header, int flags, int _num, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), _num_files(_num), double_precision(dp),
+      ramses_header(header), preproc(p)
+  {
+  }
+  
+  ~RamsesLoader()
+  {
+    delete ramses_header;
+  }
+  
+  SimuData *getHeader() {
+    return ramses_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+
+// cludgy hack until I can get baseid working in this function... the problem is with   SimuData *loadFile(int id) ... just need to learn a bit more C++ ~ Ben
+
+    string baseidstr = snapshot_name.c_str();
+    unsigned found = baseidstr.find_last_of("/");
+    baseidstr = baseidstr.substr(found-5,found);
+    baseidstr = baseidstr.substr(0,5);  // remove trailing slash
+    baseid = atoi(baseidstr.c_str());
+    
+    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);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags, SimulationPreprocessor *p)
+{
+  SimuData *d, *header;
+  int num_files = 0; // how many particle files are there?
+
+  header = loadRamsesSimu(snapshot.c_str(), baseid, 0, double_precision, 0);
+
+  if (header == 0)
+    return 0;
+
+
+// count number of CPU's for this output... kinda copy/paste of loadRamses.cpp in cosmotool/src
+
+  ostringstream ss_fname;  
+  ss_fname << snapshot.c_str() << "/info_" << setfill('0') << setw(5) << baseid << ".txt";
+  cout << "Opening info file " << ss_fname.str() << " to find cpu number" << endl;
+  ifstream infile(ss_fname.str().c_str());
+  if (!infile)
+    return 0;
+
+
+  int err;
+  regex_t unit_l_rx;
+
+  //  const char *pattern = "^unit_l[ ]*=[ ]*([0-9\\.E+\\-]+)";
+  const char *pattern = "^([A-Za-z_]+)[ ]*=[ ]*([0-9\\.E+\\-]+)";
+
+  err = regcomp (&unit_l_rx, pattern, REG_EXTENDED);
+  cout << unit_l_rx.re_nsub << endl;
+  if (err)
+    {
+      char errString[255];
+      regerror(err, &unit_l_rx, errString, sizeof(errString));
+      cout << errString << endl;
+      abort();
+    }
+
+  map<string,double> infoMap;
+  string line;
+  while (getline(infile, line))
+    {
+      regmatch_t allMatch[4];
+      if (!regexec(&unit_l_rx, line.c_str(), 4, allMatch, 0))
+	{
+	  uint32_t start0 = allMatch[1].rm_so, end0 = allMatch[1].rm_eo;
+	  uint32_t start1 = allMatch[2].rm_so, end1 = allMatch[2].rm_eo;
+
+	  string keyword = line.substr(start0, end0-start0);
+	  istringstream iss(line.substr(start1, end1-start1));
+	  double unitLength;
+
+	  iss >> unitLength;
+
+	  infoMap[keyword] = unitLength;
+	}
+    }
+  
+  regfree(&unit_l_rx);
+
+
+  num_files = infoMap["ncpu"];
+
+
+  return new RamsesLoader(snapshot, baseid, double_precision, header, flags, num_files, p);
+}
+

c_src/prep/loaders/sdf_loader.cpp

@@ -0,0 +1,308 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/sdf_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <iostream>
+#include <boost/format.hpp>
+#include <vector>
+#include <cassert>
+#include <string>
+#include "sdfloader_internal.hpp"
+#include "simulation_loader.hpp"
+#include <libSDF/SDF.h>
+
+#undef SDFgetfloatOrDie
+#undef SDFgetdoubleOrDie
+#undef SDFgetintOrDie
+
+using boost::format;
+
+using namespace std;
+using namespace CosmoTool;
+
+static const double one_kpc = 3.08567802e16; /* km */
+static const double one_Gyr = 3.1558149984e16; /* sec */
+
+static void SDFgetfloatOrDie(SDF *sdfp, const char *name, float *v)
+{
+   if( SDFgetfloat(sdfp, name, v) ) {
+     cerr << format("SDFgetfloat(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+static void SDFgetdoubleOrDie(SDF *sdfp, const char *name, double *v)
+{
+   if( SDFgetdouble(sdfp, name, v) ) {
+     cerr << format("SDFgetdouble(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+static void SDFgetintOrDie(SDF *sdfp, const char *name, int *v)
+{
+   if( SDFgetint(sdfp, name, v) ) {
+     cerr << format("SDFgetint(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+class SDFLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  double unitMpc;
+  SimuData *sdf_header;
+  SDF *sdfp;
+  SimulationPreprocessor *preproc;
+  int num_splitting;
+public:
+  SDFLoader(SDF *sdf_file, SimuData *header, int flags, int num_splitting,
+            SimulationPreprocessor *p)
+    : sdfp(sdf_file), load_flags(flags),
+      sdf_header(header), preproc(p)
+  {
+    this->num_splitting = num_splitting;
+  }
+  
+  ~SDFLoader()
+  {
+    delete sdf_header;
+  }
+  
+  SimuData *getHeader() {
+    return sdf_header;
+  }
+  
+  int num_files() {
+    return num_splitting;
+  }
+
+  int64_t getStart(int id)
+  {
+    return sdf_header->TotalNumPart * int64_t(id) / num_splitting;
+  }
+
+  int64_t getNumberInSplit(int id)
+  {
+    return getStart(id+1)-getStart(id);
+  }
+
+  void rescaleParticles(SimuData *d,
+                        double rescale_position,
+                        double rescale_velocity)
+  {
+#define INFINITY std::numeric_limits<float>::max()
+    float shift = 0.5*d->BoxSize;
+    rescale_position /= d->time;
+    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
+
+    if (d->Pos[0] != 0)
+      {
+        for (int k = 0; k < 3; k++)
+          {
+            for (int64_t i = 0; i < d->NumPart; i++)
+              {
+                d->Pos[k][i] = (d->Pos[k][i]*rescale_position + shift);
+                min_pos[k] = std::min(min_pos[k], d->Pos[k][i]);
+                max_pos[k] = std::max(max_pos[k], d->Pos[k][i]);
+              }
+            cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
+          }
+      }
+    if (d->Vel[0] != 0)
+      {
+        for (int k = 0; k < 3; k++)
+          {
+            for (int64_t i = 0; i < d->NumPart; i++)
+              {
+                d->Vel[k][i] *= rescale_velocity;
+              }
+          }
+      }
+  }
+
+
+  void enforceBoxSize(SimuData *d)
+  {
+    for (int k = 0; k < 3; k++)
+      {
+        for (int64_t i = 0; i < d->NumPart; i++)
+         {
+           while (d->Pos[k][i]<0)
+             d->Pos[k][i] += d->BoxSize;
+           while (d->Pos[k][i]>=d->BoxSize)
+             d->Pos[k][i] -= d->BoxSize;
+         }
+       }
+  }
+
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    int64_t starts[7];
+    int nobjs[7];
+    int strides[7];
+    void *addrs[7];
+    const char *names[7];
+    int ret;
+
+    if (id >= num_splitting)
+      return 0;
+
+    d = new SimuData;
+    d->BoxSize = sdf_header->BoxSize;
+    d->time = sdf_header->time;
+    d->Hubble = sdf_header->Hubble;
+    d->Omega_M = sdf_header->Omega_M;
+    d->Omega_Lambda = sdf_header->Omega_Lambda;
+    
+    d->NumPart = getNumberInSplit(id);
+
+    int64_t numPartToLoad = getNumberInSplit(id);
+    int64_t start = getStart(id);
+    int k = 0;
+#define SETUP_READ(name, vec) \
+    names[k] = name, starts[k] = start, nobjs[k] = numPartToLoad, strides[k] = sizeof(vec[0]), addrs[k] = vec, k++;
+
+    if (load_flags & NEED_POSITION)
+      {
+        const char *name_template[3] = { "x", "y", "z" };
+        for (int q = 0; q < 3; q++)
+          {
+            d->Pos[q] = new float[numPartToLoad];
+            SETUP_READ(name_template[q], d->Pos[q]);
+          }
+      }
+    if (load_flags & NEED_VELOCITY)
+      {
+        const char *name_template[3] = { "vx", "vy", "vz" };
+        for (int q = 0; q < 3; q++)
+          {
+            d->Vel[q] = new float[numPartToLoad];
+            SETUP_READ(name_template[q], d->Vel[q]);
+          }
+      }
+    if (load_flags & NEED_GADGET_ID)
+      {
+        d->Id = new int64_t[numPartToLoad];
+        SETUP_READ("ident", d->Id);
+      }
+#undef SETUP_READ
+
+    ret = SDFseekrdvecsarr(sdfp, k, (char**)names, starts, nobjs, addrs, strides);
+    if (ret != 0)
+      {
+        cerr << format("Splitting %d/%d, SDF read failure: %s") % id % num_splitting % SDFerrstring << endl;
+        abort();
+      }
+
+    if (load_flags & (NEED_POSITION | NEED_VELOCITY))
+      rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr);
+
+    enforceBoxSize(d);
+    applyTransformations(d);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+SimulationLoader *sdfLoader(const std::string& snapshot, int flags, 
+                            int num_splitting,
+                            SimulationPreprocessor *p)
+{
+  SDF *sdfp;
+  int fileversion;
+  SimuData *hdr;
+  int64_t gnobj;
+
+  sdfp = SDFopen(0, snapshot.c_str());
+  if (sdfp == 0)
+    {
+      return 0;
+    }
+   cout << "Loading SDF with artificial splitting " << num_splitting << endl;
+  hdr = new SimuData;
+  
+  SDFgetintOrDefault(sdfp, "version", &fileversion, 1);
+  double h0;
+  if (fileversion == 1)
+    {
+      SDFgetfloatOrDie(sdfp, "Omega0", &hdr->Omega_M);
+      SDFgetfloatOrDie(sdfp, "Lambda_prime", &hdr->Omega_Lambda);
+      SDFgetfloatOrDie(sdfp, "H0", &hdr->Hubble);
+      hdr->Hubble *= 1000.*(one_kpc/one_Gyr);
+      h0 = hdr->Hubble / 100.;
+    }
+  else
+    {
+      float Or, Of;
+      SDFgetfloatOrDie(sdfp, "Omega0_m", &hdr->Omega_M);
+      SDFgetfloatOrDie(sdfp, "Omega0_lambda", &hdr->Omega_Lambda);
+      SDFgetfloatOrDie(sdfp, "Omega0_r", &Or);
+      SDFgetfloatOrDie(sdfp, "Omega0_fld", &Of);
+
+      hdr->Omega_M += Or;
+      hdr->Omega_Lambda += Of;
+      SDFgetfloatOrDie(sdfp, "h_100", &hdr->Hubble);
+      hdr->Hubble *= 100.;
+      h0 = hdr->Hubble / 100.; 
+    }
+
+  if (SDFhasname("R0", sdfp))
+    {
+      double R0;
+
+      SDFgetdoubleOrDie(sdfp, "R0", &R0);
+      hdr->BoxSize = 2.0*0.001*R0*h0;
+    }
+
+  if (SDFhasname("Rx", sdfp))
+    {
+      double R0;
+
+      SDFgetdoubleOrDie(sdfp, "Rx", &R0);
+      hdr->BoxSize = 2.0*0.001*R0*h0;
+    }
+
+  if (SDFhasname("redshift", sdfp))
+    {
+      double redshift;
+
+      SDFgetdoubleOrDie(sdfp, "redshift", &redshift);
+      hdr->time = 1/(1+redshift);
+    }
+   
+  if (SDFgetint64(sdfp, "npart", &gnobj))
+    {
+      gnobj = SDFnrecs("x", sdfp);
+      cerr << format("[Warning] No 'npart' found in SDF file '%s', guessing npart=%ld from SDFnrecs") % snapshot % gnobj << endl;
+    }
+  hdr->NumPart = hdr->TotalNumPart = gnobj;
+
+  return new SDFLoader(sdfp, hdr, flags, num_splitting, p);
+}
+
+
+void sdfLoaderInit(int& argc, char **& argv)
+{
+}

c_src/prep/loaders/sdfloader_internal.hpp

@@ -0,0 +1,27 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/sdfloader_internal.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef __VOID_SDF_LOADER_INTERNAL_HPP
+#define __VOID_SDF_LOADER_INTERNAL_HPP
+
+void sdfLoaderInit(int& argc, char **&argv);
+
+#endif

c_src/prep/loaders/simulation_loader.cpp

@@ -0,0 +1,110 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/simulation_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cmath>
+#include <CosmoTool/loadSimu.hpp>
+#include "simulation_loader.hpp"
+#ifdef SDF_SUPPORT
+#include "sdfloader_internal.hpp"
+#endif
+
+using std::min;
+using namespace CosmoTool;
+
+template<typename T> void reallocArray(T *& a, long newSize, long toCopy)
+{
+  T *b = new T[newSize];
+  if (a != 0)
+    {
+      memcpy(b, a, sizeof(T)*toCopy);
+      delete[] a;
+    }
+  a = b;
+}
+
+void SimulationLoader::reallocSimu(SimuData *s, long newNumPart)
+{
+  long to_copy = min(newNumPart, s->NumPart);
+  
+  for (int j = 0; j < 3; j++)
+    {
+      reallocArray(s->Pos[j], newNumPart, to_copy);
+      reallocArray(s->Vel[j], newNumPart, to_copy);
+    }
+  reallocArray(s->Id, newNumPart, to_copy);
+  reallocArray(s->type, newNumPart, to_copy);
+}
+
+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.;
+    }
+}
+
+
+void SimulationLoader::basicPreprocessing(SimuData *d, 
+                                          SimulationPreprocessor *preproc)
+{
+  if (preproc == 0)
+    return;
+
+  long numAccepted = 0;
+  bool *accepted = new bool[d->NumPart];
+  for (long i = 0; i < d->NumPart; i++)
+    {
+      SingleParticle p;
+      
+      for (int k = 0; k < 3; k++)
+        {
+          p.Pos[k] = (d->Pos[k]) ? 0 : d->Pos[k][i];
+          p.Vel[k] = (d->Vel[k]) ? 0 : d->Vel[k][i];
+        }
+      p.ID = (d->Id) ? 0 : d->Id[i];
+      
+      accepted[i] = preproc->accept(p);
+      numAccepted += accepted[i];
+    }
+  
+  for (int k = 0; k < 3; k++)
+    {
+      filteredCopy(d->Pos[k], accepted, d->NumPart);
+      filteredCopy(d->Vel[k], accepted, d->NumPart);
+    }
+  filteredCopy(d->Id, accepted, d->NumPart);
+  filteredCopy(d->type, accepted, d->NumPart);
+  
+  filterAttribute<long>(d, "uniqueID", accepted, d->NumPart);
+  
+  d->NumPart = numAccepted;
+  delete[] accepted;
+}
+
+void simulationLoadersInit(int& argc, char **& argv)
+{
+#ifdef SDF_SUPPORT
+  sdfLoaderInit(argc, argv);
+#endif
+}

c_src/prep/loaders/simulation_loader.hpp

@@ -0,0 +1,156 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/simulation_loader.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef _MOCK_SIMULATION_LOADER_HPP
+#define _MOCK_SIMULATION_LOADER_HPP
+
+#include <string>
+#include <algorithm>
+#include <CosmoTool/loadSimu.hpp>
+
+struct SingleParticle
+{
+  float Pos[3];
+  float Vel[3];
+  long Index;
+  long ID;
+};
+
+class SimulationPreprocessor
+{
+public:
+  SimulationPreprocessor() {}
+  virtual ~SimulationPreprocessor() {}
+
+  virtual long getEstimatedPostprocessed(long numParticles) { return numParticles; };
+  virtual bool accept(const SingleParticle& p) { return true; }
+  virtual void reset() {}
+};
+
+class SimulationLoader
+{
+protected:
+  bool do_redshift;
+  int redshift_axis;
+
+  SimulationLoader() 
+  {
+    do_redshift = false;
+    redshift_axis = 2;
+  }
+
+  template<typename T> void reallocArray(T *& a, long newSize, long toCopy)
+  {
+    T *b = new T[newSize];
+    if (a != 0)
+     {
+       std::copy(a, a+toCopy, b);
+       delete[] a;
+     }
+   a = b;
+  } 
+
+
+
+  void reallocSimu(CosmoTool::SimuData *s, long newNumPart);
+ 
+
+  void basicPreprocessing(CosmoTool::SimuData *d, SimulationPreprocessor *preproc);
+  void applyTransformations(CosmoTool::SimuData *s);
+
+  void copyParticleToSimu(const SingleParticle& p, CosmoTool::SimuData *s, long index)
+  {
+    s->Pos[0][index] = p.Pos[0];
+    s->Pos[1][index] = p.Pos[1];
+    s->Pos[2][index] = p.Pos[2];
+    if (s->Vel[0])
+      s->Vel[0][index] = p.Vel[0];
+    if (s->Vel[1])
+      s->Vel[1][index] = p.Vel[1];
+    if (s->Vel[2])
+      s->Vel[2][index] = p.Vel[2];
+    s->Id[index] = p.ID;
+  }
+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 filteredCopy(T *a, bool *accepted, long N)
+  {
+    long i = 0, j = 0;
+
+    if (a == 0)
+      return;
+
+    while (i < N)
+      {
+        if (!accepted[i])
+          {
+            i++;
+            continue;
+          }
+
+        a[j] = a[i];
+        j++;
+        i++;
+      }
+  }
+  
+  template<typename T>
+  void filterAttribute(CosmoTool::SimuData *d, const std::string& attrname, bool *accepted, long NumPart)
+  {
+    if (d->attributes.find(attrname) == d->attributes.end())
+      return;
+ 
+    long *l = d->as<T>(attrname);
+    filteredCopy<T>(l, accepted, NumPart);
+  }
+};
+
+
+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, int gadgetFormat, SimulationPreprocessor *p);
+SimulationLoader *flashLoader(const std::string& snapshot, int flags, SimulationPreprocessor *p);
+SimulationLoader *multidarkLoader(const std::string& snapshot, SimulationPreprocessor *p);
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags, SimulationPreprocessor *p);
+SimulationLoader *sdfLoader(const std::string& snapshot, int flags, int num_splitting, SimulationPreprocessor *p);
+
+/* This is required for some parallel I/O handler (thus MPI beneath it) */
+void simulationLoadersInit(int& argc, char **& argv);
+
+#endif