From 414edbd28a777c269e7a4cc9a77e6d61de71bde1 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 10 Nov 2012 09:02:08 -0500 Subject: [PATCH] Changed dependency to libsharp (MPI support, better performance). Added new power spectrum abstraction --- CMakeLists.txt | 3 + sample/CMakeLists.txt | 11 +- src/fourier/base_types.hpp | 36 +++++- src/fourier/euclidian.hpp | 233 ++++++++++++++++++++++++++++++++++--- src/fourier/healpix.hpp | 153 +++++++++++++++++++----- 5 files changed, 390 insertions(+), 46 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6c4dd9b..e04a42b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,9 @@ find_library(NETCDFCPP_LIBRARY NAMES netcdf_c++ netcdf_c++4) find_library(GSL_LIBRARY gsl) find_library(GSLCBLAS_LIBRARY gslcblas) +find_library(SHARP_LIBRARY sharp) +find_path(SHARP_INCLUDE_PATH sharp.h) + set(HDF5_FIND_COMPONENTS HL CXX) if(HDF5_ROOTDIR) SET(ENV{HDF5_ROOT} ${HDF5_ROOTDIR}) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 3a18be0..a437ff9 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -1,6 +1,10 @@ SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS}) include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH}) +IF(SHARP_INCLUDE_PATH) + include_directories(BEFORE ${SHARP_INCLUDE_PATH}) +ENDIF(SHARP_INCLUDE_PATH) + add_executable(testBQueue testBQueue.cpp) target_link_libraries(testBQueue ${tolink}) @@ -41,6 +45,11 @@ add_executable(testBSP testBSP.cpp) target_link_libraries(testBSP ${tolink}) if (FFTW3_FOUND AND EIGEN3_FOUND) - add_executable(test_fft_calls test_fft_calls) + add_executable(test_fft_calls test_fft_calls.cpp) target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES}) endif (FFTW3_FOUND AND EIGEN3_FOUND) + +if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) + add_executable(test_sharp_calls test_sharp_calls.cpp) + target_link_libraries(test_sharp_calls ${tolink} ${SHARP_LIBRARY}) +endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) \ No newline at end of file diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index 9a02a7e..9f54320 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -1,12 +1,37 @@ #ifndef __BASE_FOURIER_TYPES_HPP #define __BASE_FOURIER_TYPES_HPP +#include +#include #include #include #include +#include namespace CosmoTool { + class IncompatibleMap: virtual std::exception {}; + + template class FourierMap; + + template + class SpectrumFunction + { + protected: + SpectrumFunction() {} + public: + typedef Eigen::Array VecType; + typedef Eigen::Map MapType; + typedef Eigen::Map ConstMapType; + typedef FourierMap > FourierMapType; + + virtual boost::shared_ptr + newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; + + virtual void mul(FourierMap >& m) const = 0; + }; + + template class FourierMap { @@ -14,7 +39,7 @@ namespace CosmoTool FourierMap() {} public: - typedef Eigen::Matrix VecType; + typedef Eigen::Array VecType; typedef Eigen::Map MapType; typedef Eigen::Map ConstMapType; @@ -34,6 +59,12 @@ namespace CosmoTool { return ConstMapType(data(), size()); } + + void sqrt() + { + MapType m = eigen(); + m = m.sqrt(); + } void scale(const T& factor) { @@ -71,6 +102,9 @@ namespace CosmoTool return m; } + virtual T dot_product(const FourierMap& second) const + throw(std::bad_cast) = 0; + virtual FourierMap *mimick() const = 0; }; diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index 0f40e1a..a29c623 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -1,23 +1,49 @@ #ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP #define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP +#include #include #include +#include #include "base_types.hpp" #include "fft/fftw_calls.hpp" +#include "../algo.hpp" namespace CosmoTool { + template + class EuclidianSpectrum_1D: public SpectrumFunction + { + public: + typedef boost::function1 Function; + protected: + Function f; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef boost::shared_ptr ptr_map; + + EuclidianSpectrum_1D(Function P) + : f(P) + { + } + + ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; + + void mul(FourierMap >& m) const; + }; template - class EuclidianFourierMap: public FourierMap + class EuclidianFourierMapBase: public FourierMap { + public: + typedef std::vector DimArray; private: boost::shared_ptr m_data; - std::vector m_dims; + DimArray m_dims; long m_size; public: - EuclidianFourierMap(boost::shared_ptr indata, std::vector indims) + + EuclidianFourierMapBase(boost::shared_ptr indata, const DimArray& indims) { m_data = indata; m_dims = indims; @@ -26,10 +52,12 @@ namespace CosmoTool m_size *= m_dims[i]; } - virtual ~EuclidianFourierMap() + 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; } @@ -41,33 +69,150 @@ namespace CosmoTool 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 EuclidianFourierMap( - boost::shared_ptr((T *)fftw_malloc(sizeof(T)*m_size), + return new EuclidianFourierMapReal( + boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), std::ptr_fun(fftw_free)), - m_dims); + 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; + long plane_size; + public: + typedef typename EuclidianFourierMapBase >::DimArray DimArray; + + EuclidianFourierMapComplex(ptr_t indata, + const DimArray& indims, + const std::vector& dk) + : EuclidianFourierMapBase >(indata, indims), delta_k(dk) + { + assert(dk.size() == indims.size()); + plane_size = 1; + for (int q = 1; q < indims.size(); q++) + plane_size *= indims[q]; + } + + virtual FourierMap > *mimick() const + { + return + new EuclidianFourierMapComplex( + ptr_t((std::complex *) + fftw_malloc(sizeof(std::complex)*this->size()), + std::ptr_fun(fftw_free)), + 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; + for (int q = 0; q < ik.size(); q++) + { + int dk = ik; + + 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); + } + + 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]; + std::complex result = 0; + + for (long q0 = 1; q0 < N0-1; q0++) + { + for (long p = 0; p < plane_size; p++) + { + result += 2*(conj(d1[q0+N0*p]) * d2[q0+N0*p]).real(); + } + } + 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; - EuclidianFourierMap *realMap; - EuclidianFourierMap > *fourierMap; + EuclidianFourierMapReal *realMap; + EuclidianFourierMapComplex *fourierMap; typename calls::plan_type m_analysis, m_synthesis; double volume; long N; - std::vector m_dims; + DimArray m_dims; std::vector m_L; public: - EuclidianFourierTransform(const std::vector& dims, const std::vector& L) + EuclidianFourierTransform(const DimArray& dims, const std::vector& L) { assert(L.size() == dims.size()); - + std::vector dk(L.size()); + m_dims = dims; m_L = L; @@ -77,16 +222,17 @@ namespace CosmoTool { N *= dims[i]; volume *= L[i]; + dk[i] = 2*M_PI/L[i]; } - realMap = new EuclidianFourierMap( + realMap = new EuclidianFourierMapReal( boost::shared_ptr(calls::alloc_real(N), std::ptr_fun(calls::free)), dims); - fourierMap = new EuclidianFourierMap >( + fourierMap = new EuclidianFourierMapComplex( boost::shared_ptr >((std::complex*)calls::alloc_complex(N), std::ptr_fun(calls::free)), - dims); + dims, dk); m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0], realMap->data(), (typename calls::complex_type *)fourierMap->data(), FFTW_MEASURE); @@ -194,7 +340,64 @@ namespace CosmoTool }; + template + typename EuclidianSpectrum_1D::ptr_map + EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const + { + typedef EuclidianFourierMapComplex MapT; + typedef typename MapT::DimArray DimArray; + MapT& m_c = dynamic_cast(like_map); + MapT *rand_map = m_c.mimick(); + std::complex *d = rand_map->data(); + long idx; + 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))); + d[p] = std::complex(gsl_ran_gaussian(rng, A_k), + gsl_ran_gaussian(rng, A_k)); + } + // Generate the mean value + d[0] = std::complex(std::sqrt(f(0)), 0); + // Correct the Nyquist plane + idx = dims[0]-1; // Stick to the last element of the first dimension + // 1D is special case + if (dims.size() == 1) + { + d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); + return boost::shared_ptr::FourierMapType>(rand_map); + } + + 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]); + } + long q = dims[0]; + d[q] = std::complex(d[q].real() + d[q].imag()); + } + + template + void EuclidianSpectrum_1D::mul(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] *= f(m.get_K(p)); + } }; #endif diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp index 1759a78..11ffddd 100644 --- a/src/fourier/healpix.hpp +++ b/src/fourier/healpix.hpp @@ -1,24 +1,45 @@ #ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP #define __COSMOTOOL_FOURIER_HEALPIX_HPP +#include +#include +#include #include #include -#include -#include -#include +#include +#include namespace CosmoTool { - template - class HealpixFourierMap: public FourierMap, public Healpix_Base + class HealpixSpectrum: public FourierSpectrum + { + protected: + std::vector cls; + public: + typedef typename FourierSpectrum::FourierMapType FourierMapType; + + HealpixSpectrum(long Lmax) + : cls(Lmax+1) {} + + T *data() { return &cls[0]; } + const T *data() const { return &cls[0]; } + long size() const { return cls.size(); } + + boost::shared_ptr + newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; + }; + + template + class HealpixFourierMap: public FourierMap { private: T *m_data; + long Npix, Nside; Eigen::aligned_allocator alloc; public: HealpixFourierMap(long nSide) - : Healpix_Base(RING, nSide, SET_NSIDE) + : Npix(12*nSide*nSide), Nside(nSide) { m_data = alloc.allocate(Npix); } @@ -30,27 +51,61 @@ namespace CosmoTool virtual const T* data() const { return m_data; } virtual T *data() { return m_data; } - virtual long size() const { return Npix(); } + virtual long size() const { return Npix; } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + const HealpixFourierMap& mfm = dynamic_cast&>(other); + if (Npix() != mfm.Npix()) + 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(Nside()); + return new HealpixFourierMap(Nside); } }; template - class HealpixFourierALM: public FourierMap >, public Alm_Base + class HealpixFourierALM: public FourierMap > { private: std::complex *alms; long size; + long Lmax_, Mmax_, TVal_; Eigen::aligned_allocator > alloc; public: - HealpixFourierALM(long Lmax, long Mmax) - : Alm_Base(Lmax, Mmax) + typedef unsigned long LType; + + LType Lmax() const { return Lmax_; } + LType Mmax() const { return Mmax_; } + + LType Num_Alms() const { - size = Num_Alms(Lmax, Mmax); + 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) + { + size = Num_Alms(); alms = alloc.allocate(size); } @@ -63,25 +118,41 @@ namespace CosmoTool virtual T * data() { return alms;} virtual long size() const { return size; } - virtual FourierMap *mimick() const + virtual FourierMap > *mimick() const { - return new HealpixFourierALM(Lmax(), Mmax()); + return new HealpixFourierALM(Lmax_, Mmax_); + } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const HealpixFourierALM& mfm = dynamic_cast&>(other); + std::complex S; + + if (size != mfm.size) + throw std::bad_cast(); + + MapType m1(m_data, Npix); + MapType m2(mfm.m_data, mfm.Npix); + + S = (m1.block(0,0,1,Lmax_+1).adjoint() * m2(0,0,1,Lmax_+1)).sum(); + S += 2*(m1.block(0,1+Lmax_,1,size-1-Lmax_).adjoint() * m2(0,1+Lmax_,1,size-1-Lmax_)).sum(); + return S; } }; template - class HealpixFourierTransform: public FourierTransform + class HealpixFourierTransform: public FourierTransform, public sharp_base { private: HealpixFourierMap realMap; HealpixFourierALM fourierMap; - psht_joblist jobs; public: HealpixFourierTransform(long nSide, long Lmax, long Mmax) : realMap(nSide), fourierMap(Lmax, Mmax) { - jobs.set_Healpix_geometry(nSide); - jobs.set_triangular_alm_info(Lmax, Mmax); + set_Healpix_geometry(nSide); + set_triangular_alm_info(Lmax, Mmax); } virtual ~HealpixFourierTransform() {} @@ -96,25 +167,23 @@ namespace CosmoTool virtual FourierTransform *mimick() const { - return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); + return new HealpixFourierTransform(realMap.Nside, fourierMap.Lmax, fourierMap.Mmax); } virtual void analysis() { - jobs.add_map2alm(realMap.data(), - reinterpret_cast *>(fourierMap.data()), - false); - jobs.execute(); - jobs.clear_jobs(); + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + cxxjobhelper__::val,0,0,0); } virtual void synthesis() { - jobs.add_alm2map(reinterpret_cast *>(fourierMap.data()), - realMap.data(), - false); - jobs.execute(); - jobs.clear_jobs(); + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + cxxjobhelper__::val,0,0,0); } virtual void analysis_conjugate() @@ -130,6 +199,32 @@ namespace CosmoTool } }; + + template + boost::shared_ptr::FourierMapType> + HealpixSpectrum::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const + { + const HealpixFourierALM& alms = dynamic_cast&>(like_map); + HealpixFourierALM *new_alms; + boost::shared_ptr r(new_alms = new HealpixFourierALM(alms.Lmax(), alms.Mmax())); + long lmaxGen = std::min(cls.size()-1, alms.Lmax()); + std::complex *new_data = new_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 < 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); + } + } + return r; + } }; #endif