Check healpix utility calls

sample/test_healpix_calls.cpp

@@ -7,14 +7,19 @@ using namespace std;
 int main()
 {
   HealpixFourierTransform<double> dft(128,2*128,2*128, 40);
+  HealpixUtilityFunction<double> 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<double>::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace());
+
   dft.synthesis();
   cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+
   return 0;
 
 }

src/fourier/details/healpix_spectrum.hpp

@@ -13,8 +13,9 @@ namespace CosmoTool
     typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
     typedef boost::shared_ptr<FourierMapType> ptr_map;
     typedef typename SpectrumFunction<T>::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<T> *s = new HealpixSpectrum<T>(Lmax());
       s->cls = cls;
       return SpectrumFunctionPtr(s);
     }
@@ -54,13 +56,13 @@ namespace CosmoTool
     long lmaxGen = std::min(cls.size()-1, alms.Lmax());
     std::complex<T> *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<T>& c = new_data[alms.index(l,m)];
             c.real() = gsl_ran_gaussian(rng, Al); 
@@ -75,11 +77,11 @@ namespace CosmoTool
     HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
     std::complex<T> *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<typename T>
  void HealpixSpectrum<T>::mul_sqrt(FourierMapType& like_map) const
   {
-    HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(like_map);
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
     std::complex<T> *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<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
     std::complex<T> *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<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
     std::complex<T> *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;
           }

src/fourier/details/healpix_utility.hpp

@@ -11,28 +11,31 @@ namespace CosmoTool
     typedef typename MapUtilityFunction<T>::Spectrum Spectrum;
     typedef typename MapUtilityFunction<T>::Spectrum_ptr Spectrum_ptr;
     typedef typename MapUtilityFunction<T>::FMap FMap;
+    typedef typename HealpixSpectrum<T>::LType LType; 
     
     Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const
     {
       const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(m);
-      long Lmax = alms.Lmax(), Mmax = alms.Mmax();
+      LType Lmax = alms.Lmax(), Mmax = alms.Mmax();
       HealpixSpectrum<T> *spectrum = new HealpixSpectrum<T>(Lmax);
       T *data_cls = spectrum->data();
-      std::complex<T> *data_alms = alms.data();
+      const std::complex<T> *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);