Merge branch 'master' of bitbucket.org:glavaux/cosmotool

src/fourier/details/euclidian_maps.hpp

@@ -37,10 +37,15 @@ knowledge of the CeCILL license and that you accept its terms.
 #define __DETAILS_EUCLIDIAN_MAPS
 
 #include <cmath>
+#include <boost/multi_array.hpp>
 
 namespace CosmoTool
 {
   
+  namespace details {
+    static void no_free_euclidian_map(void *) {}
+  }
+  
   template<typename T>
   class EuclidianFourierMapBase: public FourierMap<T>
   {
@@ -58,9 +63,22 @@ namespace CosmoTool
       m_dims = indims;
       m_size = 1;
       for (int i = 0; i < m_dims.size(); i++)
-	m_size *= m_dims[i];
+        m_size *= m_dims[i];
     }
 
+    template<typename ArrayType>
+    EuclidianFourierMapBase(ArrayType& indata)
+    {
+      m_data = boost::shared_ptr<T>(
+            indata.origin(), 
+            std::ptr_fun(details::no_free_euclidian_map));
+      m_dims = DimArray(indata.num_dimensions());
+      m_size = indata.num_elements();
+      for (int i = 0; i < m_dims.size(); i++)
+        m_dims[i] = indata.shape()[i];
+    }
+
+
     virtual ~EuclidianFourierMapBase()
     {
     }
@@ -71,6 +89,14 @@ namespace CosmoTool
     virtual T *data() { return m_data.get(); }
     virtual long size() const  { return m_size; } 
     
+    boost::multi_array_ref<T, 1>& array() { 
+      return boost::multi_array_ref<T, 1>(m_data.get(), boost::extents[size]);
+    }
+
+    boost::const_multi_array_ref<T, 1>& array() const { 
+      return boost::const_multi_array_ref<T, 1>(m_data.get(), boost::extents[size]);
+    }
+    
     virtual FourierMap<T> *copy() const
     {
       FourierMap<T> *m = this->mimick();
@@ -90,6 +116,12 @@ namespace CosmoTool
       : EuclidianFourierMapBase<T>(indata, indims) 
     {}
 
+    template<typename ArrayType>
+    EuclidianFourierMapReal(ArrayType& indata) 
+      : EuclidianFourierMapBase<T>(indata) 
+    {}
+
+
     virtual FourierMap<T> *mimick() const
     {
       return new EuclidianFourierMapReal<T>(
@@ -125,16 +157,32 @@ namespace CosmoTool
 			       int dim0,
 			       const DimArray& indims, 
 			       const std::vector<double>& dk)
-      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
+      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), 
+        delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
     {
       assert(dk.size() == indims.size()); 
       plane_size = 1;
       alleven = true;
-      for (int q = 1; q < indims.size(); q++)
+      for (int q = 1; q < indims.size(); q++) {
-	{
+        plane_size *= indims[q];
-	  plane_size *= indims[q];
+        alleven = alleven && ((indims[q]%2)==0);
-	  alleven = alleven && ((indims[q]%2)==0);
+      }
-	}
+    }
+
+    template<typename ArrayType>
+    EuclidianFourierMapComplex(ArrayType& indata, 
+                      	       int dim0,
+                               const std::vector<double>& dk) 
+      : EuclidianFourierMapBase<std::complex<T> >(indata),
+        delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
+    {
+      assert(dk.size() == indata.num_dimensions()); 
+      plane_size = 1;
+      alleven = true;
+      for (int q = 1; q < this->m_dims.size(); q++) {
+        plane_size *= this->m_dims[q];
+        alleven = alleven && ((this->m_dims[q]%2)==0);
+      }
     }
 
     virtual FourierMap<std::complex<T> > *mimick() const
@@ -162,14 +210,13 @@ namespace CosmoTool
       assert(kvec.size() == dims.size());
 
       kvec[0] = ik[0] * delta_k[0];
-      for (int q = 1; q < ik.size(); q++)
+      for (int q = 1; q < ik.size(); q++) {
-	{
+        int dk = ik[q];
-	  int dk = ik[q];
+        if (dk > dims[q]/2)
-	  if (dk > dims[q]/2)
+          dk = dk - dims[q];
-	    dk = dk - dims[q];
 
-          kvec[q] = dk*delta_k[q];
+        kvec[q] = dk*delta_k[q];
-        }
+      }
     }
 
     template<typename Array2>
@@ -201,15 +248,14 @@ namespace CosmoTool
       double k2 = 0;
       k2 += CosmoTool::square(ik[0]*delta_k[0]);
 
-      for (int q = 1; q < ik.size(); q++)
+      for (int q = 1; q < ik.size(); q++) {
-	{
+        int dk = ik[q];
-	  int dk = ik[q];
 
-	  if (dk > dims[q]/2)
+        if (dk > dims[q]/2)
-	    dk = dk - dims[q];
+          dk = dk - dims[q];
-	  
+        
-	  k2 += CosmoTool::square(delta_k[q]*dk);
+        k2 += CosmoTool::square(delta_k[q]*dk);
-	}
+      }
       return std::sqrt(k2);
     }
 
