Added a test for the healpix/sharp tranform. Basic code compiles and work.

CMakeLists.txt

@@ -5,11 +5,13 @@ project(CosmoToolbox)
 
 include(GetGitRevisionDescription)
 include(ExternalProject)
+include(FindOpenMP)
 
 get_git_head_revision(HEAD GIT_VER)
 
 option(BUILD_SHARED_LIBS "Build shared libraries." OFF)
 option(BUILD_STATIC_LIBS "Build static libraries." ON)
+option(ENABLE_OPENMP "Enable OpenMP support." OFF)
 
 find_path(NETCDF_INCLUDE_PATH NAMES netcdf.h)
 find_path(NETCDFCPP_INCLUDE_PATH NAMES netcdfcpp.h netcdf)
@@ -39,6 +41,9 @@ ExternalProject_Add(sharp
 )
 
 SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a)
+SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a)
+SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a)
+SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY})
 SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include)
 
 set(HDF5_FIND_COMPONENTS HL CXX)

sample/CMakeLists.txt

@@ -49,7 +49,9 @@ if (FFTW3_FOUND AND EIGEN3_FOUND)
   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)
+if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
-#  add_executable(test_sharp_calls test_sharp_calls.cpp)
+  include_directories(${SHARP_INCLUDE_PATH})
-#  target_link_libraries(test_sharp_calls ${tolink} ${SHARP_LIBRARY})
+  add_executable(test_healpix_calls test_healpix_calls.cpp)
-#endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
+  target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES})
+  set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
+endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)

sample/test_healpix_calls.cpp

@@ -0,0 +1,20 @@
+#include <iostream>
+#include "fourier/healpix.hpp"
+
+using namespace CosmoTool;
+using namespace std;
+
+int main()
+{
+  HealpixFourierTransform<double> dft(128,3*128,3*128, 40);
+  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;
+  dft.synthesis();
+  cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+  return 0;
+
+}

src/fourier/healpix.hpp

@@ -4,20 +4,26 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 #include <cmath>
+#include <vector>
 #include <boost/bind.hpp>
 #include <boost/shared_ptr.hpp>
 #include <exception>
