From b7156fcc76f0fb54d27e0c6e823877635742b4b8 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 18 Nov 2012 15:47:37 -0500 Subject: [PATCH] Splitted healpix.hpp into several files for easier maintainability --- src/fourier/details/healpix_alms.hpp | 76 +++++ src/fourier/details/healpix_map.hpp | 54 ++++ src/fourier/details/healpix_spectrum.hpp | 133 +++++++++ src/fourier/details/healpix_transform.hpp | 97 +++++++ src/fourier/healpix.hpp | 333 +--------------------- 5 files changed, 364 insertions(+), 329 deletions(-) create mode 100644 src/fourier/details/healpix_alms.hpp create mode 100644 src/fourier/details/healpix_map.hpp create mode 100644 src/fourier/details/healpix_spectrum.hpp create mode 100644 src/fourier/details/healpix_transform.hpp diff --git a/src/fourier/details/healpix_alms.hpp b/src/fourier/details/healpix_alms.hpp new file mode 100644 index 0000000..849781f --- /dev/null +++ b/src/fourier/details/healpix_alms.hpp @@ -0,0 +1,76 @@ +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP + +namespace CosmoTool +{ + template + class HealpixFourierALM: public FourierMap > + { + private: + std::complex *alms; + long m_size; + long Lmax_, Mmax_, TVal_; + Eigen::aligned_allocator > alloc; + public: + typedef unsigned long LType; + + LType Lmax() const { return Lmax_; } + LType Mmax() const { return Mmax_; } + + LType Num_Alms() const + { + return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_); + } + + LType index_l0(LType m) const + { + return ((m*(TVal_-m))/2); + } + + LType index(LType l, LType m) const + { + return index_l0(m) + l; + } + + HealpixFourierALM(LType lmax, LType mmax) + : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1) + { + m_size = Num_Alms(); + alms = alloc.allocate(m_size); + } + + virtual ~HealpixFourierALM() + { + alloc.deallocate(alms, m_size); + } + + virtual const std::complex* data() const { return alms; } + virtual std::complex * data() { return alms;} + virtual long size() const { return m_size; } + + virtual FourierMap > *mimick() const + { + return new HealpixFourierALM(Lmax_, Mmax_); + } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const HealpixFourierALM& mfm = dynamic_cast&>(other); + typedef typename FourierMap >::MapType MapType; + std::complex S; + + if (m_size != mfm.m_size) + throw std::bad_cast(); + + MapType m1(alms, m_size); + MapType m2(mfm.alms, mfm.m_size); + + S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum(); + S += std::complex(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum(); + return S; + } + }; +}; + +#endif diff --git a/src/fourier/details/healpix_map.hpp b/src/fourier/details/healpix_map.hpp new file mode 100644 index 0000000..6b96974 --- /dev/null +++ b/src/fourier/details/healpix_map.hpp @@ -0,0 +1,54 @@ +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP + +namespace CosmoTool +{ + + template + class HealpixFourierMap: public FourierMap + { + private: + T *m_data; + long Npix, m_Nside; + Eigen::aligned_allocator alloc; + public: + HealpixFourierMap(long nSide) + : Npix(12*nSide*nSide), m_Nside(nSide) + { + m_data = alloc.allocate(Npix); + } + + virtual ~HealpixFourierMap() + { + alloc.deallocate(m_data, Npix); + } + + long Nside() const { return m_Nside; } + virtual const T* data() const { return m_data; } + virtual T *data() { return m_data; } + virtual long size() const { return Npix; } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + typedef typename FourierMap::MapType MapType; + + const HealpixFourierMap& mfm = dynamic_cast&>(other); + if (Npix != mfm.size()) + throw std::bad_cast(); + + MapType m1(m_data, Npix); + MapType m2(mfm.m_data, mfm.Npix); + + return (m1*m2).sum(); + } + + virtual FourierMap *mimick() const + { + return new HealpixFourierMap(m_Nside); + } + }; + +}; + +#endif diff --git a/src/fourier/details/healpix_spectrum.hpp b/src/fourier/details/healpix_spectrum.hpp new file mode 100644 index 0000000..2031acf --- /dev/null +++ b/src/fourier/details/healpix_spectrum.hpp @@ -0,0 +1,133 @@ +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP + +namespace CosmoTool +{ + template + class HealpixSpectrum: public SpectrumFunction + { + protected: + std::vector cls; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef boost::shared_ptr ptr_map; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + + HealpixSpectrum(long Lmax) + : cls(Lmax+1) {} + + T *data() { return &cls[0]; } + const T *data() const { return &cls[0]; } + 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); + s->cls = cls; + return SpectrumFunctionPtr(s); + } + + void sqrt() { + std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun(std::sqrt)); + } + + 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 HealpixSpectrum::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(out_map); + long lmaxGen = std::min(cls.size()-1, alms.Lmax()); + std::complex *new_data = alms.data(); + + 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 <= std::min(l,alms.Mmax()); m++) + { + std::complex& c = new_data[alms.index(l,m)]; + c.real() = gsl_ran_gaussian(rng, Al); + c.imag() = gsl_ran_gaussian(rng, Al); + } + } + } + + 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 diff --git a/src/fourier/details/healpix_transform.hpp b/src/fourier/details/healpix_transform.hpp new file mode 100644 index 0000000..bc542e7 --- /dev/null +++ b/src/fourier/details/healpix_transform.hpp @@ -0,0 +1,97 @@ +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP + +namespace CosmoTool +{ + + template struct HealpixJobHelper__ {}; + + template<> struct HealpixJobHelper__ + { enum {val=1}; }; + + template<> struct HealpixJobHelper__ + { enum {val=0}; }; + + + template + class HealpixFourierTransform: public FourierTransform + { + private: + sharp_alm_info *ainfo; + sharp_geom_info *ginfo; + HealpixFourierMap realMap; + HealpixFourierALM fourierMap; + int m_iterate; + public: + HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0) + : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate) + { + sharp_make_healpix_geom_info (nSide, 1, &ginfo); + sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo); + } + + virtual ~HealpixFourierTransform() + { + sharp_destroy_geom_info(ginfo); + sharp_destroy_alm_info(ainfo); + } + + virtual const FourierMap >& fourierSpace() const { return fourierMap; } + + virtual FourierMap >& fourierSpace() { return fourierMap; } + + virtual const FourierMap& realSpace() const { return realMap; } + + virtual FourierMap& realSpace() { return realMap; } + + virtual FourierTransform *mimick() const + { + return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); + } + + virtual void analysis() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + for (int i = 0; i < m_iterate; i++) + { + HealpixFourierMap tmp_map(realMap.Nside()); + void *tmp_ptr=reinterpret_cast(tmp_map.data()); + typename HealpixFourierMap::MapType m0 = tmp_map.eigen(); + typename HealpixFourierMap::MapType m1 = realMap.eigen(); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + m0 = m1 - m0; + sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + } + } + + virtual void synthesis() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + } + + virtual void analysis_conjugate() + { + synthesis(); + realMap.scale(4*M_PI/realMap.size()); + } + + virtual void synthesis_conjugate() + { + analysis(); + fourierMap.scale(realMap.size()/(4*M_PI)); + } + + }; + +}; + +#endif diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp index 460935d..46f5356 100644 --- a/src/fourier/healpix.hpp +++ b/src/fourier/healpix.hpp @@ -15,335 +15,10 @@ #include #include -namespace CosmoTool -{ - template - class HealpixSpectrum: public SpectrumFunction - { - protected: - std::vector cls; - public: - typedef typename SpectrumFunction::FourierMapType FourierMapType; - typedef boost::shared_ptr ptr_map; - typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; +#include "details/healpix_map.hpp" +#include "details/healpix_alms.hpp" +#include "details/healpix_transform.hpp" +#include "details/healpix_spectrum.hpp" - HealpixSpectrum(long Lmax) - : cls(Lmax+1) {} - - T *data() { return &cls[0]; } - const T *data() const { return &cls[0]; } - 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); - s->cls = cls; - return SpectrumFunctionPtr(s); - } - - void sqrt() { - std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun(std::sqrt)); - } - - 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 - { - private: - T *m_data; - long Npix, m_Nside; - Eigen::aligned_allocator alloc; - public: - HealpixFourierMap(long nSide) - : Npix(12*nSide*nSide), m_Nside(nSide) - { - m_data = alloc.allocate(Npix); - } - - virtual ~HealpixFourierMap() - { - alloc.deallocate(m_data, Npix); - } - - long Nside() const { return m_Nside; } - virtual const T* data() const { return m_data; } - virtual T *data() { return m_data; } - virtual long size() const { return Npix; } - - virtual T dot_product(const FourierMap& other) const - throw(std::bad_cast) - { - typedef typename FourierMap::MapType MapType; - - const HealpixFourierMap& mfm = dynamic_cast&>(other); - if (Npix != mfm.size()) - throw std::bad_cast(); - - MapType m1(m_data, Npix); - MapType m2(mfm.m_data, mfm.Npix); - - return (m1*m2).sum(); - } - - virtual FourierMap *mimick() const - { - return new HealpixFourierMap(m_Nside); - } - }; - - - template - class HealpixFourierALM: public FourierMap > - { - private: - std::complex *alms; - long m_size; - long Lmax_, Mmax_, TVal_; - Eigen::aligned_allocator > alloc; - public: - typedef unsigned long LType; - - LType Lmax() const { return Lmax_; } - LType Mmax() const { return Mmax_; } - - LType Num_Alms() const - { - return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_); - } - - LType index_l0(LType m) const - { - return ((m*(TVal_-m))/2); - } - - LType index(LType l, LType m) const - { - return index_l0(m) + l; - } - - HealpixFourierALM(LType lmax, LType mmax) - : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1) - { - m_size = Num_Alms(); - alms = alloc.allocate(m_size); - } - - virtual ~HealpixFourierALM() - { - alloc.deallocate(alms, m_size); - } - - virtual const std::complex* data() const { return alms; } - virtual std::complex * data() { return alms;} - virtual long size() const { return m_size; } - - virtual FourierMap > *mimick() const - { - return new HealpixFourierALM(Lmax_, Mmax_); - } - - virtual std::complex dot_product(const FourierMap >& other) const - throw(std::bad_cast) - { - const HealpixFourierALM& mfm = dynamic_cast&>(other); - typedef typename FourierMap >::MapType MapType; - std::complex S; - - if (m_size != mfm.m_size) - throw std::bad_cast(); - - MapType m1(alms, m_size); - MapType m2(mfm.alms, mfm.m_size); - - S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum(); - S += std::complex(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum(); - return S; - } - }; - - template struct HealpixJobHelper__ {}; - - template<> struct HealpixJobHelper__ - { enum {val=1}; }; - - template<> struct HealpixJobHelper__ - { enum {val=0}; }; - - - template - class HealpixFourierTransform: public FourierTransform - { - private: - sharp_alm_info *ainfo; - sharp_geom_info *ginfo; - HealpixFourierMap realMap; - HealpixFourierALM fourierMap; - int m_iterate; - public: - HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0) - : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate) - { - sharp_make_healpix_geom_info (nSide, 1, &ginfo); - sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo); - } - - virtual ~HealpixFourierTransform() - { - sharp_destroy_geom_info(ginfo); - sharp_destroy_alm_info(ainfo); - } - - virtual const FourierMap >& fourierSpace() const { return fourierMap; } - - virtual FourierMap >& fourierSpace() { return fourierMap; } - - virtual const FourierMap& realSpace() const { return realMap; } - - virtual FourierMap& realSpace() { return realMap; } - - virtual FourierTransform *mimick() const - { - return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); - } - - virtual void analysis() - { - void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); - - sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - for (int i = 0; i < m_iterate; i++) - { - HealpixFourierMap tmp_map(realMap.Nside()); - void *tmp_ptr=reinterpret_cast(tmp_map.data()); - typename HealpixFourierMap::MapType m0 = tmp_map.eigen(); - typename HealpixFourierMap::MapType m1 = realMap.eigen(); - - sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - m0 = m1 - m0; - sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - } - } - - virtual void synthesis() - { - void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); - - sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - } - - virtual void analysis_conjugate() - { - synthesis(); - realMap.scale(4*M_PI/realMap.size()); - } - - virtual void synthesis_conjugate() - { - analysis(); - fourierMap.scale(realMap.size()/(4*M_PI)); - } - - }; - - template - void HealpixSpectrum::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const - { - HealpixFourierALM& alms = dynamic_cast&>(out_map); - long lmaxGen = std::min(cls.size()-1, alms.Lmax()); - std::complex *new_data = alms.data(); - - 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 <= std::min(l,alms.Mmax()); m++) - { - std::complex& c = new_data[alms.index(l,m)]; - c.real() = gsl_ran_gaussian(rng, Al); - c.imag() = gsl_ran_gaussian(rng, Al); - } - } - } - - 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