@@ -231,7 +277,8 @@ namespace CosmoTool
     virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
       throw(std::bad_cast)
     {
-      const EuclidianFourierMapComplex<T>& m2 = dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
+      const EuclidianFourierMapComplex<T>& m2 = 
+        dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
       if (this->size() != m2.size())
         throw std::bad_cast();
       
@@ -241,24 +288,21 @@ namespace CosmoTool
       int N0 = dims[0] + (even0 ? 0 : 1);
       std::complex<T> result = 0;
 
-      for (long q0 = 1; q0 < N0-1; q0++)
+      for (long q0 = 1; q0 < N0-1; q0++) {
-	{
+        for (long p = 0; p < plane_size; p++) {
-	  for (long p = 0; p < plane_size; p++)
+          long idx = q0+dims[0]*p;
-	    {
+          assert(idx < this->size());
-	      long idx = q0+dims[0]*p;
+          result += T(2)*(std::conj(d1[idx]) * d2[idx]).real();
-	      assert(idx < this->size());
+        }
-	      result += T(2)*(std::conj(d1[idx]) * d2[idx]).real();
+      }
-	    }
+      if (even0) {
-	}
+        for (long p = 0; p < plane_size; p++)
-      if (even0)
+          {
-	{
+            long q0 = N0*p, q1 = (p+1)*N0-1;
-	  for (long p = 0; p < plane_size; p++)
+            result += T(2)*std::conj(d1[q0]) * d2[q0];
-	    {
+            result += T(2)*std::conj(d1[q1]) * d2[q1];
-	      long q0 = N0*p, q1 = (p+1)*N0-1;
+          }
-	      result += T(2)*std::conj(d1[q0]) * d2[q0];
+        }
-	      result += T(2)*std::conj(d1[q1]) * d2[q1];
-	    }
-	}
       return result;
     }
 

src/fourier/details/euclidian_transform.hpp

@@ -96,16 +96,20 @@ namespace CosmoTool
                                     std::ptr_fun(calls::free)),
                                     m_dims);
       fourierMap = new EuclidianFourierMapComplex<T>(
-                boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(Nc),
+                boost::shared_ptr<std::complex<T> >(
-                                                    std::ptr_fun(calls::free)), 
+                  (std::complex<T>*)calls::alloc_complex(Nc),
+                  std::ptr_fun(calls::free)), 
                 dims[0], m_dims_hc, dk);
       {
-        m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0],
+        m_analysis = calls::plan_dft_r2c(
-                                        realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+              dims.size(), &swapped_dims[0],
-                                        FFTW_DESTROY_INPUT|FFTW_MEASURE);
+              realMap->data(), 
-        m_synthesis = calls::plan_dft_c2r(dims.size(), &swapped_dims[0],
+              (typename calls::complex_type *)fourierMap->data(),
-                                          (typename calls::complex_type *)fourierMap->data(), realMap->data(),
+              FFTW_DESTROY_INPUT|FFTW_MEASURE);
-                                          FFTW_DESTROY_INPUT|FFTW_MEASURE);
+        m_synthesis = calls::plan_dft_c2r(
+              dims.size(), &swapped_dims[0],
+              (typename calls::complex_type *)fourierMap->data(), realMap->data(),
+              FFTW_DESTROY_INPUT|FFTW_MEASURE);
       }
     }
     

src/fourier/fft/fftw_calls.hpp

@@ -68,6 +68,8 @@ public: \
   static void free(void *p) { fftw_free(p); } \
 \
   static void execute(plan_type p) { prefix ## _execute(p); } \
+  static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _execute_dft_r2c(p, in, out); } \
+  static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _execute_dft_c2r(p, in, out); } \
   static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
 		       real_type *in, complex_type *out, \
 		       unsigned flags) \
@@ -88,6 +90,13 @@ public: \
   { \
     return prefix ## _plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
   } \
