Updated again cosmotoolbox

external/cosmotool/sample/test_cosmopower.cpp

@@ -0,0 +1,29 @@
+#include <cmath>
+#include <iostream>
+#include "cosmopower.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+int main()
+{
+   CosmoPower cp;
+   CosmoPower::CosmoFunction f[] = { CosmoPower::POWER_EFSTATHIOU, CosmoPower::HU_WIGGLES, CosmoPower::HU_BARYON, CosmoPower::POWER_SUGIYAMA };
+   int num_F = sizeof(f)/sizeof(f[0]);
+
+   cp.setFunction(f[0]);
+   cp.normalize();
+   for (int ik = 0; ik < 100; ik++)
+    {
+      double k = pow(10.0, 4*ik/100.-3);
+      
+      cout << k << " ";
+      for (int q = 0; q < num_F; q++)
+        {
+          cp.setFunction(f[q]);
+          cout << cp.power(k) << " ";
+        } 
+      cout << endl;
+    }
+   return 0;
+}

external/cosmotool/sample/test_fft_calls.cpp

@@ -1,12 +1,58 @@
+#include <boost/format.hpp>
+#include "yorick.hpp"
+#include <gsl/gsl_rng.h>
+#include <iostream>
+#include <cmath>
 #include "fourier/euclidian.hpp"
 
 using namespace CosmoTool;
+using boost::format;
+using boost::str;
+using namespace std;
 
-int main()
+double spectrum_generator(double k)
 {
-  EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+  if (k==0)
+    return 0;
+  return 1/(0.01+pow(k, 3.0));
+}
+
+template<typename T>
+void test_2d(int Nx, int Ny)
+{
+  EuclidianFourierTransform_2d<T> dft(Nx,Ny,1.0,1.0);
+  EuclidianSpectrum_1D<T> spectrum(spectrum_generator);
+  double volume = Nx*Ny;
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
 
   dft.realSpace().eigen().setRandom();
+  dft.fourierSpace().scale(std::complex<T>(0,0));
   dft.analysis();
+  cout << format("Testing (%d,%d)") % Nx % Ny << endl;
+  cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+  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);
+
+  spectrum.newRandomFourier(rng, dft.fourierSpace()); 
+  dft.synthesis();
+
+  uint32_t dims[2] = { Ny, Nx };
+  CosmoTool::saveArray(str(format("generated_map_%d_%d.nc") %Nx % Ny) , dft.realSpace().data(), dims, 2);
+
+  gsl_rng_free(rng);
+}
+
+int main(int argc, char **argv)
+{
+  test_2d<double>(128,128);
+  test_2d<double>(131,128);  
+  test_2d<double>(130,128);  
+  test_2d<float>(128,128);
+  test_2d<float>(128,131);  
+  test_2d<float>(128,130);  
   return 0;
 }

external/cosmotool/src/fourier/base_types.hpp

@@ -1,12 +1,51 @@
 #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;
+    typedef boost::shared_ptr<FourierMapType> FourierMapPtr;
+    typedef boost::shared_ptr<SpectrumFunction<T> > SpectrumFunctionPtr;
+
+    virtual ~SpectrumFunction() {}
+
+    virtual void 
+       newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0;
+
+    virtual SpectrumFunctionPtr copy() const = 0;
+
+    virtual const T *data() const { return 0; }
+    virtual const int *dof() const { return 0; }
+    virtual long size() const { return -1; }
+
+    virtual void sqrt() = 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
   {
@@ -14,7 +53,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;
 
@@ -34,6 +73,18 @@ namespace CosmoTool
     {
       return ConstMapType(data(), size());
     }
+
+    FourierMap<T>& operator=(const FourierMap<T>& m)
+    {
+      eigen() = m.eigen();
+      return *this;
+    }
+
+    void sqrt()
+    {
+      MapType m = eigen();
+      m = m.sqrt();
+    }
  
     void scale(const T& factor)
     {
@@ -45,7 +96,7 @@ namespace CosmoTool
     {
       assert(size() == map2->size());
       MapType m(data(), size());
-      MapType m2(map2->data(), map2->size());
+      ConstMapType m2(map2->data(), map2->size());
       m *= m2;
     }
 
@@ -71,6 +122,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;
   };
 
@@ -96,6 +150,18 @@ namespace CosmoTool
     virtual void synthesis_conjugate() = 0;
   };
 
