Fixed merge

lib/Makefile

@@ -1,5 +1,5 @@
 SHLIBS= libCosmoTool.so
-SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp  fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp octTree.cpp
+SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp  fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp octTree.cpp loadGadget.cpp
 LIBS= -lnetcdf_c++ -lnetcdf -lgsl -lgslcblas -lm
 
 include config.mk
@@ -8,7 +8,7 @@ VPATH=../src
 
 all: $(SHLIBS)
 
-libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o octTree.o
+libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o octTree.o loadGadget.o
 
 depend: $(SOURCES)
 	@echo "[DEPENDS] $^"

sample/testkd2.cpp

@@ -2,10 +2,12 @@
 #include <cassert>
 #include <cstdlib>
 #include <iostream>
+#include <fstream>
+#define __KD_TREE_NUMNODES
 #include "mykdtree.hpp"
 
 #define NTRY 3
-#define ND 2
+#define ND 3
 
 using namespace std;
 using namespace CosmoTool;
@@ -61,6 +63,7 @@ int main()
   cout << "Check consistency..." << endl;
   MyCell **ngb = new MyCell *[12];
   
+  ofstream fngb("nearest.txt");
   for (int k = 0; k < NTRY; k++) {
     cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl;
     tree.getNearestNeighbours(xc[k], 12, ngb);
@@ -70,7 +73,7 @@ int main()
 	double d2 = 0;
 	for (int l = 0; l < 3; l++)
 	  d2 += ({double delta = xc[k][l] - ngb[i]->coord[l]; delta*delta;});
-	cout << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl;
+	fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl;
       }
   }
 

src/fortran.cpp

@@ -49,6 +49,25 @@ void UnformattedRead::setCheckpointSize(CheckpointSize cs)
   cSize = cs;
 }
     
+void UnformattedRead::skip(int64_t off)
+  throw (InvalidUnformattedAccess)
+{
+  if (checkPointAccum == 0 && checkPointRef == 0)
+    {
+      // We are not in a checked block
+      f->seekg(off, ios::cur);
+      return;
+    }
+  if (off < 0)
+    throw InvalidUnformattedAccess();
+  
+  if ((checkPointAccum+off) > checkPointRef)
+       throw InvalidUnformattedAccess();
+
+  f->seekg(off, ios::cur);
+  checkPointAccum += off;
+}
+
 void UnformattedRead::beginCheckpoint()
   throw (InvalidUnformattedAccess,EndOfFileException)
 {
@@ -65,11 +84,15 @@ void UnformattedRead::beginCheckpoint()
     throw EndOfFileException();
 }
 
