Changed dependency to libsharp (MPI support, better performance). Added new power spectrum abstraction

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})

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)

src/fourier/base_types.hpp

@@ -1,12 +1,37 @@
 #ifndef __BASE_FOURIER_TYPES_HPP
 #define __BASE_FOURIER_TYPES_HPP
 
+#include <gsl/gsl_rng.h>
+#include <boost/shared_ptr.hpp>
 #include <string>
 #include <Eigen/Dense>
 #include <complex>
+#include <exception>
 
 namespace CosmoTool
 {
+  class IncompatibleMap: virtual std::exception {};
+
+  template<typename T> class FourierMap;
+
+  template<typename T>
+  class SpectrumFunction
+  {
+  protected:
+    SpectrumFunction() {}
+  public:
+    typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
+    typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
+    typedef FourierMap<std::complex<T> > FourierMapType;
+
+    virtual boost::shared_ptr<FourierMapType> 
+       newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
+
+    virtual void mul(FourierMap<std::complex<T> >& m) const = 0;
+  };
+
+
   template<typename T>
   class FourierMap
   {
@@ -14,7 +39,7 @@ namespace CosmoTool
     FourierMap() {}
     
   public:
-    typedef Eigen::Matrix<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
     typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
     typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
 
@@ -35,6 +60,12 @@ namespace CosmoTool
       return ConstMapType(data(), size());
     }
 
+    void sqrt()
+    {
+      MapType m = eigen();
+      m = m.sqrt();
+    }
+ 
     void scale(const T& factor)
     {
       MapType m(data(), size());
@@ -71,6 +102,9 @@ namespace CosmoTool
       return m;
     }
 
+    virtual T dot_product(const FourierMap<T>& second) const
+       throw(std::bad_cast) = 0;
+
     virtual FourierMap<T> *mimick() const = 0;
   };
 

src/fourier/euclidian.hpp

@@ -1,23 +1,49 @@
 #ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 #define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 
+#include <boost/function.hpp>
 #include <vector>
 #include <boost/shared_ptr.hpp>
+#include <gsl/gsl_randist.h>
 #include "base_types.hpp"
 #include "fft/fftw_calls.hpp"
+#include "../algo.hpp"
 
 namespace CosmoTool
 {
+  template<typename T>
+  class EuclidianSpectrum_1D: public SpectrumFunction<T>
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+  protected:
+    Function f;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    EuclidianSpectrum_1D(Function P)
+      : f(P)
+    {
+    }
+
+    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const;
+
+    void mul(FourierMap<std::complex<T> >& m) const;
+  };
   
   template<typename T>
-  class EuclidianFourierMap: public FourierMap<T>
+  class EuclidianFourierMapBase: public FourierMap<T>
   {
+  public:
+    typedef std::vector<int> DimArray;
   private:
     boost::shared_ptr<T> m_data;
-    std::vector<int> m_dims;
+    DimArray m_dims;
     long m_size;
   public:
-    EuclidianFourierMap(boost::shared_ptr<T> indata, std::vector<int> indims)
+    
+    EuclidianFourierMapBase(boost::shared_ptr<T> 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,32 +69,149 @@ namespace CosmoTool
       return m;
     }
 
+  };
+
+  template<typename T>
+  class EuclidianFourierMapReal: public EuclidianFourierMapBase<T>
+  {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
+
+    EuclidianFourierMapReal(boost::shared_ptr<T> indata, const DimArray& indims)
+      : EuclidianFourierMapBase<T>(indata, indims) 
+    {}
+
     virtual FourierMap<T> *mimick() const
     {
-      return new EuclidianFourierMap<T>(
+      return new EuclidianFourierMapReal<T>(
-		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*m_size), 
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*this->size()), 
 					   std::ptr_fun(fftw_free)),
-		       m_dims);
+		      this->getDims());
+    }
+
+    virtual T dot_product(const FourierMap<T>& other) const
+      throw(std::bad_cast)
+    {
+      const EuclidianFourierMapReal<T>& m2 = dynamic_cast<const EuclidianFourierMapReal<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+
+      return (this->eigen()*m2.eigen()).sum();
     }
   };
 
