Implemented power spectrum support in Euclidian/Healpix

sample/test_fft_calls.cpp

@@ -1,13 +1,23 @@
+#include "yorick.hpp"
+#include <gsl/gsl_rng.h>
 #include <iostream>
+#include <cmath>
 #include "fourier/euclidian.hpp"
 
 using namespace CosmoTool;
 using namespace std;
 
+double spectrum_generator(double k)
+{
+  return 1/(0.1+pow(k, 3.0));
+}
+
 int main()
 {
   EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+  EuclidianSpectrum_1D<double> spectrum(spectrum_generator);
   double volume = 128*128;
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
 
   dft.realSpace().eigen().setRandom();
   dft.analysis();
@@ -15,6 +25,17 @@ int main()
   cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl;
   dft.synthesis();
   cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+
+  dft.realSpace().scale(2.0);
+  dft.fourierSpace().scale(0.2);
+
+  SpectrumFunction<double>::FourierMapPtr m = spectrum.newRandomFourier(rng, dft.fourierSpace()); 
+  dft.fourierSpace() = *m.get();
+  dft.synthesis();
+
+  uint32_t dims[2] = { 128, 128 };
+  CosmoTool::saveArray("generated_map.nc", dft.realSpace().data(), dims, 2);
+  
   return 0;
 
 }

sample/test_healpix_calls.cpp

@@ -6,7 +6,7 @@ using namespace std;
 
 int main()
 {
-  HealpixFourierTransform<double> dft(128,3*128,3*128, 40);
+  HealpixFourierTransform<double> dft(128,2*128,2*128, 40);
   long Npix = dft.realSpace().size();
 
   dft.realSpace().eigen().setRandom();

src/fourier/base_types.hpp

@@ -24,14 +24,17 @@ namespace CosmoTool
     typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
     typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
     typedef FourierMap<std::complex<T> > FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> FourierMapPtr;
 
-    virtual boost::shared_ptr<FourierMapType> 
+    virtual FourierMapPtr 
        newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
 
-    virtual void mul(FourierMap<std::complex<T> >& m) const = 0;
+    virtual void mul(FourierMapType& m) const = 0;
+    virtual void mul_sqrt(FourierMapType& m) const = 0;
+    virtual void mul_inv(FourierMapType& m) const = 0;
+    virtual void mul_inv_sqrt(FourierMapType& m) const = 0;
   };
 
-
   template<typename T>
   class FourierMap
   {
@@ -60,6 +63,12 @@ namespace CosmoTool
       return ConstMapType(data(), size());
     }
 
+    FourierMap<T>& operator=(const FourierMap<T>& m)
+    {
+      eigen() = m.eigen();
+      return *this;
+    }
+
     void sqrt()
     {
       MapType m = eigen();
@@ -130,6 +139,14 @@ namespace CosmoTool
     virtual void synthesis_conjugate() = 0;
   };
 
+  template<typename T>
+  class MapUtilityFunction
+  {
+    typedef SpectrumFunction<T> Spectrum;
+    typedef Spectrum *Spectrum_ptr;
+    typedef FourierMap<std::complex<T> > FMap;
+    typedef Spectrum_ptr (*SpectrumFromMap)(const FMap& m);
+  };
 
 };
 

src/fourier/euclidian.hpp

@@ -30,6 +30,9 @@ namespace CosmoTool
     ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const;
 
     void mul(FourierMap<std::complex<T> >& m) const;
+    void mul_sqrt(FourierMap<std::complex<T> >& m) const;
+    void mul_inv(FourierMap<std::complex<T> >& m) const;
+    void mul_inv_sqrt(FourierMap<std::complex<T> >& m) const;
   };
   
   template<typename T>
@@ -140,7 +143,7 @@ namespace CosmoTool
       double k2 = 0;
       for (int q = 0; q < ik.size(); q++)
 	{
-	  int dk = ik;
+	  int dk = ik[q];
 
 	  if (dk > dims[q]/2)
 	    dk = dk - dims[q];
@@ -349,18 +352,21 @@ namespace CosmoTool
    EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
   {
     typedef EuclidianFourierMapComplex<T> MapT;
+    typedef typename EuclidianSpectrum_1D<T>::ptr_map ptr_map;
     typedef typename MapT::DimArray DimArray;
 
-    MapT& m_c = dynamic_cast<MapT&>(like_map);
+    const MapT& m_c = dynamic_cast<const MapT&>(like_map);
-    MapT *rand_map = m_c.mimick();
+    ptr_map ret_map = ptr_map(m_c.mimick());
-    std::complex<T> *d = rand_map->data();
+    MapT& rand_map = dynamic_cast<MapT&>(*ret_map.get());
+
+    std::complex<T> *d = rand_map.data();
     long idx;
-    const DimArray& dims = rand_map->getDims();
+    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)));
+	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));
       }
@@ -372,7 +378,7 @@ namespace CosmoTool
     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);
+        return ret_map;
       }
 
     plane_size = 1;
@@ -391,6 +397,7 @@ namespace CosmoTool
       }
     long q = dims[0];
     d[q] = std::complex<T>(d[q].real() + d[q].imag());
+    return ret_map;
   }
 
   template<typename T>
@@ -400,8 +407,51 @@ namespace CosmoTool
     std::complex<T> *d = m.data();
 
     for (long p = 0; p < m_c.size(); p++)
-      d[p] *= f(m.get_K(p));
+      d[p] *= f(m_c.get_K(p));
   }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_sqrt(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] *= std::sqrt(f(m_c.get_K(p)));
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv(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++)
+     {
+        T A = f(m_c.get_K(p));
+        if (A==0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+     }
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv_sqrt(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++)
+      {
+        T A = std::sqrt(f(m_c.get_K(p)));
+        if (A == 0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+      }
+  }
+
 };
 
 #endif

src/fourier/healpix.hpp

@@ -32,7 +32,12 @@ namespace CosmoTool
     const T *data() const { return &cls[0]; }
     long size() const { return cls.size(); }
 
-    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
+    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; 
+
+    void mul(FourierMapType& m) const;
+    void mul_sqrt(FourierMapType& m) const;
+    void mul_inv(FourierMapType& m) const;
+    void mul_inv_sqrt(FourierMapType& m) const;
   };
 
   template<typename T>
@@ -248,13 +253,13 @@ namespace CosmoTool
     long lmaxGen = std::min(cls.size()-1, alms.Lmax());
     std::complex<T> *new_data = new_alms->data();
 
-    for (long l = 0; l < lmaxGen; l++)
+    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++)
+        for (long 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); 
@@ -263,6 +268,75 @@ namespace CosmoTool
       }
     return r;
   } 
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (long l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = cls[l];
+
+        for (long m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (long l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = std::sqrt(cls[l]);
+
+        for (long m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (long 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++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (long 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++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
 };
 
 #endif