#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 <exception>
#include <sharp_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>
  {
  private:
    T *m_data;
    long Npix, Nside;
    Eigen::aligned_allocator<T> alloc;
  public:
    HealpixFourierMap(long nSide)
      : Npix(12*nSide*nSide), Nside(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 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);
    }
  };


  template<typename T>
  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:
    typedef unsigned long LType;

    LType Lmax() const { return Lmax_; }
    LType Mmax() const { return Mmax_; }

    LType Num_Alms() const
    {
      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);
    }

    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<std::complex<T> > *mimick() const
    {
      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>, public sharp_base
  {
  private:
    HealpixFourierMap<T> realMap;
    HealpixFourierALM<T> fourierMap;
  public:
    HealpixFourierTransform(long nSide, long Lmax, long Mmax)
      : realMap(nSide), fourierMap(Lmax, Mmax)
    {
      set_Healpix_geometry(nSide);
      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()
    {
      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());

      sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
        cxxjobhelper__<T>::val,0,0,0);
    }

    virtual void synthesis()
    {
      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());

      sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
        cxxjobhelper__<T>::val,0,0,0);
    }

    virtual void analysis_conjugate()
    {
      synthesis();
      realMap.scale(4*M_PI/realMap.Npix());
    }

    virtual void synthesis_conjugate()
    {
      analysis();
      fourierMap.scale(realMap.Npix()/(4*M_PI));
    }
    
  };

  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