+  static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
+                                   complex_type *in, real_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+\
   static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
                                 complex_type *out, unsigned flags) \
   { \

src/hdf5_array.hpp

@@ -97,13 +97,13 @@ namespace CosmoTool {
 //!     \author leo Goodstadt (04 March 2013)
 //!     
 //!_______________________________________________________________________________________
-  template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
+  template<typename ArrayType,  typename hdf5_data_type>
   void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, 
-                        const boost::multi_array<T, DIMENSIONS>& data, 
+                        const ArrayType& data, 
                         const hdf5_data_type& datatype)
   {
-      std::vector<hsize_t> dimensions(data.shape(), data.shape() + DIMENSIONS);
+      std::vector<hsize_t> dimensions(data.shape(), data.shape() + data.num_dimensions());
-      H5::DataSpace dataspace(DIMENSIONS, dimensions.data());
+      H5::DataSpace dataspace(data.num_dimensions(), dimensions.data());
 
       H5::DataSet dataset = fg.createDataSet(data_set_name, datatype, dataspace);
 
@@ -134,19 +134,27 @@ namespace CosmoTool {
     }
   };
 
-  template<typename T, std::size_t DIMENSIONS>
+  template<> struct get_hdf5_data_type<std::complex<float> > { 
-  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<T, DIMENSIONS>& data )
+    static H5::DataType type() { 
+      return hdf5_ComplexType<float>::ctype()->type; 
+    }
+  };
+
+  template<> struct get_hdf5_data_type<std::complex<double> > { 
+    static H5::DataType type() { 
+      return hdf5_ComplexType<double>::ctype()->type; 
+    }
+  };
+
+
+  template<typename ArrayType>
+  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const ArrayType& data )
   {
+      typedef typename ArrayType::element T;
       get_hdf5_data_type<T> hdf_data_type;
+      
       hdf5_write_array(fg, data_set_name, data, hdf_data_type.type());
   }
-
-  template<typename T, std::size_t DIMENSIONS>
-  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<std::complex<T>, DIMENSIONS>& data )
-  {
-      hdf5_write_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
-  }
-
   
   // HDF5 array reader
   //
@@ -172,39 +180,77 @@ namespace CosmoTool {
       return boost::extents;
     }
   };
+  
+  
+  template<typename Array> class array_has_resize {
+    struct Fallback { int resize; };
+    struct Derived: Array, Fallback {};
 
-  template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
+    typedef char yes[1];
+    typedef char no[2];
+
+    template<typename U, U> struct Check;
+
+    template<typename U> 
+    static yes& func(Check<int Fallback::*, &U::resize> *);
+
+    template<typename U>
+    static no& func(...);
+  public:
+    typedef array_has_resize type;
+    enum { value = sizeof(func<Derived>(0)) == sizeof(no) };
+  };
+
+  
+
+  template<typename ArrayType>
+    typename boost::enable_if<
+      array_has_resize<ArrayType> 
+  >::type 
+    hdf5_resize_array(ArrayType& data, std::vector<hsize_t>& dims) {
+    data.resize(
+        hdf5_extent_gen<ArrayType::dimensionality>::build(dims.data())
+    );
+  }
+
+  template<typename ArrayType>
+    typename boost::disable_if<
+      array_has_resize<ArrayType> 
+  >::type 
+    hdf5_resize_array(ArrayType& data, std::vector<hsize_t>& dims) {
+    for (long i = 0; i < data.num_dimensions(); i++) {
+        assert (data.shape()[i] == dims[i]);
+    }
+  }
+  
+  template<typename ArrayType, typename hdf5_data_type>
   void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, 
-                       boost::multi_array<T, DIMENSIONS>& data, 
+                       ArrayType& data, 
                        const hdf5_data_type& datatype)
   {
       H5::DataSet dataset = fg.openDataSet(data_set_name);
       H5::DataSpace dataspace = dataset.getSpace();
-      std::vector<hsize_t> dimensions(DIMENSIONS);
+      std::vector<hsize_t> dimensions(data.num_dimensions());
       
-      if (dataspace.getSimpleExtentNdims() != DIMENSIONS)
+      if (dataspace.getSimpleExtentNdims() != data.num_dimensions())
         {
           throw InvalidDimensions();
         }
             
       dataspace.getSimpleExtentDims(dimensions.data());
-      data.resize(hdf5_extent_gen<DIMENSIONS>::build(dimensions.data()));
+      hdf5_resize_array(data, dimensions);
       dataset.read(data.data(), datatype);
   }
 
-  template<typename T, std::size_t DIMENSIONS>
+  template<typename ArrayType>
-  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<T, DIMENSIONS>& data )
+  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, ArrayType& data )
   {
+      typedef typename ArrayType::element T;
+      
       get_hdf5_data_type<T> hdf_data_type;
       hdf5_read_array(fg, data_set_name, data, hdf_data_type.type());
   }
 
-  template<typename T, std::size_t DIMENSIONS>
-  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<std::complex<T>, DIMENSIONS>& data )
-  {
-      hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
-  }
-
 
 #define CTOOL_HDF5_NAME(STRUCT) BOOST_PP_CAT(hdf5_,STRUCT)
 #define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \