From bb20fcee059fc5f49c5c775c97ba25f1c1a34249 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 18 Nov 2012 16:36:18 -0500 Subject: [PATCH] Check healpix utility calls --- sample/test_healpix_calls.cpp | 5 ++++ src/fourier/details/healpix_spectrum.hpp | 30 +++++++++++++----------- src/fourier/details/healpix_utility.hpp | 19 ++++++++------- 3 files changed, 32 insertions(+), 22 deletions(-) diff --git a/sample/test_healpix_calls.cpp b/sample/test_healpix_calls.cpp index d1f1bf3..7cc4272 100644 --- a/sample/test_healpix_calls.cpp +++ b/sample/test_healpix_calls.cpp @@ -7,14 +7,19 @@ using namespace std; int main() { HealpixFourierTransform dft(128,2*128,2*128, 40); + HealpixUtilityFunction utils; long Npix = dft.realSpace().size(); dft.realSpace().eigen().setRandom(); dft.analysis(); cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*Npix/(4*M_PI) << endl; + cout << "Spectrum analysis" << endl; + HealpixUtilityFunction::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace()); + dft.synthesis(); cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; + return 0; } diff --git a/src/fourier/details/healpix_spectrum.hpp b/src/fourier/details/healpix_spectrum.hpp index e797c91..cfeebc9 100644 --- a/src/fourier/details/healpix_spectrum.hpp +++ b/src/fourier/details/healpix_spectrum.hpp @@ -13,8 +13,9 @@ namespace CosmoTool typedef typename SpectrumFunction::FourierMapType FourierMapType; typedef boost::shared_ptr ptr_map; typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef unsigned long LType; - HealpixSpectrum(long Lmax) + HealpixSpectrum(LType Lmax) : cls(Lmax+1), m_dof(new int[Lmax+1]) { for (int l = 0; l <= Lmax; l++) @@ -25,14 +26,15 @@ namespace CosmoTool const T *data() const { return &cls[0]; } const int *dof() const { return m_dof; } - void set_dof(long l, int dof) { m_dof[l] = dof; } + void set_dof(LType l, int dof) { m_dof[l] = dof; } + LType Lmax() const { return cls.size()-1; } long size() const { return cls.size(); } void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; SpectrumFunctionPtr copy() const { - HealpixSpectrum *s = new HealpixSpectrum(cls.size()-1); + HealpixSpectrum *s = new HealpixSpectrum(Lmax()); s->cls = cls; return SpectrumFunctionPtr(s); } @@ -54,13 +56,13 @@ namespace CosmoTool long lmaxGen = std::min(cls.size()-1, alms.Lmax()); std::complex *new_data = alms.data(); - for (long l = 0; l <= lmaxGen; l++) + for (LType 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 <= std::min(l,alms.Mmax()); m++) + for (LType 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); @@ -75,11 +77,11 @@ namespace CosmoTool HealpixFourierALM& alms = dynamic_cast&>(like_map); std::complex *data = alms.data(); - for (long l = 0; l <= alms.Lmax(); l++) + for (LType l = 0; l <= alms.Lmax(); l++) { double Al = cls[l]; - for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) { data[alms.index(l,m)] *= Al; } @@ -89,14 +91,14 @@ namespace CosmoTool template void HealpixSpectrum::mul_sqrt(FourierMapType& like_map) const { - HealpixFourierALM& alms = dynamic_cast&>(like_map); + HealpixFourierALM& alms = dynamic_cast&>(like_map); std::complex *data = alms.data(); - for (long l = 0; l <= alms.Lmax(); l++) + for (LType l = 0; l <= alms.Lmax(); l++) { double Al = std::sqrt(cls[l]); - for (long m = 0; m <= std::min(l,alms.Mmax()); m++) + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) { data[alms.index(l,m)] *= Al; } @@ -109,11 +111,11 @@ namespace CosmoTool HealpixFourierALM& alms = dynamic_cast&>(like_map); std::complex *data = alms.data(); - for (long l = 0; l <= alms.Lmax(); l++) + for (LType 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++) + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) { data[alms.index(l,m)] *= Al; } @@ -126,11 +128,11 @@ namespace CosmoTool HealpixFourierALM& alms = dynamic_cast&>(like_map); std::complex *data = alms.data(); - for (long l = 0; l <= alms.Lmax(); l++) + for (LType 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++) + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) { data[alms.index(l,m)] *= Al; } diff --git a/src/fourier/details/healpix_utility.hpp b/src/fourier/details/healpix_utility.hpp index b527366..004208b 100644 --- a/src/fourier/details/healpix_utility.hpp +++ b/src/fourier/details/healpix_utility.hpp @@ -11,28 +11,31 @@ namespace CosmoTool typedef typename MapUtilityFunction::Spectrum Spectrum; typedef typename MapUtilityFunction::Spectrum_ptr Spectrum_ptr; typedef typename MapUtilityFunction::FMap FMap; - + typedef typename HealpixSpectrum::LType LType; + Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const { const HealpixFourierALM& alms = dynamic_cast&>(m); - long Lmax = alms.Lmax(), Mmax = alms.Mmax(); + LType Lmax = alms.Lmax(), Mmax = alms.Mmax(); HealpixSpectrum *spectrum = new HealpixSpectrum(Lmax); T *data_cls = spectrum->data(); - std::complex *data_alms = alms.data(); + const std::complex *data_alms = alms.data(); // make an estimate of the cls std::fill(data_cls, data_cls+Lmax+1, 0); - for (long m = 0, q = 0; m <= Mmax; ++m ) + LType q = 0; + for (LType m = 0; m <= Mmax; ++m ) { - for (long l = 0; l <= Lmax; ++l) + for (LType l = m; l <= Lmax; ++l) { // Triangular storage - data_cls[l] += std::norm(alms[l+q]); + data_cls[l] += std::norm(data_alms[l+q]); } - q += Lmax+1; + q += Lmax+1-m; } + assert(q==alms.size()); - for (long l = 0; l <= Lmax; ++l) + for (LType l = 0; l <= Lmax; ++l) { int dof = 1 + std::min(l, Mmax); spectrum->set_dof(l, dof);