diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index b1ee549..ae62f54 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -27,11 +27,17 @@ namespace CosmoTool typedef boost::shared_ptr FourierMapPtr; typedef boost::shared_ptr > SpectrumFunctionPtr; + virtual ~SpectrumFunction() {} + virtual void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0; virtual SpectrumFunctionPtr copy() const = 0; + virtual const T *data() const { return 0; } + virtual const T *dof() const { return 0; } + virtual long size() const { return -1; } + virtual void sqrt() = 0; virtual void mul(FourierMapType& m) const = 0; @@ -147,10 +153,14 @@ namespace CosmoTool template class MapUtilityFunction { + public: typedef SpectrumFunction Spectrum; - typedef Spectrum *Spectrum_ptr; + typedef boost::shared_ptr Spectrum_ptr; typedef FourierMap > FMap; - typedef Spectrum_ptr (*SpectrumFromMap)(const FMap& m); + + virtual Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const = 0; + virtual Spectrum_ptr newSpectrumFromRaw(T *data, long size, + Spectrum_ptr like_spec) const = 0; }; }; diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp new file mode 100644 index 0000000..c890471 --- /dev/null +++ b/src/fourier/details/euclidian_maps.hpp @@ -0,0 +1,188 @@ +#ifndef __DETAILS_EUCLIDIAN_MAPS +#define __DETAILS_EUCLIDIAN_MAPS + + +namespace CosmoTool +{ + + template + class EuclidianFourierMapBase: public FourierMap + { + public: + typedef std::vector DimArray; + private: + boost::shared_ptr m_data; + DimArray m_dims; + long m_size; + public: + + EuclidianFourierMapBase(boost::shared_ptr indata, const DimArray& indims) + { + m_data = indata; + m_dims = indims; + m_size = 1; + for (int i = 0; i < m_dims.size(); i++) + m_size *= m_dims[i]; + } + + virtual ~EuclidianFourierMapBase() + { + } + + const DimArray& getDims() const { return m_dims; } + + virtual const T *data() const { return m_data.get(); } + virtual T *data() { return m_data.get(); } + virtual long size() const { return m_size; } + + virtual FourierMap *copy() const + { + FourierMap *m = this->mimick(); + m->eigen() = this->eigen(); + return m; + } + + }; + + template + class EuclidianFourierMapReal: public EuclidianFourierMapBase + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + + EuclidianFourierMapReal(boost::shared_ptr indata, const DimArray& indims) + : EuclidianFourierMapBase(indata, indims) + {} + + virtual FourierMap *mimick() const + { + return new EuclidianFourierMapReal( + boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), + std::ptr_fun(fftw_free)), + this->getDims()); + } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapReal& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + return (this->eigen()*m2.eigen()).sum(); + } + }; + + template + class EuclidianFourierMapComplex: public EuclidianFourierMapBase > + { + protected: + typedef boost::shared_ptr > ptr_t; + std::vector delta_k; + int m_dim0; + bool even0, alleven; + long plane_size; + public: + typedef typename EuclidianFourierMapBase >::DimArray DimArray; + + EuclidianFourierMapComplex(ptr_t indata, + int dim0, + const DimArray& indims, + const std::vector& dk) + : 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); + } + } + + virtual FourierMap > *mimick() const + { + return + new EuclidianFourierMapComplex( + ptr_t((std::complex *) + fftw_malloc(sizeof(std::complex)*this->size()), + std::ptr_fun(fftw_free)), + m_dim0, + this->getDims(), + this->delta_k); + } + + template + double get_K(const Array& ik) + { + const DimArray& dims = this->getDims(); + assert(ik.size() == dims.size()); + double k2 = 0; + k2 += CosmoTool::square(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]; + + k2 += CosmoTool::square(delta_k[q]*dk); + } + return std::sqrt(k2); + } + + double get_K(long p) + { + const DimArray& dims = this->getDims(); + DimArray d(delta_k.size()); + for (int q = 0; q < d.size(); q++) + { + d[q] = p%dims[q]; + p = (p-d[q])/dims[q]; + } + return get_K(d); + } + + bool allDimensionsEven() const { return alleven; } + bool firstDimensionEven() const { return even0; } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + const std::complex *d1 = this->data(); + const std::complex *d2 = m2.data(); + const DimArray& dims = this->getDims(); + 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 += 2*(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 += conj(d1[q0]) * d2[q0]; + result += conj(d1[q1]) * d2[q1]; + } + } + return result; + } + + }; +}; + +#endif diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp new file mode 100644 index 0000000..b12864b --- /dev/null +++ b/src/fourier/details/euclidian_spectrum_1d.hpp @@ -0,0 +1,172 @@ +#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D +#define __DETAILS_EUCLIDIAN_SPECTRUM_1D + + +namespace CosmoTool +{ + template + class EuclidianOperator + { + public: + typedef boost::function1 Function; + + Function base, op; + T operator()(T k) { + return op(base(k)); + } + }; + + template + class EuclidianSpectrum_1D: public SpectrumFunction + { + public: + typedef boost::function1 Function; + protected: + Function f; + + static T msqrt(T a) { return std::sqrt(a); } + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr ptr_map; + + EuclidianSpectrum_1D(Function P) + : f(P) + { + } + + void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const; + + SpectrumFunctionPtr copy() const { + return SpectrumFunctionPtr(new EuclidianSpectrum_1D(f)); + } + + void sqrt() { + EuclidianOperator o; + o.base = f; + o.op = &EuclidianSpectrum_1D::msqrt; + f = (Function(o)); + } + + void mul(FourierMapType& m) const; + void mul_sqrt(FourierMapType& m) const; + void mul_inv(FourierMapType& m) const; + void mul_inv_sqrt(FourierMapType& m) const; + }; + + + template + void EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + typedef EuclidianFourierMapComplex MapT; + typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; + typedef typename MapT::DimArray DimArray; + + MapT& rand_map = dynamic_cast(out_map); + + std::complex *d = rand_map.data(); + long idx; + const DimArray& dims = rand_map.getDims(); + long plane_size; + bool alleven = rand_map.allDimensionsEven(); + + for (long p = 1; p < rand_map.size(); p++) + { + double A_k = std::sqrt(0.5*f(rand_map.get_K(p))); + d[p] = std::complex(gsl_ran_gaussian(rng, A_k), + gsl_ran_gaussian(rng, A_k)); + } + // Generate the mean value + d[0] = std::complex(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0); + + if (!rand_map.firstDimensionEven()) + return; + + // Correct the Nyquist plane + idx = dims[0]-1; // Stick to the last element of the first dimension + d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); + // 1D is special case + if (dims.size() == 1) + return; + + plane_size = 1; + for (int q = 1; q < dims.size(); q++) + { + plane_size *= dims[q]; + } + + for (long p = 1; p < plane_size/2; p++) + { + long q = (p+1)*dims[0]-1; + long q2 = (plane_size-p+1)*dims[0]-1; + assert(q < plane_size*dims[0]); + assert(q2 < plane_size*dims[0]); + d[q] = conj(d[q2]); + } + + if (alleven) + { + long q = 0; + for (int i = dims.size()-1; i >= 1; i--) + q = dims[i]*q + dims[i]/2; + q += dims[0]-1; + d[q] = std::complex(d[q].real()+d[q].imag(),0); + } + } + + template + void EuclidianSpectrum_1D::mul(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= f(m_c.get_K(p)); + } + + template + void EuclidianSpectrum_1D::mul_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= std::sqrt(f(m_c.get_K(p))); + } + + template + void EuclidianSpectrum_1D::mul_inv(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + { + T A = f(m_c.get_K(p)); + if (A==0) + d[p] = 0; + else + d[p] /= A; + } + } + + template + void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + { + T A = std::sqrt(f(m_c.get_K(p))); + if (A == 0) + d[p] = 0; + else + d[p] /= A; + } + } + +}; + + +#endif diff --git a/src/fourier/details/euclidian_spectrum_1d_bin.hpp b/src/fourier/details/euclidian_spectrum_1d_bin.hpp new file mode 100644 index 0000000..ff1c98d --- /dev/null +++ b/src/fourier/details/euclidian_spectrum_1d_bin.hpp @@ -0,0 +1,65 @@ +#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN +#define __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN + +#include +#include + +namespace CosmoTool +{ + template + class EuclidianSpectrum_1D_Binned: public EuclidianSpectrum_1D + { + protected: + T *m_data; + long m_size; + T m_kmin, m_kmax, m_logdeltak; + std::vector m_dof; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr ptr_map; + + T interpolate_spectrum(T k) + { + T q = std::log(k/m_kmin)/m_logdeltak; + long ik = std::floor(q); + + if (ik >= m_size-1) + return m_data[m_size-1]; + else if (ik < 0) + return k/m_kmin*m_data[0]; + + return std::exp((q-ik)*m_data[ik+1] + (1-(q-ik))*m_data[ik]); + } + + EuclidianSpectrum_1D_Binned(int numBin, T kmin, T kmax) + : EuclidianSpectrum_1D(boost::bind(&EuclidianSpectrum_1D_Binned::interpolate_spectrum, this, _1)) + { + m_data = new T[numBin]; + m_kmin = kmin; + m_kmax = kmax; + m_size = numBin; + m_logdeltak = std::log(m_kmax/m_kmin); + } + + SpectrumFunctionPtr copy() const + { + EuclidianSpectrum_1D_Binned *s = new EuclidianSpectrum_1D_Binned(m_size, m_kmin, m_kmax); + std::copy(m_data, m_data+m_size, s->m_data); + return SpectrumFunctionPtr(s); + } + + void set_dof(std::vector& dof_array) + { + assert(m_size == dof_array.size()); + m_dof = dof_array; + } + + const T *data() const { return m_data; } + long size() const { return m_size; } + const T *dof() const { return &m_dof[0]; } + }; + +}; + +#endif diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp new file mode 100644 index 0000000..b533227 --- /dev/null +++ b/src/fourier/details/euclidian_transform.hpp @@ -0,0 +1,181 @@ +#ifndef __DETAILS_EUCLIDIAN_TRANSFORM +#define __DETAILS_EUCLIDIAN_TRANSFORM + +namespace CosmoTool +{ + + template + class EuclidianFourierTransform: public FourierTransform + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + private: + typedef FFTW_Calls calls; + EuclidianFourierMapReal *realMap; + EuclidianFourierMapComplex *fourierMap; + typename calls::plan_type m_analysis, m_synthesis; + double volume; + long N, Nc; + DimArray m_dims, m_dims_hc; + std::vector m_L; + public: + EuclidianFourierTransform(const DimArray& dims, const std::vector& L) + { + assert(L.size() == dims.size()); + std::vector dk(L.size()); + std::vector swapped_dims(dims.size()); + + m_dims = dims; + m_dims_hc = dims; + m_dims_hc[0] = dims[0]/2+1; + m_L = L; + + N = 1; + Nc = 1; + volume = 1; + for (int i = 0; i < dims.size(); i++) + { + N *= dims[i]; + Nc *= m_dims_hc[i]; + volume *= L[i]; + dk[i] = 2*M_PI/L[i]; + swapped_dims[dims.size()-1-i] = dims[i]; + } + + realMap = new EuclidianFourierMapReal( + boost::shared_ptr(calls::alloc_real(N), + std::ptr_fun(calls::free)), + m_dims); + fourierMap = new EuclidianFourierMapComplex( + 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); + } + + virtual ~EuclidianFourierTransform() + { + delete realMap; + delete fourierMap; + calls::destroy_plan(m_synthesis); + calls::destroy_plan(m_analysis); + } + + void synthesis() + { + calls::execute(m_synthesis); + realMap->scale(1/volume); + } + + void analysis() + { + calls::execute(m_analysis); + fourierMap->scale(volume/N); + } + + void synthesis_conjugate() + { + calls::execute(m_analysis); + fourierMap->scale(1/volume); + } + + void analysis_conjugate() + { + calls::execute(m_synthesis); + realMap->scale(volume/N); + } + + const FourierMap >& fourierSpace() const + { + return *fourierMap; + } + + FourierMap >& fourierSpace() + { + return *fourierMap; + } + + const FourierMap& realSpace() const + { + return *realMap; + } + + FourierMap& realSpace() + { + return *realMap; + } + + FourierTransform *mimick() const + { + return new EuclidianFourierTransform(m_dims, m_L); + } + }; + + template + class EuclidianFourierTransform_1d: public EuclidianFourierTransform + { + private: + template + static std::vector make_1d_vector(T2 a) + { + T2 arr[2] = { a}; + return std::vector(&arr[0],&arr[1]); + } + public: + EuclidianFourierTransform_1d(int Nx, double Lx) + : EuclidianFourierTransform(make_1d_vector(Nx), make_1d_vector(Lx)) + { + } + + virtual ~EuclidianFourierTransform_1d() {} + + }; + + template + class EuclidianFourierTransform_2d: public EuclidianFourierTransform + { + private: + template + static std::vector make_2d_vector(T2 a, T2 b) + { + T2 arr[2] = { a, b}; + return std::vector(&arr[0], &arr[2]); + } + public: + EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) + : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) + { + } + + virtual ~EuclidianFourierTransform_2d() {} + + }; + + template + class EuclidianFourierTransform_3d: public EuclidianFourierTransform + { + private: + template + static std::vector make_3d_vector(T2 a, T2 b, T2 c) + { + T2 arr[2] = { a, b, c}; + return std::vector(&arr[0], &arr[3]); + } + + public: + EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) + : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) + { + } + + virtual ~EuclidianFourierTransform_3d() {} + }; + +}; + +#endif diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index 1a8c53f..45ec3dc 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -9,522 +9,9 @@ #include "base_types.hpp" #include "fft/fftw_calls.hpp" #include "../algo.hpp" - -namespace CosmoTool -{ - template - class EuclidianOperator - { - public: - typedef boost::function1 Function; - - Function base, op; - T operator()(T k) { - return op(base(k)); - } - }; - - template - class EuclidianSpectrum_1D: public SpectrumFunction - { - public: - typedef boost::function1 Function; - protected: - Function f; - - static T msqrt(T a) { return std::sqrt(a); } - public: - typedef typename SpectrumFunction::FourierMapType FourierMapType; - typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; - typedef boost::shared_ptr ptr_map; - - EuclidianSpectrum_1D(Function P) - : f(P) - { - } - - void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const; - - SpectrumFunctionPtr copy() const { - return SpectrumFunctionPtr(new EuclidianSpectrum_1D(f)); - } - - void sqrt() { - EuclidianOperator o; - o.base = f; - o.op = &EuclidianSpectrum_1D::msqrt; - f = (Function(o)); - } - - void mul(FourierMapType& m) const; - void mul_sqrt(FourierMapType& m) const; - void mul_inv(FourierMapType& m) const; - void mul_inv_sqrt(FourierMapType& m) const; - }; - - template - class EuclidianFourierMapBase: public FourierMap - { - public: - typedef std::vector DimArray; - private: - boost::shared_ptr m_data; - DimArray m_dims; - long m_size; - public: - - EuclidianFourierMapBase(boost::shared_ptr indata, const DimArray& indims) - { - m_data = indata; - m_dims = indims; - m_size = 1; - for (int i = 0; i < m_dims.size(); i++) - m_size *= m_dims[i]; - } - - virtual ~EuclidianFourierMapBase() - { - } - - const DimArray& getDims() const { return m_dims; } - - virtual const T *data() const { return m_data.get(); } - virtual T *data() { return m_data.get(); } - virtual long size() const { return m_size; } - - virtual FourierMap *copy() const - { - FourierMap *m = this->mimick(); - m->eigen() = this->eigen(); - return m; - } - - }; - - template - class EuclidianFourierMapReal: public EuclidianFourierMapBase - { - public: - typedef typename EuclidianFourierMapBase::DimArray DimArray; - - EuclidianFourierMapReal(boost::shared_ptr indata, const DimArray& indims) - : EuclidianFourierMapBase(indata, indims) - {} - - virtual FourierMap *mimick() const - { - return new EuclidianFourierMapReal( - boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), - std::ptr_fun(fftw_free)), - this->getDims()); - } - - virtual T dot_product(const FourierMap& other) const - throw(std::bad_cast) - { - const EuclidianFourierMapReal& m2 = dynamic_cast&>(other); - if (this->size() != m2.size()) - throw std::bad_cast(); - - return (this->eigen()*m2.eigen()).sum(); - } - }; - - template - class EuclidianFourierMapComplex: public EuclidianFourierMapBase > - { - protected: - typedef boost::shared_ptr > ptr_t; - std::vector delta_k; - int m_dim0; - bool even0, alleven; - long plane_size; - public: - typedef typename EuclidianFourierMapBase >::DimArray DimArray; - - EuclidianFourierMapComplex(ptr_t indata, - int dim0, - const DimArray& indims, - const std::vector& dk) - : 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); - } - } - - virtual FourierMap > *mimick() const - { - return - new EuclidianFourierMapComplex( - ptr_t((std::complex *) - fftw_malloc(sizeof(std::complex)*this->size()), - std::ptr_fun(fftw_free)), - m_dim0, - this->getDims(), - this->delta_k); - } - - template - double get_K(const Array& ik) - { - const DimArray& dims = this->getDims(); - assert(ik.size() == dims.size()); - double k2 = 0; - k2 += CosmoTool::square(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]; - - k2 += CosmoTool::square(delta_k[q]*dk); - } - return std::sqrt(k2); - } - - double get_K(long p) - { - const DimArray& dims = this->getDims(); - DimArray d(delta_k.size()); - for (int q = 0; q < d.size(); q++) - { - d[q] = p%dims[q]; - p = (p-d[q])/dims[q]; - } - return get_K(d); - } - - bool allDimensionsEven() const { return alleven; } - bool firstDimensionEven() const { return even0; } - - virtual std::complex dot_product(const FourierMap >& other) const - throw(std::bad_cast) - { - const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); - if (this->size() != m2.size()) - throw std::bad_cast(); - - const std::complex *d1 = this->data(); - const std::complex *d2 = m2.data(); - const DimArray& dims = this->getDims(); - 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 += 2*(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 += conj(d1[q0]) * d2[q0]; - result += conj(d1[q1]) * d2[q1]; - } - } - return result; - } - - }; - - template - class EuclidianFourierTransform: public FourierTransform - { - public: - typedef typename EuclidianFourierMapBase::DimArray DimArray; - private: - typedef FFTW_Calls calls; - EuclidianFourierMapReal *realMap; - EuclidianFourierMapComplex *fourierMap; - typename calls::plan_type m_analysis, m_synthesis; - double volume; - long N, Nc; - DimArray m_dims, m_dims_hc; - std::vector m_L; - public: - EuclidianFourierTransform(const DimArray& dims, const std::vector& L) - { - assert(L.size() == dims.size()); - std::vector dk(L.size()); - std::vector swapped_dims(dims.size()); - - m_dims = dims; - m_dims_hc = dims; - m_dims_hc[0] = dims[0]/2+1; - m_L = L; - - N = 1; - Nc = 1; - volume = 1; - for (int i = 0; i < dims.size(); i++) - { - N *= dims[i]; - Nc *= m_dims_hc[i]; - volume *= L[i]; - dk[i] = 2*M_PI/L[i]; - swapped_dims[dims.size()-1-i] = dims[i]; - } - - realMap = new EuclidianFourierMapReal( - boost::shared_ptr(calls::alloc_real(N), - std::ptr_fun(calls::free)), - m_dims); - fourierMap = new EuclidianFourierMapComplex( - 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); - } - - virtual ~EuclidianFourierTransform() - { - delete realMap; - delete fourierMap; - calls::destroy_plan(m_synthesis); - calls::destroy_plan(m_analysis); - } - - void synthesis() - { - calls::execute(m_synthesis); - realMap->scale(1/volume); - } - - void analysis() - { - calls::execute(m_analysis); - fourierMap->scale(volume/N); - } - - void synthesis_conjugate() - { - calls::execute(m_analysis); - fourierMap->scale(1/volume); - } - - void analysis_conjugate() - { - calls::execute(m_synthesis); - realMap->scale(volume/N); - } - - const FourierMap >& fourierSpace() const - { - return *fourierMap; - } - - FourierMap >& fourierSpace() - { - return *fourierMap; - } - - const FourierMap& realSpace() const - { - return *realMap; - } - - FourierMap& realSpace() - { - return *realMap; - } - - FourierTransform *mimick() const - { - return new EuclidianFourierTransform(m_dims, m_L); - } - }; - - template - class EuclidianFourierTransform_1d: public EuclidianFourierTransform - { - private: - template - static std::vector make_1d_vector(T2 a) - { - T2 arr[2] = { a}; - return std::vector(&arr[0],&arr[1]); - } - public: - EuclidianFourierTransform_1d(int Nx, double Lx) - : EuclidianFourierTransform(make_1d_vector(Nx), make_1d_vector(Lx)) - { - } - - virtual ~EuclidianFourierTransform_1d() {} - - }; - - template - class EuclidianFourierTransform_2d: public EuclidianFourierTransform - { - private: - template - static std::vector make_2d_vector(T2 a, T2 b) - { - T2 arr[2] = { a, b}; - return std::vector(&arr[0], &arr[2]); - } - public: - EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) - : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) - { - } - - virtual ~EuclidianFourierTransform_2d() {} - - }; - - template - class EuclidianFourierTransform_3d: public EuclidianFourierTransform - { - private: - template - static std::vector make_3d_vector(T2 a, T2 b, T2 c) - { - T2 arr[2] = { a, b, c}; - return std::vector(&arr[0], &arr[3]); - } - - public: - EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) - : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) - { - } - - virtual ~EuclidianFourierTransform_3d() {} - }; - - - template - void EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const - { - typedef EuclidianFourierMapComplex MapT; - typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; - typedef typename MapT::DimArray DimArray; - - MapT& rand_map = dynamic_cast(out_map); - - std::complex *d = rand_map.data(); - long idx; - const DimArray& dims = rand_map.getDims(); - long plane_size; - bool alleven = rand_map.allDimensionsEven(); - - for (long p = 1; p < rand_map.size(); p++) - { - double A_k = std::sqrt(0.5*f(rand_map.get_K(p))); - d[p] = std::complex(gsl_ran_gaussian(rng, A_k), - gsl_ran_gaussian(rng, A_k)); - } - // Generate the mean value - d[0] = std::complex(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0); - - if (!rand_map.firstDimensionEven()) - return; - - // Correct the Nyquist plane - idx = dims[0]-1; // Stick to the last element of the first dimension - d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); - // 1D is special case - if (dims.size() == 1) - return; - - plane_size = 1; - for (int q = 1; q < dims.size(); q++) - { - plane_size *= dims[q]; - } - - for (long p = 1; p < plane_size/2; p++) - { - long q = (p+1)*dims[0]-1; - long q2 = (plane_size-p+1)*dims[0]-1; - assert(q < plane_size*dims[0]); - assert(q2 < plane_size*dims[0]); - d[q] = conj(d[q2]); - } - - if (alleven) - { - long q = 0; - for (int i = dims.size()-1; i >= 1; i--) - q = dims[i]*q + dims[i]/2; - q += dims[0]-1; - d[q] = std::complex(d[q].real()+d[q].imag(),0); - } - } - - template - void EuclidianSpectrum_1D::mul(FourierMapType& m) const - { - EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); - std::complex *d = m.data(); - - for (long p = 0; p < m_c.size(); p++) - d[p] *= f(m_c.get_K(p)); - } - - template - void EuclidianSpectrum_1D::mul_sqrt(FourierMapType& m) const - { - EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); - std::complex *d = m.data(); - - for (long p = 0; p < m_c.size(); p++) - d[p] *= std::sqrt(f(m_c.get_K(p))); - } - - template - void EuclidianSpectrum_1D::mul_inv(FourierMapType& m) const - { - EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); - std::complex *d = m.data(); - - for (long p = 0; p < m_c.size(); p++) - { - T A = f(m_c.get_K(p)); - if (A==0) - d[p] = 0; - else - d[p] /= A; - } - } - - template - void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMapType& m) const - { - EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); - std::complex *d = m.data(); - - for (long p = 0; p < m_c.size(); p++) - { - T A = std::sqrt(f(m_c.get_K(p))); - if (A == 0) - d[p] = 0; - else - d[p] /= A; - } - } - -}; +#include "details/euclidian_maps.hpp" +#include "details/euclidian_spectrum_1d.hpp" +#include "details/euclidian_spectrum_1d_bin.hpp" +#include "details/euclidian_transform.hpp" #endif