+  template<typename T>
+  class EuclidianFourierMapComplex: public EuclidianFourierMapBase<std::complex<T> >
+  {
+  protected:
+    typedef boost::shared_ptr<std::complex<T> > ptr_t;
+    std::vector<double> delta_k;
+    long plane_size;
+  public:
+    typedef typename EuclidianFourierMapBase<std::complex<T> >::DimArray DimArray;
+
+    EuclidianFourierMapComplex(ptr_t indata,
+			       const DimArray& indims, 
+			       const std::vector<double>& dk)
+      : EuclidianFourierMapBase<std::complex<T> >(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<std::complex<T> > *mimick() const
+    {
+      return 
+	new EuclidianFourierMapComplex<T>(
+		   ptr_t((std::complex<T> *)
+			 fftw_malloc(sizeof(std::complex<T>)*this->size()), 
+			 std::ptr_fun(fftw_free)),
+		   this->getDims(),
+		   this->delta_k);
+    }
+
+    template<typename Array>
+    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<T> dot_product(const FourierMap<std::complex<T> >& other) const
+      throw(std::bad_cast)
+    {
+      const EuclidianFourierMapComplex<T>& m2 = dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+      
+      const std::complex<T> *d1 = this->data();
+      const std::complex<T> *d2 = m2.data();
+      const DimArray& dims = this->getDims();
+      int N0 = dims[0];
+      std::complex<T> 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<typename T>
   class EuclidianFourierTransform: public FourierTransform<T>
   {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
   private:
     typedef FFTW_Calls<T> calls;
-    EuclidianFourierMap<T> *realMap;
+    EuclidianFourierMapReal<T> *realMap;
-    EuclidianFourierMap<std::complex<T> > *fourierMap;
+    EuclidianFourierMapComplex<T> *fourierMap;
     typename calls::plan_type m_analysis, m_synthesis;
     double volume;
     long N;
-    std::vector<int> m_dims;
+    DimArray m_dims;
     std::vector<double> m_L;
   public:
-    EuclidianFourierTransform(const std::vector<int>& dims, const std::vector<double>& L)
+    EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L)
     {
       assert(L.size() == dims.size());
+      std::vector<double> 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<T>(
+      realMap = new EuclidianFourierMapReal<T>(
 		   boost::shared_ptr<T>(calls::alloc_real(N),
 					std::ptr_fun(calls::free)),
 		   dims);
-      fourierMap = new EuclidianFourierMap<std::complex<T> >(
+      fourierMap = new EuclidianFourierMapComplex<T>(
 		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)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 T>
+  typename EuclidianSpectrum_1D<T>::ptr_map
+   EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
+  {
+    typedef EuclidianFourierMapComplex<T> MapT;
+    typedef typename MapT::DimArray DimArray;
 
+    MapT& m_c = dynamic_cast<MapT&>(like_map);
+    MapT *rand_map = m_c.mimick();
+    std::complex<T> *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<T>(gsl_ran_gaussian(rng, A_k),
+			       gsl_ran_gaussian(rng, A_k));
+      }
+    // Generate the mean value
+    d[0] = std::complex<T>(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<T>(d[idx].real() + d[idx].imag(), 0);
+	return boost::shared_ptr<EuclidianSpectrum_1D<T>::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<T>(d[q].real() + d[q].imag());
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul(FourierMap<std::complex<T> >& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+      d[p] *= f(m.get_K(p));
+  }
 };
 
 #endif

src/fourier/healpix.hpp

@@ -1,24 +1,45 @@
 #ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
 #define __COSMOTOOL_FOURIER_HEALPIX_HPP
 
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <cmath>
 #include <boost/bind.hpp>
 #include <boost/shared_ptr.hpp>
-#include <alm.h>
+#include <exception>
-#include <healpix_base.h>
+#include <sharp_cxx.h>
-#include <psht_cxx.h>
 
 namespace CosmoTool
 {
+  template<typename T>
+  class HealpixSpectrum: public FourierSpectrum<T>
+  {
+  protected:
+    std::vector<T> cls;
+  public:
+    typedef typename FourierSpectrum<T>::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<FourierMapType> 
+       newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
+  };    
 
   template<typename T>
-  class HealpixFourierMap: public FourierMap<T>, public Healpix_Base
+  class HealpixFourierMap: public FourierMap<T>
   {
   private:
     T *m_data;
+    long Npix, Nside;
     Eigen::aligned_allocator<T> 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<T>& other) const
+      throw(std::bad_cast)
+    {
+      const HealpixFourierMap<T>& mfm = dynamic_cast<const HealpixFourierMap<T>&>(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<T> *mimick() const
     {
-      return new HealpixFourierMap<T>(Nside());
+      return new HealpixFourierMap<T>(Nside);
     }
   };
 
 
   template<typename T>
-  class HealpixFourierALM: public FourierMap<std::complex<T> >, public Alm_Base
+  class HealpixFourierALM: public FourierMap<std::complex<T> >
   {
   private:
     std::complex<T> *alms;
     long size;
+    long Lmax_, Mmax_, TVal_;
     Eigen::aligned_allocator<std::complex<T> > alloc;
   public:
-    HealpixFourierALM(long Lmax, long Mmax)
+    typedef unsigned long LType;
-      : Alm_Base(Lmax, Mmax)
+
+    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<T> *mimick() const
+    virtual FourierMap<std::complex<T> > *mimick() const
     {
-      return new HealpixFourierALM<T>(Lmax(), Mmax());
+      return new HealpixFourierALM<T>(Lmax_, Mmax_);
+    }
+
+    virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
+      throw(std::bad_cast)
+    {
+      const HealpixFourierALM<T>& mfm = dynamic_cast<const HealpixFourierALM<T>&>(other);
+      std::complex<T> 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<typename T>
-  class HealpixFourierTransform: public FourierTransform<T>
+  class HealpixFourierTransform: public FourierTransform<T>, public sharp_base
   {
   private:
     HealpixFourierMap<T> realMap;
     HealpixFourierALM<T> fourierMap;
-    psht_joblist<T> jobs;
   public:
     HealpixFourierTransform(long nSide, long Lmax, long Mmax)
       : realMap(nSide), fourierMap(Lmax, Mmax)
     {
-      jobs.set_Healpix_geometry(nSide);
+      set_Healpix_geometry(nSide);
-      jobs.set_triangular_alm_info(Lmax, Mmax);
+      set_triangular_alm_info(Lmax, Mmax);
     }
 
     virtual ~HealpixFourierTransform() {}
@@ -96,25 +167,23 @@ 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()
     {
-      jobs.add_map2alm(realMap.data(), 
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
-		       reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+
-		       false);
+      sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
-      jobs.execute();
+        cxxjobhelper__<T>::val,0,0,0);
-      jobs.clear_jobs();
     }
 
     virtual void synthesis()
     {
-      jobs.add_alm2map(reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
-		       realMap.data(),
+
-		       false);
+      sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
-      jobs.execute();
+        cxxjobhelper__<T>::val,0,0,0);
-      jobs.clear_jobs();
     }
 
     virtual void analysis_conjugate()
@@ -130,6 +199,32 @@ namespace CosmoTool
     }
     
   };
+
+  template<typename T>
+  boost::shared_ptr<HealpixSpectrum<T>::FourierMapType> 
+    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()));
+    long lmaxGen = std::min(cls.size()-1, alms.Lmax());
+    std::complex<T> *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<T>& c = new_data[alms.index(l,m)];
+            c.real() = gsl_ran_gaussian(rng, Al); 
+            c.imag() = gsl_ran_gaussian(rng, Al); 
+          }
+      }
+    return r;
+  } 
 };
 
 #endif