From 8af3dc4359a1ce7d0eb315332353d239315c073a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 11 Feb 2015 15:10:09 +0100 Subject: [PATCH 1/6] Made fourier library more flexible. Compatible with boost::multi_array now. --- src/fourier/details/euclidian_maps.hpp | 126 +++++++++++++------- src/fourier/details/euclidian_transform.hpp | 20 ++-- src/fourier/fft/fftw_calls.hpp | 2 + 3 files changed, 99 insertions(+), 49 deletions(-) diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp index 5b808b2..1e82a90 100644 --- a/src/fourier/details/euclidian_maps.hpp +++ b/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 +#include namespace CosmoTool { + namespace details { + static void no_free_euclidian_map(void *) {} + } + template class EuclidianFourierMapBase: public FourierMap { @@ -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 + EuclidianFourierMapBase(ArrayType& indata) + { + m_data = boost::shared_ptr( + 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& array() { + return boost::multi_array_ref(m_data.get(), boost::extents[size]); + } + + boost::const_multi_array_ref& array() const { + return boost::const_multi_array_ref(m_data.get(), boost::extents[size]); + } + virtual FourierMap *copy() const { FourierMap *m = this->mimick(); @@ -90,6 +116,12 @@ namespace CosmoTool : EuclidianFourierMapBase(indata, indims) {} + template + EuclidianFourierMapReal(ArrayType& indata) + : EuclidianFourierMapBase(indata) + {} + + virtual FourierMap *mimick() const { return new EuclidianFourierMapReal( @@ -125,16 +157,32 @@ namespace CosmoTool int dim0, const DimArray& indims, const std::vector& dk) - : EuclidianFourierMapBase >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0) + : EuclidianFourierMapBase >(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++) - { - plane_size *= indims[q]; - alleven = alleven && ((indims[q]%2)==0); - } + for (int q = 1; q < indims.size(); q++) { + plane_size *= indims[q]; + alleven = alleven && ((indims[q]%2)==0); + } + } + + template + EuclidianFourierMapComplex(ArrayType& indata, + int dim0, + const std::vector& dk) + : EuclidianFourierMapBase >(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 < indims.size(); q++) { + plane_size *= m_dims[q]; + alleven = alleven && ((indims[q]%2)==0); + } } virtual FourierMap > *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++) - { - int dk = ik[q]; - if (dk > dims[q]/2) - dk = dk - dims[q]; + for (int q = 1; q < ik.size(); q++) { + int dk = ik[q]; + if (dk > dims[q]/2) + dk = dk - dims[q]; - kvec[q] = dk*delta_k[q]; - } + kvec[q] = dk*delta_k[q]; + } } template @@ -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++) - { - int dk = ik[q]; + for (int q = 1; q < ik.size(); q++) { + int dk = ik[q]; - if (dk > dims[q]/2) - dk = dk - dims[q]; - - k2 += CosmoTool::square(delta_k[q]*dk); - } + if (dk > dims[q]/2) + dk = dk - dims[q]; + + k2 += CosmoTool::square(delta_k[q]*dk); + } return std::sqrt(k2); } @@ -231,7 +277,8 @@ namespace CosmoTool virtual std::complex dot_product(const FourierMap >& other) const throw(std::bad_cast) { - const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); + const EuclidianFourierMapComplex& m2 = + dynamic_cast&>(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 result = 0; - for (long q0 = 1; q0 < N0-1; q0++) - { - for (long p = 0; p < plane_size; p++) - { - long idx = q0+dims[0]*p; - assert(idx < this->size()); - result += T(2)*(std::conj(d1[idx]) * d2[idx]).real(); - } - } - if (even0) - { - for (long p = 0; p < plane_size; p++) - { - 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]; - } - } + for (long q0 = 1; q0 < N0-1; q0++) { + for (long p = 0; p < plane_size; p++) { + long idx = q0+dims[0]*p; + assert(idx < this->size()); + result += T(2)*(std::conj(d1[idx]) * d2[idx]).real(); + } + } + if (even0) { + for (long p = 0; p < plane_size; p++) + { + 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; } diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp index ac16fa9..4754361 100644 --- a/src/fourier/details/euclidian_transform.hpp +++ b/src/fourier/details/euclidian_transform.hpp @@ -96,16 +96,20 @@ namespace CosmoTool std::ptr_fun(calls::free)), m_dims); fourierMap = new EuclidianFourierMapComplex( - boost::shared_ptr >((std::complex*)calls::alloc_complex(Nc), - std::ptr_fun(calls::free)), + boost::shared_ptr >( + (std::complex*)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], - realMap->data(), (typename calls::complex_type *)fourierMap->data(), - 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); + m_analysis = calls::plan_dft_r2c( + dims.size(), &swapped_dims[0], + realMap->data(), + (typename calls::complex_type *)fourierMap->data(), + 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); } } diff --git a/src/fourier/fft/fftw_calls.hpp b/src/fourier/fft/fftw_calls.hpp index 7d1138f..b5f37db 100644 --- a/src/fourier/fft/fftw_calls.hpp +++ b/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) \ From 1dbe46a3587324fdc6276de32195e6b89a2fe01f Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 11 Feb 2015 18:23:28 +0100 Subject: [PATCH 2/6] Added missing FFTW call --- src/fourier/fft/fftw_calls.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/fourier/fft/fftw_calls.hpp b/src/fourier/fft/fftw_calls.hpp index b5f37db..f61f94b 100644 --- a/src/fourier/fft/fftw_calls.hpp +++ b/src/fourier/fft/fftw_calls.hpp @@ -90,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) \ { \ From ebaf13e33670e5b049cc7c4df68f968bc9b597f2 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 11 Feb 2015 21:23:44 +0100 Subject: [PATCH 3/6] HDF5 work --- src/hdf5_array.hpp | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp index abcdc79..68c3dcb 100644 --- a/src/hdf5_array.hpp +++ b/src/hdf5_array.hpp @@ -97,13 +97,13 @@ namespace CosmoTool { //! \author leo Goodstadt (04 March 2013) //! //!_______________________________________________________________________________________ - template + template void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, - const boost::multi_array& data, + const ArrayType& data, const hdf5_data_type& datatype) { - std::vector dimensions(data.shape(), data.shape() + DIMENSIONS); - H5::DataSpace dataspace(DIMENSIONS, dimensions.data()); + std::vector dimensions(data.shape(), data.shape() + data.num_dimensions()); + 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 - void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array& data ) + template<> struct get_hdf5_data_type > { + static H5::DataType type() { + return hdf5_ComplexType::ctype()->type; + } + }; + + template<> struct get_hdf5_data_type > { + static H5::DataType type() { + return hdf5_ComplexType::ctype()->type; + } + }; + + + template + void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const ArrayType& data ) { + typedef typename ArrayType::value_type T; get_hdf5_data_type hdf_data_type; + hdf5_write_array(fg, data_set_name, data, hdf_data_type.type()); } - - template - void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array, DIMENSIONS>& data ) - { - hdf5_write_array(fg, data_set_name, data, hdf5_ComplexType::ctype()->type); - } - // HDF5 array reader // From a26de4a2faf52a12fab2496974d2aafe43674c7a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 11 Feb 2015 21:27:27 +0100 Subject: [PATCH 4/6] Fixed compilation error in HDF5 --- src/fourier/details/euclidian_maps.hpp | 6 +++--- src/hdf5_array.hpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp index 1e82a90..2e6d5bb 100644 --- a/src/fourier/details/euclidian_maps.hpp +++ b/src/fourier/details/euclidian_maps.hpp @@ -179,9 +179,9 @@ namespace CosmoTool assert(dk.size() == indata.num_dimensions()); plane_size = 1; alleven = true; - for (int q = 1; q < indims.size(); q++) { - plane_size *= m_dims[q]; - alleven = alleven && ((indims[q]%2)==0); + for (int q = 1; q < this->m_dims.size(); q++) { + plane_size *= this->m_dims[q]; + alleven = alleven && ((this->m_dims[q]%2)==0); } } diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp index 68c3dcb..98ea28f 100644 --- a/src/hdf5_array.hpp +++ b/src/hdf5_array.hpp @@ -150,7 +150,7 @@ namespace CosmoTool { template void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const ArrayType& data ) { - typedef typename ArrayType::value_type T; + typedef typename ArrayType::element T; get_hdf5_data_type hdf_data_type; hdf5_write_array(fg, data_set_name, data, hdf_data_type.type()); From c853abe354aa35ef809362d538b646f29bb95275 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 11 Feb 2015 21:38:43 +0100 Subject: [PATCH 5/6] Simplified template substitution. Also more general. --- src/hdf5_array.hpp | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp index 98ea28f..c6ce5c9 100644 --- a/src/hdf5_array.hpp +++ b/src/hdf5_array.hpp @@ -181,38 +181,34 @@ namespace CosmoTool { } }; - template + template void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, - boost::multi_array& data, + ArrayType& data, const hdf5_data_type& datatype) { H5::DataSet dataset = fg.openDataSet(data_set_name); H5::DataSpace dataspace = dataset.getSpace(); - std::vector dimensions(DIMENSIONS); + std::vector 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::build(dimensions.data())); + data.resize(hdf5_extent_gen::build(dimensions.data())); dataset.read(data.data(), datatype); } - template - void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array& data ) + template + void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, ArrayType& data ) { + typedef typename ArrayType::element T; + get_hdf5_data_type hdf_data_type; hdf5_read_array(fg, data_set_name, data, hdf_data_type.type()); } - template - void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array, DIMENSIONS>& data ) - { - hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType::ctype()->type); - } - #define CTOOL_HDF5_NAME(STRUCT) BOOST_PP_CAT(hdf5_,STRUCT) #define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \ From 1ab31e6c197b68b7ea4f35fd0feb2fd49a84e3ed Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 12 Feb 2015 08:50:07 +0100 Subject: [PATCH 6/6] hdf5_read_array now supports array witout resize. --- src/hdf5_array.hpp | 44 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp index c6ce5c9..318181f 100644 --- a/src/hdf5_array.hpp +++ b/src/hdf5_array.hpp @@ -180,7 +180,49 @@ namespace CosmoTool { return boost::extents; } }; + + + template class array_has_resize { + struct Fallback { int resize; }; + struct Derived: Array, Fallback {}; + typedef char yes[1]; + typedef char no[2]; + + template struct Check; + + template + static yes& func(Check *); + + template + static no& func(...); + public: + typedef array_has_resize type; + enum { value = sizeof(func(0)) == sizeof(no) }; + }; + + + + template + typename boost::enable_if< + array_has_resize + >::type + hdf5_resize_array(ArrayType& data, std::vector& dims) { + data.resize( + hdf5_extent_gen::build(dims.data()) + ); + } + + template + typename boost::disable_if< + array_has_resize + >::type + hdf5_resize_array(ArrayType& data, std::vector& dims) { + for (long i = 0; i < data.num_dimensions(); i++) { + assert (data.shape()[i] == dims[i]); + } + } + template void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, ArrayType& data, @@ -196,7 +238,7 @@ namespace CosmoTool { } dataspace.getSimpleExtentDims(dimensions.data()); - data.resize(hdf5_extent_gen::build(dimensions.data())); + hdf5_resize_array(data, dimensions); dataset.read(data.data(), datatype); }