From cc7996175cc1a5989b3a4f16fb6efe2e2108d19f Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 10 Nov 2012 17:35:27 -0500 Subject: [PATCH] Implemented power spectrum support in Euclidian/Healpix --- sample/test_fft_calls.cpp | 21 +++++++++ sample/test_healpix_calls.cpp | 2 +- src/fourier/base_types.hpp | 23 ++++++++-- src/fourier/euclidian.hpp | 66 ++++++++++++++++++++++++---- src/fourier/healpix.hpp | 82 +++++++++++++++++++++++++++++++++-- 5 files changed, 178 insertions(+), 16 deletions(-) diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp index 06c1c9a..e774a53 100644 --- a/sample/test_fft_calls.cpp +++ b/sample/test_fft_calls.cpp @@ -1,13 +1,23 @@ +#include "yorick.hpp" +#include #include +#include #include "fourier/euclidian.hpp" using namespace CosmoTool; using namespace std; +double spectrum_generator(double k) +{ + return 1/(0.1+pow(k, 3.0)); +} + int main() { EuclidianFourierTransform_2d dft(128,128,1.0,1.0); + EuclidianSpectrum_1D spectrum(spectrum_generator); double volume = 128*128; + gsl_rng *rng = gsl_rng_alloc(gsl_rng_default); dft.realSpace().eigen().setRandom(); dft.analysis(); @@ -15,6 +25,17 @@ int main() cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl; dft.synthesis(); cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; + + dft.realSpace().scale(2.0); + dft.fourierSpace().scale(0.2); + + SpectrumFunction::FourierMapPtr m = spectrum.newRandomFourier(rng, dft.fourierSpace()); + dft.fourierSpace() = *m.get(); + dft.synthesis(); + + uint32_t dims[2] = { 128, 128 }; + CosmoTool::saveArray("generated_map.nc", dft.realSpace().data(), dims, 2); + return 0; } diff --git a/sample/test_healpix_calls.cpp b/sample/test_healpix_calls.cpp index ae1844e..d1f1bf3 100644 --- a/sample/test_healpix_calls.cpp +++ b/sample/test_healpix_calls.cpp @@ -6,7 +6,7 @@ using namespace std; int main() { - HealpixFourierTransform dft(128,3*128,3*128, 40); + HealpixFourierTransform dft(128,2*128,2*128, 40); long Npix = dft.realSpace().size(); dft.realSpace().eigen().setRandom(); diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index 9f54320..3a10798 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -24,14 +24,17 @@ namespace CosmoTool typedef Eigen::Map MapType; typedef Eigen::Map ConstMapType; typedef FourierMap > FourierMapType; + typedef boost::shared_ptr FourierMapPtr; - virtual boost::shared_ptr + virtual FourierMapPtr newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; - virtual void mul(FourierMap >& m) const = 0; + virtual void mul(FourierMapType& m) const = 0; + virtual void mul_sqrt(FourierMapType& m) const = 0; + virtual void mul_inv(FourierMapType& m) const = 0; + virtual void mul_inv_sqrt(FourierMapType& m) const = 0; }; - template class FourierMap { @@ -60,6 +63,12 @@ namespace CosmoTool return ConstMapType(data(), size()); } + FourierMap& operator=(const FourierMap& m) + { + eigen() = m.eigen(); + return *this; + } + void sqrt() { MapType m = eigen(); @@ -130,6 +139,14 @@ namespace CosmoTool virtual void synthesis_conjugate() = 0; }; + template + class MapUtilityFunction + { + typedef SpectrumFunction Spectrum; + typedef Spectrum *Spectrum_ptr; + typedef FourierMap > FMap; + typedef Spectrum_ptr (*SpectrumFromMap)(const FMap& m); + }; }; diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index fac27f6..605013c 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -30,6 +30,9 @@ namespace CosmoTool ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; void mul(FourierMap >& m) const; + void mul_sqrt(FourierMap >& m) const; + void mul_inv(FourierMap >& m) const; + void mul_inv_sqrt(FourierMap >& m) const; }; template @@ -140,7 +143,7 @@ namespace CosmoTool double k2 = 0; for (int q = 0; q < ik.size(); q++) { - int dk = ik; + int dk = ik[q]; if (dk > dims[q]/2) dk = dk - dims[q]; @@ -349,18 +352,21 @@ namespace CosmoTool EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const { typedef EuclidianFourierMapComplex MapT; + typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; typedef typename MapT::DimArray DimArray; - MapT& m_c = dynamic_cast(like_map); - MapT *rand_map = m_c.mimick(); - std::complex *d = rand_map->data(); + const MapT& m_c = dynamic_cast(like_map); + ptr_map ret_map = ptr_map(m_c.mimick()); + MapT& rand_map = dynamic_cast(*ret_map.get()); + + std::complex *d = rand_map.data(); long idx; - const DimArray& dims = rand_map->getDims(); + const DimArray& dims = rand_map.getDims(); long plane_size; for (long p = 1; p < m_c.size(); p++) { - double A_k = std::sqrt(0.5*f(rand_map->get_K(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)); } @@ -372,7 +378,7 @@ namespace CosmoTool if (dims.size() == 1) { d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); - return boost::shared_ptr::FourierMapType>(rand_map); + return ret_map; } plane_size = 1; @@ -391,6 +397,7 @@ namespace CosmoTool } long q = dims[0]; d[q] = std::complex(d[q].real() + d[q].imag()); + return ret_map; } template @@ -400,8 +407,51 @@ namespace CosmoTool std::complex *d = m.data(); for (long p = 0; p < m_c.size(); p++) - d[p] *= f(m.get_K(p)); + d[p] *= f(m_c.get_K(p)); } + + template + void EuclidianSpectrum_1D::mul_sqrt(FourierMap >& 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(FourierMap >& 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(FourierMap >& 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/healpix.hpp b/src/fourier/healpix.hpp index 8aac379..06da0fc 100644 --- a/src/fourier/healpix.hpp +++ b/src/fourier/healpix.hpp @@ -32,8 +32,13 @@ namespace CosmoTool const T *data() const { return &cls[0]; } long size() const { return cls.size(); } - ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; - }; + ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; + + 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 HealpixFourierMap: public FourierMap @@ -248,13 +253,13 @@ namespace CosmoTool long lmaxGen = std::min(cls.size()-1, alms.Lmax()); std::complex *new_data = new_alms->data(); - for (long l = 0; l < lmaxGen; l++) + for (long l = 0; l <= lmaxGen; l++) { double Al = std::sqrt(cls[l]); new_data[alms.index(l,0)] = gsl_ran_gaussian(rng, Al); Al *= M_SQRT1_2; - for (long m = 1; m < alms.Mmax(); m++) + for (long m = 1; m <= std::min(l,alms.Mmax()); m++) { std::complex& c = new_data[alms.index(l,m)]; c.real() = gsl_ran_gaussian(rng, Al); @@ -263,6 +268,75 @@ namespace CosmoTool } return r; } + + template + void HealpixSpectrum::mul(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (long l = 0; l <= alms.Lmax(); l++) + { + double Al = cls[l]; + + for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (long l = 0; l <= alms.Lmax(); l++) + { + double Al = std::sqrt(cls[l]); + + for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_inv(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (long l = 0; l <= alms.Lmax(); l++) + { + double Al = (cls[l] <= 0) ? 0 : (1/cls[l]); + + for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_inv_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (long l = 0; l <= alms.Lmax(); l++) + { + double Al = (cls[l] <= 0) ? 0 : std::sqrt(1/cls[l]); + + for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + }; #endif