+  template<typename T>
+  class MapUtilityFunction
+  {
+  public:
+    typedef SpectrumFunction<T> Spectrum;
+    typedef boost::shared_ptr<Spectrum> Spectrum_ptr;
+    typedef FourierMap<std::complex<T> > FMap;
+
+    virtual Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const = 0;
+    virtual Spectrum_ptr newSpectrumFromRaw(T *data, long size,
+                                            const Spectrum& like_spec) const = 0;
+  };
 
 };
 

external/cosmotool/src/fourier/euclidian.hpp

@@ -1,200 +1,17 @@
 #ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 #define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 
+#include <cmath>
+#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"
-
-namespace CosmoTool
-{
-  
-  template<typename T>
-  class EuclidianFourierMap: public FourierMap<T>
-  {
-  private:
-    boost::shared_ptr<T> m_data;
-    std::vector<int> m_dims;
-    long m_size;
-  public:
-    EuclidianFourierMap(boost::shared_ptr<T> indata, std::vector<int> indims)
-    {
-      m_data = indata;
-      m_dims = indims;
-      m_size = 1;
-      for (int i = 0; i < m_dims.size(); i++)
-	m_size *= m_dims[i];
-    }
-
-    virtual ~EuclidianFourierMap()
-    {
-    }
-
-    virtual const T *data() const { return m_data.get(); }
-    virtual T *data() { return m_data.get(); }
-    virtual long size() const  { return m_size; } 
-    
-    virtual FourierMap<T> *copy() const
-    {
-      FourierMap<T> *m = this->mimick();
-      m->eigen() = this->eigen();
-      return m;
-    }
-
-    virtual FourierMap<T> *mimick() const
-    {
-      return new EuclidianFourierMap<T>(
-		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*m_size), 
-					   std::ptr_fun(fftw_free)),
-		       m_dims);
-    }
-  };
-
-  
-  template<typename T>
-  class EuclidianFourierTransform: public FourierTransform<T>
-  {
-  private:
-    typedef FFTW_Calls<T> calls;
-    EuclidianFourierMap<T> *realMap;
-    EuclidianFourierMap<std::complex<T> > *fourierMap;
-    typename calls::plan_type m_analysis, m_synthesis;
-    double volume;
-    long N;
-    std::vector<int> m_dims;
-    std::vector<double> m_L;
-  public:
-    EuclidianFourierTransform(const std::vector<int>& dims, const std::vector<double>& L)
-    {
-      assert(L.size() == dims.size());
-
-      m_dims = dims;
-      m_L = L;
-
-      N = 1;
-      volume = 1;
-      for (int i = 0; i < dims.size(); i++)
-	{
-	  N *= dims[i];
-	  volume *= L[i];
-	}
-
-      realMap = new EuclidianFourierMap<T>(
-		   boost::shared_ptr<T>(calls::alloc_real(N),
-					std::ptr_fun(calls::free)),
-		   dims);
-      fourierMap = new EuclidianFourierMap<std::complex<T> >(
-		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
-						       std::ptr_fun(calls::free)), 
-		   dims);
-      m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
-				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
-				       FFTW_MEASURE);
-      m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0],
-					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
-					FFTW_MEASURE);
-    }
-    
-    virtual ~EuclidianFourierTransform()
-    {
-      delete realMap;
-      delete fourierMap;
-      calls::destroy_plan(m_synthesis);
-      calls::destroy_plan(m_analysis);
-    }
-
-    void synthesis()
-    {
-      calls::execute(m_synthesis);
-      realMap->scale(1/volume);
-    }
-    
-    void analysis()
-    {
-      calls::execute(m_analysis);
-      fourierMap->scale(volume/N);
-    }
-    
-    void synthesis_conjugate()
-    {
-      calls::execute(m_analysis);
-      fourierMap->scale(1/volume);
-    }
-    
-    void analysis_conjugate()
-    {
-      calls::execute(m_synthesis);
-      realMap->scale(volume/N);
-    }
-
-    const FourierMap<std::complex<T> >& fourierSpace() const
-    {
-      return *fourierMap;
-    }
-
-    FourierMap<std::complex<T> >& fourierSpace()
-    {
-      return *fourierMap;
-    }
-
-    const FourierMap<T>& realSpace() const
-    {
-      return *realMap;
-    }
-
-    FourierMap<T>& realSpace()
-    {
-      return *realMap;
-    }
-
-    FourierTransform<T> *mimick() const
-    {
-      return new EuclidianFourierTransform(m_dims, m_L);
-    }
-  };
-
-  template<typename T>
-  class EuclidianFourierTransform_2d: public EuclidianFourierTransform<T>
-  {
-  private:
-    template<typename T2>
-    static std::vector<T2> make_2d_vector(T2 a, T2 b)
-    {
-      T2 arr[2] = { a, b};
-      return std::vector<T2>(&arr[0], &arr[2]);
-    }
-  public:
-    EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly)
-      : EuclidianFourierTransform<T>(make_2d_vector<int>(Nx, Ny), make_2d_vector<double>(Lx, Ly))
-    {
-    }
-
-    virtual ~EuclidianFourierTransform_2d() {}
-
-  };
-
-  template<typename T>
-  class EuclidianFourierTransform_3d: public EuclidianFourierTransform<T>
-  {
-  private:
-    template<typename T2>
-    static std::vector<T2> make_3d_vector(T2 a, T2 b, T2 c)
-    {
-      T2 arr[2] = { a, b, c};
-      return std::vector<T2>(&arr[0], &arr[3]);
-    }
-
-  public:
-    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
-      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
-    {
-    }
-
-    virtual ~EuclidianFourierTransform_3d() {}
-  };
-
-
-
-};
+#include "../algo.hpp"
+#include "details/euclidian_maps.hpp"
+#include "details/euclidian_spectrum_1d.hpp"
+#include "details/euclidian_spectrum_1d_bin.hpp"
+#include "details/euclidian_transform.hpp"
 
 #endif