-void UnformattedRead::endCheckpoint()
+void UnformattedRead::endCheckpoint(bool autodrop)
   throw (InvalidUnformattedAccess)
 {
   if (checkPointRef != checkPointAccum)
-    throw InvalidUnformattedAccess();
+    {
+      if (!autodrop || checkPointAccum > checkPointRef)
+	throw InvalidUnformattedAccess();
+      f->seekg(checkPointRef-checkPointAccum, ios::cur);
+    }
 
   int64_t oldCheckPoint = checkPointRef;
 

src/fortran.hpp

@@ -41,7 +41,7 @@ namespace CosmoTool
     
     void beginCheckpoint()
       throw (InvalidUnformattedAccess,EndOfFileException);
-    void endCheckpoint()
+    void endCheckpoint(bool autodrop = false)
       throw (InvalidUnformattedAccess);
     
     double readReal64()
@@ -53,6 +53,9 @@ namespace CosmoTool
     int64_t readInt64()
       throw (InvalidUnformattedAccess);
 
+    void skip(int64_t off)
+      throw (InvalidUnformattedAccess);
+
   protected:
     bool swapOrdering;
     CheckpointSize cSize;

src/interpolate.hpp

@@ -1,5 +1,5 @@
-#ifndef __MAK_INTERPOLATE_HPP
-#define __MAK_INTERPOLATE_HPP
+#ifndef __CTOOL_INTERPOLATE_HPP
+#define __CTOOL_INTERPOLATE_HPP
 
 #include "config.hpp"
 #include <inttypes.h>
@@ -26,6 +26,7 @@ namespace CosmoTool
     uint32_t getNumPoints() const;
     void fillWithXY(double *x, double *y) const;
     double getMaxX() const;
+    double getXi(int i) const { return spline->x[i]; }
   protected:
     gsl_interp_accel *accel_interp;
     gsl_spline *spline;

src/kdtree_leaf.hpp

@@ -0,0 +1,105 @@
+#ifndef __LEAF_KDTREE_HPP
+#define __LEAF_KDTREE_HPP
+
+#include <cmath>
+#include "config.hpp"
+#include "bqueue.hpp"
+
+namespace CosmoTool {
+
+  template<int N, typename CType = ComputePrecision> 
+  struct KDLeafDef
+  {
+    typedef CType CoordType;
+    typedef float KDLeafCoordinates[N];
+  };
+  
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDLeafCell
+  {
+    bool active;
+    ValType val;
+    typename KDLeafDef<N,CType>::KDLeafCoordinates coord;
+  };
+
+  class NotEnoughCells: public Exception
+  {
+  public:
+    NotEnoughCells() : Exception() {}
+    ~NotEnoughCells() throw () {}
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDLeafTreeNode
+  {
+    bool leaf;
+    union {
+      KDLeafCell<N,ValType,CType> *value;
+      KDLeafTreeNode<N,ValType,CType> *children[2];
+    };
+    typename KDLeafDef<N,CType>::KDLeafCoordinates minBound, maxBound;
+#ifdef __KDLEAF_TREE_NUMNODES
+    uint32_t numNodes;
+#endif
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  class KDLeafTree
+  {
+  public:
+    typedef typename KDLeafDef<N,CType>::CoordType CoordType;
+    typedef typename KDLeafDef<N>::KDLeafCoordinates coords;
+    typedef KDLeafCell<N,ValType,CType> Cell;
+    typedef KDLeafTreeNode<N,ValType,CType> Node;
+
+    KDLeafTree(Cell *cells, uint32_t Ncells);
+    ~KDLeafTree();
+
+    Node *getRoot() { return root; }
+
+    void optimize();
+
+    Node *getAllNodes() { return nodes; }
+    uint32_t getNumNodes() const { return lastNode; }
+
+    uint32_t countActives() const;
+
+    CoordType computeDistance(const Cell *cell, const coords& x) const;
+
+#ifdef __KDLEAF_TREE_NUMNODES
+    uint32_t getNumberInNode(const Node *n) const { return n->numNodes; }
+#else
+    uint32_t getNumberInNode(const Node *n) const {
+      if (n == 0) 
+        return 0;
+      return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
+    }
+#endif
+
+    double countInRange(CType sLo, CType sHi, Node *root1 = 0, Node *root2 = 0) const;
+
+  protected:
+    Node *nodes;
+    uint32_t numNodes, numCells;
+    uint32_t lastNode;
+
+    Node *root;
+    Cell **sortingHelper;
+
+    Node *buildTree(Cell **cell0,
+		    uint32_t NumCells,
+		    uint32_t depth,
+		    coords minBound,
+		    coords maxBound);
+    
+    double recursiveCountInRange(Node *na, Node *nb, CType sLo, CType sHi) const;
+  };
+
+  template<int N, typename T, typename CType>
+  uint32_t gatherActiveCells(KDLeafCell<N,T,CType> **cells, uint32_t numCells);
+
+};
+
+#include "kdtree_leaf.tcc"
+
+#endif

src/kdtree_leaf.tcc

@@ -0,0 +1,308 @@
+#include <cstring>
+#include <algorithm>
+#include <limits>
+#include <iostream>
+#include <cassert>
+
+namespace CosmoTool {
+
+  template<int N, typename ValType, typename CType>
+  class CellCompare
+  {
+  public:
+    CellCompare(int k)
+    {
+      rank = k;
+    }
+
+    bool operator()(const KDLeafCell<N,ValType,CType> *a, const KDLeafCell<N,ValType,CType> *b) const
+    {
+      return (a->coord[rank] < b->coord[rank]);
+    }
+  protected:
+    int rank;
+  };
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTree<N,ValType,CType>::~KDLeafTree()
+  {
+  }
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTree<N,ValType,CType>::KDLeafTree(Cell *cells, uint32_t Ncells)
+  {
+    numNodes = Ncells*3;
+    numCells = Ncells;
+    nodes = new Node[numNodes];
+
+    sortingHelper = new Cell *[Ncells];
+    for (uint32_t i = 0; i < Ncells; i++)
+	sortingHelper[i] = &cells[i];    
+
+    optimize();
+  }
+
+  template<int N, typename ValType, typename CType>
+  void KDLeafTree<N,ValType,CType>::optimize()
+  {
+    coords absoluteMin, absoluteMax;
+
+    std::cout << "Optimizing the tree..." << std::endl;
+    uint32_t activeCells = gatherActiveCells(sortingHelper, numCells);
+    std::cout << "  number of active cells = " << activeCells << std::endl;
+
+    lastNode = 0;
+    for (int i = 0; i < N; i++)
+      {
+	absoluteMin[i] = std::numeric_limits<typeof (absoluteMin[0])>::max();
+	absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max();
+      }
+    // Find min and max corner
+    for (uint32_t i = 0; i < activeCells; i++)
+      {
+        KDLeafCell<N,ValType,CType> *cell = sortingHelper[i];
+
+        for (int k = 0; k < N; k++) {
+          if (cell->coord[k] < absoluteMin[k])
+            absoluteMin[k] = cell->coord[k];
+          if (cell->coord[k] > absoluteMax[k])
+            absoluteMax[k] = cell->coord[k];
+        }
+      }
+    
+    std::cout << "   rebuilding the tree..." << std::endl;
+    root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);    
+    std::cout << "   done." << std::endl;
+  }
+
+  template<int N, typename ValType, typename CType>
+  uint32_t gatherActiveCells(KDLeafCell<N,ValType,CType> **cells,
+			     uint32_t Ncells)
+  {
+    uint32_t swapId = Ncells-1;
+    uint32_t i = 0;
+
+    while (!cells[swapId]->active && swapId > 0)
+      swapId--;
+
+    while (i < swapId)
+      {
+	if (!cells[i]->active)
+	  {
+	    std::swap(cells[i], cells[swapId]);
+	    while (!cells[swapId]->active && swapId > i)
+	      {
+		swapId--;
+	      }
+	  }
+	i++;
+      }
+    return swapId+1;
+  }
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTreeNode<N,ValType,CType> *
+        KDLeafTree<N,ValType,CType>::buildTree(Cell **cell0,
+					       uint32_t Ncells,
+					       uint32_t depth,
+					       coords minBound,
+					       coords maxBound)
+  {
+    if (Ncells == 0)
+      return 0;
+
+    int axis = depth % N;
+    assert(lastNode != numNodes);
+    Node *node = &nodes[lastNode++];
+    uint32_t mid = Ncells/2;
+    coords tmpBound;
+    
+    // Isolate the environment
+    {
+      CellCompare<N,ValType,CType> compare(axis);
+      std::sort(cell0, cell0+Ncells, compare);
+    }
+
+    node->leaf = false;
+    memcpy(&node->minBound[0], &minBound[0], sizeof(coords));
+    memcpy(&node->maxBound[0], &maxBound[0], sizeof(coords));
+
+    if (Ncells == 1)
+      {
+	node->leaf = true;
+	node->value = *cell0;
+
+#ifdef __KDLEAF_TREE_NUMNODES
+	node->numNodes = 1;
+#endif
+	return node;
+      }
+
+    memcpy(tmpBound, maxBound, sizeof(coords));
+    tmpBound[axis] = (*(cell0+mid))->coord[axis];
+    depth++;
+    node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
+    
+    memcpy(tmpBound, minBound, sizeof(coords));
+    tmpBound[axis] = (*(cell0+mid))->coord[axis];
+    node->children[1] = buildTree(cell0+mid, Ncells-mid, depth, 
+				  tmpBound, maxBound);
+
+#ifdef __KDLEAF_TREE_NUMNODES
+    node->numNodes = (node->children[0] != 0) ? node->children[0]->numNodes : 0;
+    node->numNodes += (node->children[1] != 0) ? node->children[1]->numNodes : 0;
+#endif
+
+    return node;
+  }
+
+  template<int N, typename ValType, typename CType>
+  uint32_t KDLeafTree<N,ValType,CType>::countActives() const
+  {
+    uint32_t numActive = 0;
+    for (uint32_t i = 0; i < lastNode; i++)
+      {
+	if (nodes[i].value->active)
+	  numActive++;
+      }
+    return numActive;
+  }
+
+  template<int N, typename ValType, typename CType>
+  typename KDLeafDef<N,CType>::CoordType
+  KDLeafTree<N,ValType,CType>::computeDistance(const Cell *cell, const coords& x) const
+  {
+    CoordType d2 = 0;
+
+    for (int i = 0; i < N; i++)
+      {
+	CoordType delta = cell->coord[i] - x[i];
+	d2 += delta*delta;
+      }
+    return d2;
+  }
+
+ template<int N, typename ValType, typename CType>
+ double KDLeafTree<N,ValType,CType>::countInRange(CType sLo, CType sHigh, Node *root1, Node *root2) const
+ {
+   double result = recursiveCountInRange((root1 == 0) ? root : root1,
+				(root2 == 0) ? root : root2, 
+				sLo*sLo, sHigh*sHigh); 
+   return result;
+ }
+
+ 
+ template<int N, typename ValType, typename CType>
+ double KDLeafTree<N,ValType,CType>::recursiveCountInRange(Node *na, Node *nb,
+                                                         CType sLo, CType sHi) const
+ {
+   assert(nb != 0);
+   if (na == 0)
+     {
+       return 0;
+     }
+   
+   uint32_t numNa = getNumberInNode(na);
+   uint32_t numNb = getNumberInNode(nb);
+   double Cleft, Cright;
+   CType minDist, maxDist;
+
+   if (numNa == 1 && numNb == 1)
+     {
+       assert(na->leaf && nb->leaf);
+       CType ab_dist = computeDistance(na->value, nb->value->coord);
+       if (ab_dist >= sLo && ab_dist < sHi)
+	 return 1;
+       else
+	 return 0;
+     }
+   assert(numNa > 1 || numNb > 1);
+   
+   bool overlapping_a = true, overlapping_b = true;
+   for (int k = 0; k < N; k++) 
+     {
+       bool min_a_in_B = 
+	 ((na->minBound[k] >= nb->minBound[k] && 
+	   na->minBound[k] <= nb->maxBound[k]));
+       bool max_a_in_B =
+	 ((na->maxBound[k] >= nb->minBound[k] && 
+	   na->maxBound[k] <= nb->maxBound[k]));
+       bool min_b_in_A = 
+	 ((nb->minBound[k] >= na->minBound[k] && 
+	   nb->minBound[k] <= na->maxBound[k]));
+       bool max_b_in_A = 
+	 ((nb->maxBound[k] >= na->minBound[k] && 
+	   nb->maxBound[k] <= na->maxBound[k]));
+
+      if (!min_a_in_B && !max_a_in_B)
+         overlapping_a = false;
+      if (!min_b_in_A && !max_b_in_A)
+         overlapping_b = false;
+     }
+
+   if (overlapping_a || overlapping_b)
+     {
+        minDist = 0;
+        maxDist = 0;
+        for (int k = 0; k < N; k++)
+        {
+          CType delta = max(nb->maxBound[k]-na->minBound[k],na->maxBound[k]-nb->minBound[k]);
+          maxDist += delta*delta;
+        }
+     }
+   else
+     {
+       minDist = maxDist = 0;
+       for (int k = 0; k < N; k++)
+         {
+           CType delta2;
+           delta2 = max(nb->maxBound[k]-na->minBound[k],
+			na->maxBound[k]-nb->minBound[k]);
+           maxDist += delta2*delta2;
+         }
+       // mins and maxs
+       CType minmax[N][2];
+       for (int k = 0; k < N; k++)
+	 {
+	   if (na->minBound[k] < nb->minBound[k])
+	     {
+	       minmax[k][1] = na->maxBound[k];
+	       minmax[k][0] = nb->minBound[k];
+	     }
+	   else
+	     {
+	       minmax[k][1] = nb->maxBound[k];
+	       minmax[k][0] = na->minBound[k];
+	     }
+	 }
+       for (int k = 0; k < N; k++)
+	 {
+	   CType delta = max(minmax[k][0]-minmax[k][1], 0.);
+	   minDist += delta*delta;
+	 }
+     }
+   
+   if (minDist >= sHi)
+     return 0;
+   if (maxDist < sLo)
+     return 0;
+   
+   if (sLo <= minDist && maxDist < sHi)
+     return ((double)numNa)*numNb;
+   
+   if (numNa < numNb) 
+     {
+       assert(!nb->leaf);
+       Cleft = recursiveCountInRange(nb->children[0], na, sLo, sHi);
+       Cright = recursiveCountInRange(nb->children[1], na, sLo, sHi);
+     }
+   else
+     {
+       assert(!na->leaf);
+       Cleft = recursiveCountInRange(na->children[0], nb, sLo, sHi);
+       Cright = recursiveCountInRange(na->children[1], nb, sLo, sHi);
+     }
+   return Cleft+Cright;
+ }
+
+};

src/loadGadget.cpp

@@ -0,0 +1,217 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "load_data.hpp"
+#include "loadGadget.hpp"
+#include "fortran.hpp"
+
+using namespace CosmoTool;
+
+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();
+  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;
+}
+
+
+
+
+SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
+{
+  SimuData *data;
+  int p, n;
+  UnformattedRead *f;
+  GadgetHeader h;
+
+  if (id >= 0) {
+    int numDigits = 1;
+    int mul = 1;
+    while (mul < id)
+      {
+	mul *= 10;
+	numDigits++;
+      }
+    size_t len = numDigits+2+strlen(fname);
+    char *out_fname = new char[numDigits+2+strlen(fname)];
+
+    if (snprintf(out_fname, len, "%s.%d", fname, id) != len)
+      abort();
+     
+    f = new UnformattedRead(out_fname);
+    if (f == 0)
+      return 0;
+
+    delete out_fname;
+
+  } else {
+
+    f = new UnformattedRead(fname);
+    if (f == 0)
+      return 0;
+
+  }
+
+  data = new SimuData;
+  if (data == 0) {
+    delete f;
+    return 0;
+  }
+
+  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);
+
+  long NumPart = 0, NumPartTotal = 0;
+  for(int k=0; k<5; k++)
+    {
+      NumPart += h.npart[k];
+      NumPartTotal += h.npartTotal[k];
+    }
+  data->NumPart = NumPart;
+  data->TotalNumPart = NumPartTotal;
+
+  if (loadflags & NEED_POSITION) {
+      
+    for (int i = 0; i < 3; i++) {
+      data->Pos[i] = new float[data->NumPart];
+      if (data->Pos[i] == 0) {
+	delete data;
+	return 0;
+      }
+    }
+
+    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();
+    
+  } else {
+    // Skip positions
+    f->skip(NumPart * 3 * sizeof(float) + 2*4);
+  }
+  
+  if (loadflags & NEED_VELOCITY) {
+    for (int i = 0; i < 3; i++)
+      {
+	data->Vel[i] = new float[data->NumPart];
+	if (data->Vel[i] == 0) 
+	  {
+	    delete data;
+	    return 0;
+	  }
+      }
+
+    f->beginCheckpoint();
+    for(int k = 0, p = 0; k < 5; k++) {
+      for(int n = 0; n < h.npart[k]; n++) {
+	data->Vel[p][0] = f->readReal32();
+	data->Vel[p][1] = f->readReal32();
+	data->Vel[p][2] = f->readReal32();
+	p++;
+      }
+    }
+    f->endCheckpoint();
+    
+    // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!!
+  } else {
+    // Skip velocities
+    f->skip(NumPart*3*sizeof(float)+2*4);
+  }
+
+  // Skip ids
+  if (loadflags & NEED_GADGET_ID) {
+    f->beginCheckpoint();
+    data->Id = new int[data->NumPart];
+    if (data->Id == 0)
+      {
+	delete data;
+	return 0;
+      }
+
+    for(int k = 0, p = 0; k < 6; k++)
+      {
+	for(int n = 0; n < h.npart[k]; n++)
+	  {
+	    data->Id[p] = f->readInt32();
+	    p++;
+	  }
+      }
+    f->endCheckpoint();
+  } else {
+    f->skip(2*4);
+    for (int k = 0; k < 6; k++)
+      f->skip(h.npart[k]*4);
+  }
+
+  delete f;
+
+  return data;
+}
+
+