-#include <sharp_cxx.h>
+#include "base_types.hpp"
+#include <sharp_lowlevel.h>
+#include <sharp_geomhelpers.h>
+#include <sharp_almhelpers.h>
 
 namespace CosmoTool
 {
   template<typename T>
-  class HealpixSpectrum: public FourierSpectrum<T>
+  class HealpixSpectrum: public SpectrumFunction<T>
   {
   protected:
     std::vector<T> cls;
   public:
-    typedef typename FourierSpectrum<T>::FourierMapType FourierMapType;
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
 
     HealpixSpectrum(long Lmax)
      : cls(Lmax+1) {}
@@ -26,8 +32,7 @@ namespace CosmoTool
     const T *data() const { return &cls[0]; }
     long size() const { return cls.size(); }
 
-    boost::shared_ptr<FourierMapType> 
+    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
-       newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
   };    
 
   template<typename T>
@@ -35,20 +40,21 @@ namespace CosmoTool
   {
   private:
     T *m_data;
-    long Npix, Nside;
+    long Npix, m_Nside;
     Eigen::aligned_allocator<T> alloc;
   public:
     HealpixFourierMap(long nSide)
-      : Npix(12*nSide*nSide), Nside(nSide)
+      : Npix(12*nSide*nSide), m_Nside(nSide)
     {
       m_data = alloc.allocate(Npix);
     }
 
     virtual ~HealpixFourierMap()
     {
-      alloc.deallocate(m_data);
+      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; }
@@ -56,8 +62,10 @@ namespace CosmoTool
     virtual T dot_product(const FourierMap<T>& other) const
       throw(std::bad_cast)
     {
+      typedef typename FourierMap<T>::MapType MapType;
+
       const HealpixFourierMap<T>& mfm = dynamic_cast<const HealpixFourierMap<T>&>(other);
-      if (Npix() != mfm.Npix())
+      if (Npix != mfm.size())
         throw std::bad_cast();
 
       MapType m1(m_data, Npix);
@@ -68,7 +76,7 @@ namespace CosmoTool
     
     virtual FourierMap<T> *mimick() const
     {
-      return new HealpixFourierMap<T>(Nside);
+      return new HealpixFourierMap<T>(m_Nside);
     }
   };
 
@@ -78,7 +86,7 @@ namespace CosmoTool
   {
   private:
     std::complex<T> *alms;
-    long size;
+    long m_size;
     long Lmax_, Mmax_, TVal_;
     Eigen::aligned_allocator<std::complex<T> > alloc;
   public:
@@ -105,18 +113,18 @@ namespace CosmoTool
     HealpixFourierALM(LType lmax, LType mmax)
       : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1)
     {
-      size = Num_Alms();
+      m_size = Num_Alms();
-      alms = alloc.allocate(size);
+      alms = alloc.allocate(m_size);
     }
 
     virtual ~HealpixFourierALM()
     {
-      alloc.deallocate(alms);
+      alloc.deallocate(alms, m_size);
     }
 
-    virtual const T* data() const { return alms; }
+    virtual const std::complex<T>* data() const { return alms; }
-    virtual T * data() { return alms;} 
+    virtual std::complex<T> * data() { return alms;} 
-    virtual long size() const { return size; }
+    virtual long size() const { return m_size; }
 
     virtual FourierMap<std::complex<T> > *mimick() const
     {
@@ -127,35 +135,52 @@ namespace CosmoTool
       throw(std::bad_cast)
     {
       const HealpixFourierALM<T>& mfm = dynamic_cast<const HealpixFourierALM<T>&>(other);
+      typedef typename FourierMap<std::complex<T> >::MapType MapType;
       std::complex<T> S;
 
-      if (size != mfm.size)
+      if (m_size != mfm.m_size)
         throw std::bad_cast();
 
-      MapType m1(m_data, Npix);
+      MapType m1(alms, m_size);
-      MapType m2(mfm.m_data, mfm.Npix);
+      MapType m2(mfm.alms, mfm.m_size);
       
-      S = (m1.block(0,0,1,Lmax_+1).adjoint() * m2(0,0,1,Lmax_+1)).sum();
+      S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(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();
+      S += std::complex<T>(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<typename T> struct HealpixJobHelper__ {};
+
+  template<> struct HealpixJobHelper__<double>
+    { enum {val=1}; };
+
+  template<> struct HealpixJobHelper__<float>
+    { enum {val=0}; };
+
+
   template<typename T>
-  class HealpixFourierTransform: public FourierTransform<T>, public sharp_base
+  class HealpixFourierTransform: public FourierTransform<T>
   {
   private:
+    sharp_alm_info *ainfo;
+    sharp_geom_info *ginfo;
     HealpixFourierMap<T> realMap;
     HealpixFourierALM<T> fourierMap;
+    int m_iterate;
   public:
-    HealpixFourierTransform(long nSide, long Lmax, long Mmax)
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0)
-      : realMap(nSide), fourierMap(Lmax, Mmax)
+      : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate)
     {
-      set_Healpix_geometry(nSide);
+      sharp_make_healpix_geom_info (nSide, 1, &ginfo);
-      set_triangular_alm_info(Lmax, Mmax);
+      sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo);
     }
 
-    virtual ~HealpixFourierTransform() {}
+    virtual ~HealpixFourierTransform()
+    {
+      sharp_destroy_geom_info(ginfo);
+      sharp_destroy_alm_info(ainfo);
+    }
 
     virtual const FourierMap<std::complex<T> >& fourierSpace() const { return fourierMap; }
 
@@ -167,7 +192,7 @@ namespace CosmoTool
     
     virtual FourierTransform<T> *mimick() const 
     {
-      return new HealpixFourierTransform<T>(realMap.Nside, fourierMap.Lmax, fourierMap.Mmax);
+      return new HealpixFourierTransform<T>(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax());
     }
 
     virtual void analysis()
@@ -175,7 +200,20 @@ namespace CosmoTool
       void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
 
       sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
-        cxxjobhelper__<T>::val,0,0,0);
+        HealpixJobHelper__<T>::val,0,0,0);
+      for (int i = 0; i < m_iterate; i++)
+        {
+          HealpixFourierMap<T> tmp_map(realMap.Nside());          
+          void *tmp_ptr=reinterpret_cast<void *>(tmp_map.data());
+          typename HealpixFourierMap<T>::MapType m0 = tmp_map.eigen();
+          typename HealpixFourierMap<T>::MapType m1 = realMap.eigen();
+
+          sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+          m0 = m1 - m0;
+          sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+        }
     }
 
     virtual void synthesis()
@@ -183,30 +221,30 @@ namespace CosmoTool
       void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
 
       sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
-        cxxjobhelper__<T>::val,0,0,0);
+        HealpixJobHelper__<T>::val,0,0,0);
     }
 
     virtual void analysis_conjugate()
     {
       synthesis();
-      realMap.scale(4*M_PI/realMap.Npix());
+      realMap.scale(4*M_PI/realMap.size());
     }
 
     virtual void synthesis_conjugate()
     {
       analysis();
-      fourierMap.scale(realMap.Npix()/(4*M_PI));
+      fourierMap.scale(realMap.size()/(4*M_PI));
     }
     
   };
 
   template<typename T>
-  boost::shared_ptr<HealpixSpectrum<T>::FourierMapType> 
+  typename HealpixSpectrum<T>::ptr_map 
     HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
   { 
     const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(like_map);
     HealpixFourierALM<T> *new_alms;
-    boost::shared_ptr<FourierMapType> r(new_alms = new HealpixFourierALM<T>(alms.Lmax(), alms.Mmax()));
+    ptr_map r(new_alms = new HealpixFourierALM<T>(alms.Lmax(), alms.Mmax()));
     long lmaxGen = std::min(cls.size()-1, alms.Lmax());
     std::complex<T> *new_data = new_alms->data();