From 4c482daf9b8a1ee5e69f7675a361aae0c7c93878 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Sun, 18 Nov 2012 16:18:18 -0500
Subject: [PATCH] Fixed dof interface. Implement utility class for healpix

---
 src/fourier/base_types.hpp               |  2 +-
 src/fourier/details/healpix_spectrum.hpp | 11 ++++-
 src/fourier/details/healpix_utility.hpp  | 61 ++++++++++++++++++++++++
 src/fourier/healpix.hpp                  |  2 +-
 4 files changed, 73 insertions(+), 3 deletions(-)
 create mode 100644 src/fourier/details/healpix_utility.hpp

diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp
index ae62f54..5e93e48 100644
--- a/src/fourier/base_types.hpp
+++ b/src/fourier/base_types.hpp
@@ -35,7 +35,7 @@ namespace CosmoTool
     virtual SpectrumFunctionPtr copy() const = 0;
 
     virtual const T *data() const { return 0; }
-    virtual const T *dof() const { return 0; }
+    virtual const int *dof() const { return 0; }
     virtual long size() const { return -1; }
 
     virtual void sqrt() = 0;
diff --git a/src/fourier/details/healpix_spectrum.hpp b/src/fourier/details/healpix_spectrum.hpp
index 2031acf..e797c91 100644
--- a/src/fourier/details/healpix_spectrum.hpp
+++ b/src/fourier/details/healpix_spectrum.hpp
@@ -8,16 +8,25 @@ namespace CosmoTool
   {
   protected:
     std::vector<T> cls;
+    int *m_dof;
   public:
     typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
     typedef boost::shared_ptr<FourierMapType> ptr_map;
     typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
 
     HealpixSpectrum(long Lmax)
-     : cls(Lmax+1) {}
+      : cls(Lmax+1), m_dof(new int[Lmax+1]) 
+    {
+      for (int l = 0; l <= Lmax; l++)
+	m_dof[l] = l + 1;
+    }
 
     T *data() { return &cls[0]; }
     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; }
+
     long size() const { return cls.size(); }
 
     void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; 
diff --git a/src/fourier/details/healpix_utility.hpp b/src/fourier/details/healpix_utility.hpp
new file mode 100644
index 0000000..b527366
--- /dev/null
+++ b/src/fourier/details/healpix_utility.hpp
@@ -0,0 +1,61 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+
+namespace CosmoTool
+{
+
+  template<typename T>
+  class HealpixUtilityFunction: public MapUtilityFunction<T>
+  {
+  public:
+    typedef typename MapUtilityFunction<T>::Spectrum Spectrum;
+    typedef typename MapUtilityFunction<T>::Spectrum_ptr Spectrum_ptr;
+    typedef typename MapUtilityFunction<T>::FMap FMap;
+
+    Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const
+    {
+      const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(m);
+      long 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();
+
+      // make an estimate of the cls
+      std::fill(data_cls, data_cls+Lmax+1, 0);
+      for (long m = 0, q = 0; m <= Mmax; ++m )
+	{
+	  for (long l = 0; l <= Lmax; ++l)
+	    {
+	      // Triangular storage
+	      data_cls[l] += std::norm(alms[l+q]);
+	    }
+	  q += Lmax+1;
+	}
+      
+      for (long l = 0; l <= Lmax; ++l)
+	{
+	  int dof = 1 + std::min(l, Mmax);
+	  spectrum->set_dof(l, dof);
+	  data_cls[l] /= dof;
+	}
+
+      return Spectrum_ptr(spectrum);
+    }
+
+    Spectrum_ptr newSpectrumFromRaw(T *data, long size,
+				    Spectrum_ptr like_spec) const
+    {
+      Spectrum *s = like_spec.get();
+      HealpixSpectrum<T>& in_spec = dynamic_cast<HealpixSpectrum<T>&>(*s);
+      HealpixSpectrum<T> *new_spectrum = new HealpixSpectrum<T>(in_spec.Lmax());
+      T *out_d = new_spectrum->data();
+
+      std::copy(data, data + min(size,new_spectrum->size()), out_d);
+
+      return Spectrum_ptr(new_spectrum);
+    }
+  };
+
+};
+
+#endif
diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp
index 46f5356..36f7979 100644
--- a/src/fourier/healpix.hpp
+++ b/src/fourier/healpix.hpp
@@ -19,6 +19,6 @@
 #include "details/healpix_alms.hpp"
 #include "details/healpix_transform.hpp"
 #include "details/healpix_spectrum.hpp"
-
+#include "details/healpix_utility.hpp"
 
 #endif