external/cosmotool/src/fourier/fft/fftw_calls.hpp

@@ -28,8 +28,8 @@ public: \
   typedef prefix ## _complex complex_type; \
   typedef prefix ## _plan plan_type; \
   \
-  static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \
-  static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \
+  static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \
   static void free(void *p) { fftw_free(p); } \
 \
   static void execute(plan_type p) { prefix ## _execute(p); } \

external/cosmotool/src/fourier/healpix.hpp

@@ -1,135 +1,24 @@
 #ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
 #define __COSMOTOOL_FOURIER_HEALPIX_HPP
 
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <cmath>
+#include <vector>
 #include <boost/bind.hpp>
 #include <boost/shared_ptr.hpp>
-#include <alm.h>
-#include <healpix_base.h>
-#include <psht_cxx.h>
+#include <exception>
+#include "base_types.hpp"
+#include <sharp_lowlevel.h>
+#include <sharp_geomhelpers.h>
+#include <sharp_almhelpers.h>
+#include <algorithm>
+#include <functional>
 
-namespace CosmoTool
-{
-    
-  template<typename T>
-  class HealpixFourierMap: public FourierMap<T>, public Healpix_Base
-  {
-  private:
-    T *m_data;
-    Eigen::aligned_allocator<T> alloc;
-  public:
-    HealpixFourierMap(long nSide)
-      : Healpix_Base(RING, nSide, SET_NSIDE)
-    {
-      m_data = alloc.allocate(Npix);
-    }
-
-    virtual ~HealpixFourierMap()
-    {
-      alloc.deallocate(m_data);
-    }
-
-    virtual const T* data() const { return m_data; }
-    virtual T *data() { return m_data; }
-    virtual long size() const { return Npix(); }
-    
-    virtual FourierMap<T> *mimick() const
-    {
-      return new HealpixFourierMap<T>(Nside());
-    }
-  };
-
-
-  template<typename T>
-  class HealpixFourierALM: public FourierMap<std::complex<T> >, public Alm_Base
-  {
-  private:
-    std::complex<T> *alms;
-    long size;
-    Eigen::aligned_allocator<std::complex<T> > alloc;
-  public:
-    HealpixFourierALM(long Lmax, long Mmax)
-      : Alm_Base(Lmax, Mmax)
-    {
-      size = Num_Alms(Lmax, Mmax);
-      alms = alloc.allocate(size);
-    }
-
-    virtual ~HealpixFourierALM()
-    {
-      alloc.deallocate(alms);
-    }
-
-    virtual const T* data() const { return alms; }
-    virtual T * data() { return alms;} 
-    virtual long size() const { return size; }
-
-    virtual FourierMap<T> *mimick() const
-    {
-      return new HealpixFourierALM<T>(Lmax(), Mmax());
-    }
-  };
-
-  template<typename T>
-  class HealpixFourierTransform: public FourierTransform<T>
-  {
-  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);
-      jobs.set_triangular_alm_info(Lmax, Mmax);
-    }
-
-    virtual ~HealpixFourierTransform() {}
-
-    virtual const FourierMap<std::complex<T> >& fourierSpace() const { return fourierMap; }
-
-    virtual FourierMap<std::complex<T> >& fourierSpace() { return fourierMap; }
-
-    virtual const FourierMap<T>& realSpace() const { return realMap; }
-
-    virtual FourierMap<T>& realSpace() { return realMap; }
-    
-    virtual FourierTransform<T> *mimick() const 
-    {
-      return new HealpixFourierTransform<T>(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax());
-    }
-
-    virtual void analysis()
-    {
-      jobs.add_map2alm(realMap.data(), 
-		       reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
-		       false);
-      jobs.execute();
-      jobs.clear_jobs();
-    }
-
-    virtual void synthesis()
-    {
-      jobs.add_alm2map(reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
-		       realMap.data(),
-		       false);
-      jobs.execute();
-      jobs.clear_jobs();
-    }
-
-    virtual void analysis_conjugate()
-    {
-      synthesis();
-      realMap.scale(4*M_PI/realMap.Npix());
-    }
-
-    virtual void synthesis_conjugate()
-    {
-      analysis();
-      fourierMap.scale(realMap.Npix()/(4*M_PI));
-    }
-    
-  };
-};
+#include "details/healpix_map.hpp"
+#include "details/healpix_alms.hpp"
+#include "details/healpix_transform.hpp"
+#include "details/healpix_spectrum.hpp"
+#include "details/healpix_utility.hpp"
 
 #endif