src/loadGadget.hpp

@@ -0,0 +1,15 @@
+#ifndef __COSMO_LOAD_GADGET_HPP
+#define __COSMO_LOAD_GADGET_HPP
+
+#include "load_data.hpp"
+#include "loadSimu.hpp"
+
+namespace CosmoTool {
+  
+  PurePositionData *loadGadgetPosition(const char *fname);
+
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags);
+  
+};
+
+#endif

src/loadSimu.hpp

@@ -0,0 +1,41 @@
+#ifndef __COSMOTOOLBOX_HPP
+#define __COSMOTOOLBOX_HPP
+
+
+namespace CosmoTool
+{
+  static const int NEED_GADGET_ID = 1;
+  static const int NEED_POSITION = 2;
+  static const int NEED_VELOCITY = 4;
+  
+
+  class SimuData 
+  {
+  public:
+    float BoxSize;
+    float time;
+
+    long NumPart;
+    long TotalNumPart;
+    int *Id;
+    float *Pos[3];
+    float *Vel[3];
+  public:
+    SimuData() : Id(0),NumPart(0)  { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
+    ~SimuData() 
+    {
+      for (int j = 0; j < 3; j++)
+	{
+	  if (Pos[j])
+	    delete[] Pos[j];
+	  if (Vel[j])
+	    delete[] Vel[j];
+	}
+      if (Id)
+	delete[] Id;
+    }    
+  };
+
+};
+
+#endif

src/mykdtree.hpp

@@ -35,6 +35,9 @@ namespace CosmoTool {
     KDCell<N,ValType,CType> *value;
     KDTreeNode<N,ValType,CType> *children[2];
     typename KDDef<N,CType>::KDCoordinates minBound, maxBound;
+#ifdef __KD_TREE_NUMNODES
+    uint32_t numNodes;
+#endif
   };
 
   template<int N, typename ValType, typename CType = ComputePrecision>
@@ -105,6 +108,17 @@ namespace CosmoTool {
 
     uint32_t countActives() const;
 
+#ifdef __KD_TREE_NUMNODES
+    uint32_t getNumberInNode(const Node *n) const { return n->numNodes; }
+#else
+    uint32_t getNumberInNode(const Node *n) const {
+      if (n == 0) 
+        return 0;
+      return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
+    }
+#endif
+
+
   protected:
     Node *nodes;
     uint32_t numNodes;
@@ -124,7 +138,7 @@ namespace CosmoTool {
 				    int level)
       throw (NotEnoughCells);
 
-    CoordType computeDistance(Cell *cell, const coords& x);
+    CoordType computeDistance(const Cell *cell, const coords& x) const;
     void recursiveNearest(Node *node,
 			  int level,
 			  const coords& x,
@@ -132,6 +146,7 @@ namespace CosmoTool {
 			  Cell*& cell);
     void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
 				  int level);    
+
   };
 
   template<int N, typename T, typename CType>

src/mykdtree.tcc

@@ -2,6 +2,7 @@
 #include <algorithm>
 #include <limits>
 #include <iostream>
+#include <cassert>
 
 namespace CosmoTool {
 
@@ -53,8 +54,20 @@ namespace CosmoTool {
     lastNode = 0;
     for (int i = 0; i < N; i++)
       {
-	absoluteMin[i] = -std::numeric_limits<typeof (absoluteMin[0])>::max();
-	absoluteMax[i] = std::numeric_limits<typeof (absoluteMax[0])>::max();
+	absoluteMin[i] = std::numeric_limits<typeof (absoluteMin[0])>::max();
+	absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max();
+      }
+    // Find min and max corner
+    for (uint32_t i = 0; i < activeCells; i++)
+      {
+        KDCell<N,ValType,CType> *cell = sortingHelper[i];
+
+        for (int k = 0; k < N; k++) {
+          if (cell->coord[k] < absoluteMin[k])
+            absoluteMin[k] = cell->coord[k];
+          if (cell->coord[k] > absoluteMax[k])
+            absoluteMax[k] = cell->coord[k];
+        }
       }
     
     std::cout << "   rebuilding the tree..." << std::endl;
@@ -208,6 +221,12 @@ namespace CosmoTool {
     node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth, 
 				  tmpBound, maxBound);
 
+#ifdef __KD_TREE_NUMNODES
+    node->numNodes = (node->children[0] != 0) ? node->children[0]->numNodes : 0;
+    node->numNodes += (node->children[1] != 0) ? node->children[1]->numNodes : 0;
+    node->numNodes++;
+#endif
+
     return node;
   }
 
@@ -225,7 +244,7 @@ namespace CosmoTool {
 
   template<int N, typename ValType, typename CType>
   typename KDDef<N,CType>::CoordType
-  KDTree<N,ValType,CType>::computeDistance(Cell *cell, const coords& x)
+  KDTree<N,ValType,CType>::computeDistance(const Cell *cell, const coords& x) const
   {
     CoordType d2 = 0;
 
@@ -410,4 +429,5 @@ namespace CosmoTool {
    //   std::cout << "Traversed = " << info.traversed << std::endl;
  }
 
+
 };

src/octTree.cpp

@@ -157,3 +157,4 @@ void OctTree::insertParticle(octPtr node,
 		 particleId, maxAbsoluteDepth-1);  
   cells[node].children[octPos] = newNode;
 }
+

src/octTree.hpp

@@ -29,7 +29,6 @@ namespace CosmoTool
   class OctTree
   {
   public:
-    //Coordinates of particles must be in the [0:1] range
     OctTree(const FCoordinates *particles, octPtr numParticles,
 	    uint32_t maxTreeDepth,  uint32_t maxAbsoluteDepth,
 	    uint32_t threshold = 1);
@@ -47,6 +46,11 @@ namespace CosmoTool
       return cells[0].numberLeaves;
     }
 
+    static bool unconditioned(const FCoordinates&, octPtr, float, bool)
+    {
+      return true;
+    }
+
     template<typename FunT>
     void walkTree(FunT f)
     {
@@ -61,6 +65,8 @@ namespace CosmoTool
       walkTreeElements(f, condition, 0, rootCenter, octCoordCenter);
     }
 
+    
+
   protected:
     const FCoordinates *particles;
     octPtr numParticles;
@@ -72,6 +78,7 @@ namespace CosmoTool
     float xMin[3];
 
 
+<<<<<<< HEAD
     static bool unconditioned()
     {
       return true;
@@ -80,6 +87,10 @@ namespace CosmoTool
     template<typename FunT, typename CondT>
     void walkTreeElements(FunT f, CondT condition,
 			  octPtr node, 
+=======
+    template<typename FunT,typename CondT>
+    void walkTreeElements(FunT f, CondT condition, octPtr node, 
+>>>>>>> 37b41b5ac9b32213b865cbeddd63102f3fa0935a
 			  const OctCoords& icoord,
 			  octCoordType halfNodeLength)
     {
@@ -95,12 +106,21 @@ namespace CosmoTool
       f(realCenter, cells[node].numberLeaves, lenNorm*halfNodeLength/(float)octCoordCenter,
 	cells[node].children[0] == invalidOctCell, // True if this is a meta-node
 	false);
+<<<<<<< HEAD
 
       if (!condition(realCenter, cells[node].numberLeaves, 
 		     lenNorm*halfNodeLength/(float)octCoordCenter,
 		     cells[node].children[0] == invalidOctCell))
 	return;
       
+=======
+     
+      if (!condition(realCenter, cells[node].numberLeaves,
+		     lenNorm*halfNodeLength/(float)octCoordCenter,
+		     cells[node].children[0] == invalidOctCell))
+	return;
+       
+>>>>>>> 37b41b5ac9b32213b865cbeddd63102f3fa0935a
       for (int i = 0; i < 8; i++)
 	{
 	  octPtr newNode = cells[node].children[i];
@@ -125,8 +145,12 @@ namespace CosmoTool
 		false, true);
 	      continue;
 	    }
+<<<<<<< HEAD
 	  walkTreeElements(f, condition,
 			   cells[node].children[i], newCoord, halfNodeLength/2);
+=======
+	  walkTreeElements(f, condition, cells[node].children[i], newCoord, halfNodeLength/2);
+>>>>>>> 37b41b5ac9b32213b865cbeddd63102f3fa0935a
 	}
       
     }