Sparse grid support

src/sparseGrid.hpp

@@ -0,0 +1,87 @@
+#ifndef __VOID_SPARSE_GRID_HPP
+#define __VOID_SPARSE_GRID_HPP
+
+#include <cmath>
+#include "config.hpp"
+#include <list>
+
+namespace CosmoTool 
+{
+  
+  template <typename T>
+  struct bracketAccessor
+  {
+    typedef typename T::value_type result_type;
+
+    result_type
+    operator()(T const& V, size_t const N) const throw ()
+    {
+      return V[N];
+    }
+  };
+
+  template<typename T>
+  struct GridElement
+  {
+    std::list<T> eList;
+  };
+
+  template<typename T, typename CType = double>
+  struct GridDesc
+  {
+    GridElement<T> *grid;
+    uint32_t gridRes;
+    CType boxSize;
+  };
+
+  template<typename T, int N, typename CType = double>
+  class SparseIterator
+  {
+  protected:
+    typedef CType coord[N];
+
+    int32_t index_min[N], index_max[N], index_cur[N];
+    GridElement<T> *eCur;
+    GridDesc<T> desc;
+    typename std::list<T>::iterator eIterator;
+
+    GridElement<T> *getElement();
+  public:
+    SparseIterator(GridDesc<T>& d, coord& c, CType R);
+    SparseIterator();
+    ~SparseIterator();
+    
+    SparseIterator& operator++();
+
+    bool last();
+
+    T& operator*() { return *eIterator; }
+    T *operator->() { return eIterator.operator->(); }
+  };
+
+  template<typename T, int N, typename CType = double,
+	   typename CAccessor = bracketAccessor<T> >
+  class SparseGrid
+  {
+  public:
+    typedef CType coord[N];
+    typedef SparseIterator<T, N, CType> iterator;
+    CAccessor acc;
+
+    SparseGrid(CType L, uint32_t gridRes);
+    ~SparseGrid();
+    
+    void addElementInGrid(T& e);
+
+    iterator locateElements(coord& c, CType R);
+    void clear();
+
+  protected:
+    GridDesc<T> grid;
+  };
+
+};
+
+#include "sparseGrid.tcc"
+
+#endif

src/sparseGrid.tcc

@@ -0,0 +1,145 @@
+#include <cassert>
+#include <cstring>
+
+namespace CosmoTool {
+   
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseGrid<T,N,CType,CAccessor>::SparseGrid(CType L, uint32_t gridRes)
+   {
+     uint32_t Ntot = 1;
+     for (int i = 0; i < N; i++)
+       Ntot *= gridRes;
+
+     grid.grid = new GridElement<T>[Ntot];
+     grid.boxSize = L;
+     grid.gridRes = gridRes;
+   }
+
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseGrid<T,N,CType,CAccessor>::~SparseGrid()
+   {
+     delete[] grid.grid;
+   }
+    
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   void SparseGrid<T,N,CType,CAccessor>::addElementInGrid(T& e)
+   {
+     uint32_t index[N];
+     uint32_t idx = 0;
+
+     for (int i = N-1; i >= 0; i--)
+       {
+	 CType x = acc(e, i) * grid.gridRes / grid.boxSize;
+	 
+	 index[i] = (int32_t)std::floor((double)x);
+	 assert(index[i] >= 0);
+	 assert(index[i] < grid.gridRes);
+	 idx = (idx*grid.gridRes) + index[i];
+       }
+     assert(idx < grid.gridRes*grid.gridRes*grid.gridRes);
+
+     grid.grid[idx].eList.push_back(e);
+   }
+
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseIterator<T,N,CType>
+      SparseGrid<T,N,CType,CAccessor>::locateElements(coord& c, CType R)
+   {
+     return SparseIterator<T,N,CType>(grid, c, R);
+   }
+
+  template<typename T, int N, typename CType,
+	   typename CAccessor>
+  void SparseGrid<T,N,CType,CAccessor>::clear()
+  {
+    uint32_t N3 = grid.gridRes*grid.gridRes*grid.gridRes;
+
+    for (uint32_t i = 0; i < N3; i++)
+      grid.grid[i].eList.clear();
+  }
+
+
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>::SparseIterator(GridDesc<T>& d, coord& c, CType R)
+  {
+    desc = d;
+    
+    for (uint32_t i = 0; i < N; i++)
+      {
+	CType x_min = (c[i] - R) * d.gridRes / d.boxSize;
+	CType x_max = (c[i] + R) * d.gridRes / d.boxSize;
+	
+	index_min[i] = (int32_t) std::floor((double)x_min);
+	index_max[i] = (int32_t) std::ceil((double)x_max) + 1;
+	index_cur[i] = index_min[i];
+      }
+
+    eCur = getElement();
+    eIterator = eCur->eList.begin();
+    if (eIterator == eCur->eList.end())
+      operator++();
+  }
+  
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>::~SparseIterator()
+  {
+  }
+  
+  template<typename T, int N, typename CType>
+  bool SparseIterator<T,N,CType>::last()
+  {
+    return (index_cur[N-1] == index_max[N-1]);
+  }
+
+ 
+  template<typename T, int N, typename CType>
+  GridElement<T> *SparseIterator<T,N,CType>::getElement()
+  {
+    uint32_t idx = 0;
+    for (int i = (N-1); i >= 0; i--)
+      {
+	int32_t k = (index_cur[i] + desc.gridRes) % desc.gridRes;	
+	idx = (idx*desc.gridRes) + k;
+      }
+    assert(idx < desc.gridRes*desc.gridRes*desc.gridRes);
+    
+    return &desc.grid[idx];
+  }
+
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>& SparseIterator<T,N,CType>::operator++()
+  {
+    if (last())
+      return *this;
+    ++eIterator;
+    if (eIterator != eCur->eList.end())
+      return *this;
+      
+    do
+      {		
+	index_cur[0]++;
+	for (int i = 0; i < (N-1); i++)
+	  if (index_cur[i] == index_max[i])
+	    {
+	      index_cur[i] = index_min[i];
+	      index_cur[i+1]++;
+	    }
+	
+	if (last())
+	  return *this;
+	
+	eCur = getElement();
+	eIterator = eCur->eList.begin();
+      }
+    while (eIterator == eCur->eList.end());
+
+    return *this;
+  }
+ 
+
+
+};