external/cosmotool/src/h5_readFlash.cpp

@@ -185,7 +185,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id)
+                        long    *id)
 {
 
   herr_t   status;
@@ -341,7 +341,7 @@ void h5_read_flash3_particles (H5File* file,
 
     if (id) {
     for(p=0; p < (pcount); p++) {
-      id[p+poffset] = (int)  *(partBuffer+iptag-1+p*numProps);
+      id[p+poffset] = (long)  *(partBuffer+iptag-1+p*numProps);
     } }
 
     if (pos1 && pos2 && pos3) {

external/cosmotool/src/h5_readFlash.hpp

@@ -29,7 +29,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id);
+                        long    *id);
 
 void h5_read_flash3_header_info(H5File* file,
 				double* time,                   /* simulation time */

external/cosmotool/src/loadFlash.cpp

@@ -71,7 +71,7 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
     } }
  
     if (loadflags & NEED_GADGET_ID) {
-    data->Id = new int[data->NumPart];
+    data->Id = new long[data->NumPart];
     if (data->Id == 0) {
       delete data;
       return 0;

external/cosmotool/src/loadGadget.cpp

@@ -195,7 +195,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
     try
       {
 	f->beginCheckpoint();
-	data->Id = new int[data->NumPart];
+	data->Id = new long[data->NumPart];
 	if (data->Id == 0)
 	  {
 	    delete f;

external/cosmotool/src/loadRamses.cpp

@@ -384,7 +384,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
 	}
       if (flags & NEED_GADGET_ID)
 	{
-	  data->Id = new int[nPar];
+	  data->Id = new long[nPar];
 	}
 
       for (int k = 0; k < 3; k++)

external/cosmotool/src/loadSimu.hpp

@@ -38,7 +38,7 @@ namespace CosmoTool
 
     long NumPart;
     long TotalNumPart;
-    int *Id;
+    long *Id;
     float *Pos[3];
     float *Vel[3];
     int *type;

external/cosmotool/src/powerSpectrum.cpp

@@ -1,3 +1,4 @@
+#include <cassert>
 #include <vector>
 #include <algorithm>
 #include <iostream>
@@ -21,7 +22,7 @@ using namespace std;
 #define POWER_BDM 7
 #define POWER_TEST 8
 
-#define POWER_SPECTRUM HU_WIGGLES
+#define POWER_SPECTRUM POWER_EFSTATHIOU
 
 namespace Cosmology {
 
@@ -651,4 +652,17 @@ double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2)
 #endif
   }
 
+  double vvCorrection(double P_deltadelta, double k)
+  {
+    static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796;
+    if (k > 0.3)
+      return 0;
+    double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta);
+    assert(P_deltadelta > 0);
+
+    if (r < 0)
+      return 0;
+    return r;
+  }
+
 };

external/cosmotool/src/powerSpectrum.hpp

@@ -30,6 +30,8 @@ namespace Cosmology {
   double computeVarianceZero(double powNorm);
   double computeCorrel(double powNorm,  double topHatRad1);
   double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2);
+
+  double vvCorrection(double P_deltadelta, double k);
 };
 
 #endif