#include "laguerre2bessel.h"
#include "lagsht_bessel.h"
#include "laguerreTransform.h"
#include "lagsht_exceptions.h"



#include "walltimer.h"    //timing
#include <fstream>       //dump
#include <sstream>
#include <numeric>


#define DEBUG 0

namespace LagSHT {

  /* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre 
     mais ca n'a pas l'air de converger (28 Sept 15)
     Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
     recursive
     Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a  Clenshaw-Curtis quadrature and global integration strategy. 
     Version 3: calcul via une quadrature de  Clenshaw-Curtis quadrature and global integration en controlant les bornes
                d'integration
  */
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
				 LaguerreTransform* lagTrans,
				 int Nmax, int Pmax, r_8 Rmax,
				 string clenshawDir,
				 string jlnpDir,
				 bool recomputeJlnp
				 ):
  sphere_(sphere),lagTrans_(lagTrans),
  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
  //JEC 1/2/16  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmaax), PPmax_(Pmax) {

  //  PPmax_ = 4096;

  
  if(sphere_) {
    //store localy usefull parameters for the algorithms Analysis/Synthesis
    Npix_ = sphere_->Npix();
    Nalm_ = sphere_->Nalm();
    Lmax_ = sphere_->Lmax();
    Mmax_ = sphere_->Mmax();
  } else {
    throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
  }
  
  //  jnu_ = new Bessel(Lmax_,std::max(std::max(Pmax_,PPmax_),Nmax_));  //the max is here just in case PPmax =/= Pmax
  jnu_ = new Bessel(Lmax_,std::max(Pmax_,Nmax_));  //the max is here just in case Nmax =/= Pmax
  
  if(!recomputeJlnp){
    
    //read the Rmax-independant values of the jlnp and rescale
    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
    
    string fname(jlnpDir+"/"+ss.str());
    std::ifstream ifs(fname.c_str(), std::ifstream::in);
    if(!ifs.is_open())
      throw LagSHTError("ERROR: jlnp file open failed.:"+ss.str());
    
    r_8 tau = lagTrans_->GetTau();
    r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
    int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
    Jlnp_.resize(Nmax_*Jstride);
    
    int pcur, lcur, ncur;
    r_8 val;
    for(int p=0; p<Pmax_; p++){
      for(int l=0; l<Lmax_; l++){
	int ploff7 = p + l*Pmax_;
	for (int n=0; n<Nmax_; n++){
	  ifs >> lcur >> ncur >> pcur >> val;         //read
// 	  if( (l!=lcur) || (n!=ncur) || (p!=pcur))                         //check
// 	    throw LagSHTError("ERROR: jlnp file read failed.:"+ss.str());
	  Jlnp_[ploff7 + n*Jstride] =  tau3sqrt * val;                     //rescaling
	}//n
      }//l
    }//p
    
  } else {
    ComputeJInteg(clenshawDir,jlnpDir);
  }//jlnp omputation/load
  
}//Ctor


/*!
  Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
  with
  K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
  k_{lp}*tau = q_{lp}/r_{Nmax-1}
  Apply a Discrete Spherical transform Lanusse et al A&A 540, A92 (2012) Eq.40 with l'=l
  Jlnp ~ Kapa_l^(-3) tau^(3/2) Sqrt(2 Pi) Sum_{p'=0}^{PPmax-1} Kn(r_{lp'}) jl(q_{lp'}q_{lp}/q_{l,PPmax-1})/[jl+1(q_{lp})]^2
 with r_{lp'} = q_{lp'}/q_{l,PPmax-1}r_{Nmax-1}
 notice PPmax as nothing to do in principle with Pmax the Fourier-Bessel max value of the 3rd index
 Kapa_l = q_{l,PPmax-1}/r_{Nmax-1}



 JEC: 12/11/15 Add the contribution of Integral_{r_{Nmax-1}}^{r_{Nmax-1} + 3(r_{Nmax-1}-r_{Nmax-2})} by  Clenshaw-Curtis quadrature

void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){

  r_8 tau = lagTrans_->GetTau();
  r_8 tau3sqrt = pow(tau,(r_8)(3./2.));

//   cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
  
  //JEC 12/11/15 Start
  r_8 rNm1 = Rmax_ / tau; // R = r[N-1] voir acces direct  = (lagTrans_->GetNodes()).back()
  r_8 rNm2 = (lagTrans_->GetNodes())[Nmax_-2];
  r_8 deltaR = (rNm1-rNm2);
  r_8 integLowerBound = rNm1;
  r_8 integUpperBound = rNm1+3.0*deltaR; //for the integral 
  r_8 tol = 1e-10;
  typedef  ClenshawCurtisQuadrature<double> Integrator_t;
  string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
  Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
  
  Quadrature<r_8,Integrator_t>::values_t integ_val;
  //JEC 12/11/15 end
  
  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
  Jlnp_.resize(Nmax_*Jstride);  

  //JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
  int alpha = lagTrans_->GetAlpha();
  r_8 alphaFact = 1; // alpha!
  for(int i=1;i<=alpha; i++) alphaFact *= i;
  alphaFact  = sqrt((r_8)alphaFact); //sqrt(alpha!)
  r_8 invalphaFact = 1.0/alphaFact; // 1/sqrt(alpha)!                   
  
  vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) en iteratif
  facts[0] = invalphaFact;
  for(int n=1;n<Nmax_;n++) {
    facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );
  }
  //END
  
    
  int PPmaxm1 = PPmax_ -1;

  tstack_push("J Integ compute ");

  stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
  std::ofstream ofs;
  string fname(jlnpDir+"/"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);


//   cout << "(JEC) Jlpn recompute..." << endl;
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){
      
      int ploff7 = p + l*Pmax_;
      r_8 qlp = (*jnu_)(l,p);                //ql,p


      r_8 klp = qlp/rNm1; //JEC 25/1/16
      
      r_8 qlpx = (*jnu_)(l,PPmaxm1);         //ql,PPmax-1
      
      r_8 Kapa_l = qlpx/rNm1;                  // Kapa_l = (ql,PPmax-1)/R
      r_8 Kapa_l3 =  Kapa_l *  Kapa_l *  Kapa_l;
      
      r_8 coeff = sqrt(2.0 *M_PI)/Kapa_l3;  // Sqrt(2 Pi)/Kapa_l^3
      
      for (int n=0; n<Nmax_; n++){
	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);


	
	
	r_8 sum = 0.;
	//contribution of the integral: [0,rNm1] estimated
	//using Discrete Fourier-Bessel Transform

	for (int pp=0; pp<PPmax_; pp++){
	  r_8 qlpp = (*jnu_)(l,pp);           //ql,pp
	  r_8 jlp1 = Bessel::jn(l+1,qlpp);       //j(l+1,ql,pp)
	  
	  r_8 rlpp = qlpp/qlpx * rNm1;          // rl,pp = (ql,pp/ql,PPmax-1) * rNm1
	  r_8 jlqq = Bessel::jn(l,qlpp * qlp/qlpx);
	  
	  sum +=  jlqq/(jlp1*jlp1) * Ln->Value(rlpp);

	}//pp-loop
	

	//	Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
	//JEC 23/11/15 save the Rmax-independant part
	Jlnp_[ploff7 + n*Jstride] = coeff * facts[n] * sum;
	
	//contribution [rNm1, rNm1+deltaR] estimated with Clenshaw-Curtis quadrature (JEC 12/11/15) Start
	IntFunc1D f(l,qlp/rNm1,Ln);
	theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
	integ_val = theQuad.GlobalAdapStrat(tol,2,1);
	//	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
	//JEC 23/11/15 save the Rmax-independant part
	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
	

	//Limber approx JEC 25/1/16 the sqrt(2/Pi) is canceled by the definition of SphericalBesselJ function used in this approx
	r_8 limber = pow(klp,(r_8)(-3.0)) * pow((r_8)(l+0.5),(r_8)1.5) * Ln->Value(((r_8)(l)+0.5)/klp)*facts[n];

	//New Approx JEC 27/1/16
	theQuad.SetFuncBounded(f,0.8*r_8(l),rNm1+4.0*deltaR);
	integ_val = theQuad.GlobalAdapStrat(tol,2,1);
	

	cout <<l<<" "
	     <<n<<" "
	     <<p<< " "
	     << setprecision(30) << Jlnp_[ploff7 + n*Jstride]  
	     << " " << limber
	     << " " << integ_val.first *  (facts[n]*sqrt(2.0/M_PI))
	     << endl;
		
		
	//JEC 23/11/15 save the Rmax-independant part
	ofs <<l<<" "
	    <<n<<" "
	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
		
	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;

	delete Ln;
	
      }//n
    }//l
  }//p
  
  ofs.close();

  tstack_pop("J Integ compute ");
  
}//ComputeJlpnInteg
*/

/*!
  Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
  with
  K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
  k_{lp}*tau = q_{lp}/r_{Nmax-1}
  
  Use a Clenshaw-Curtis quadrature with 2*40-1 pts between LowBound and UppBound
  - LowBound is (80% * l_value)
  - UppBound is the Max between : 
      - for n>=2 nodes[n-1] + 5*(nodes[n-1]-nodes[n-2]) this is an estimate of Kn cut-off
      - q(l,1)/(k_{lp}*tau) this is an estimate of first oscillation of Bessel function
*/
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){


  if(Nmax_<2)throw LagSHTError("Laguerre2Bessel::ComputeJInteg Nmax<2");

  tstack_push("J Integ compute ");


  r_8 tau = lagTrans_->GetTau();
  r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
  r_8 rNm1 = Rmax_ / tau; // R = r[N-1] with N the original transform value. Todo (lagTrans_->GetNodes()).back()
  r_8 tol = 1e-10;
  
  typedef  ClenshawCurtisQuadrature<double> Integrator_t;
  string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-40.txt";
  Integrator_t theQuad(40,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
  
  Quadrature<r_8,Integrator_t>::values_t integ_val;
  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
  Jlnp_.resize(Nmax_*Jstride);  

  r_8 alpha = lagTrans_->GetAlpha();
  vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) the normalization
  facts[0] = 1.0/sqrt(boost::math::tgamma(alpha+1.0)); //1/sqrt(alpha)!       
  for(int n=1;n<Nmax_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );

 
  vector<r_8> integUpperBound(Nmax_);
  integUpperBound[0] = rNm1; //this is the case of L_0^{\alpha}(x) =1 with no root
  integUpperBound[1] = rNm1; //this is the case of L_1^{\alpha}(x) =1+alpha-x with single root = 1/(1+alpha)
  for (int n=2;n<Nmax_;n++){
    LaguerreFuncQuad lag(n,alpha);
    vector<r_8> nodes;
    vector<r_8> weights;
    lag.QuadWeightNodes(nodes,weights);
    integUpperBound[n] = nodes[n-1] + 5*(nodes[n-1]-nodes[n-2]); //upper bound from Laguerre function end point, may be modified by Bessel part later
  }


  stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
  std::ofstream ofs;
  string fname(jlnpDir+"/"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);
  
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){

      int ploff7 = p + l*Pmax_;
      r_8 qlp = (*jnu_)(l,p);                //ql,p
      r_8 klptau = qlp/rNm1;
      r_8 intLowBound = 0.8*l;

      for (int n=0; n<Nmax_; n++){


	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
	IntFunc1D f(l,klptau,Ln);


	//upper bound of integration is the Max between the 2nd BesselRoot rescaled and the Laguerre Function "end point" 
	r_8 intUppBound = std::max(integUpperBound[n], (*jnu_)(l,1)/klptau);
	

	theQuad.SetFuncBounded(f,intLowBound,intUppBound);
	integ_val = theQuad.GlobalAdapStrat(tol);
	Jlnp_[ploff7 + n*Jstride] = integ_val.first * (facts[n]*sqrt(2.0/M_PI));

 	cout <<l<<" "
	     <<n<<" "
	     <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
	ofs <<l<<" "
	    <<n<<" "
 	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
	
	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;

	delete Ln;

      }//n
    }//l
  }//p


  ofs.close();

  tstack_pop("J Integ compute ");


}//ComputeJlpnInteg


/*!
  FBlmp = FB_{lm}(k_{lp}) = Sum_n FL_{lmn} J_{ln}(k_{lp}) = Sum_n FL_{lmn} J_{lnp}
  the J_{lnp} includes the sqrt(2/pi) contrary to McEwen paper Exact Wavelet on the Ball
  This is the Fourrier-Bessel coefficients related to McEwen & Leistedh definition
  In case one wants to use Lanusse & Starck definition then
  FBlmp_Lanusse = k_{lp} * FBlmp_McEwen
*/
void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
				    vector< complex<r_8> >& FBlmp){

  
  int Jstride = Lmax_*Pmax_;

#if DEBUG>=1
  std::ofstream ofs;
  ofs.open ("FBlmp.txt", std::ofstream::out);
#endif

  //try here a pedestrian computation
  for(int p=0; p<Pmax_; p++){
    for(int l=0; l<Lmax_; l++){

      int ploff7 = p + l*Pmax_;

      for (int m=0; m<=l; m++) {
	int almoff7 = l+m*Lmax_-m*(m+1)/2;
	int idBess = p*Nalm_ + almoff7;
	complex<r_8> sum = 0.0;

	for (int n=0; n<Nmax_; n++){
	  int idLag =  n*Nalm_ + almoff7;
	  sum += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
	}//n
	FBlmp[idBess] = sum;

#if DEBUG>=1
       	ofs <<l<<" "
	    <<m<<" "
	    <<p<< " " << setprecision(12) << 	FBlmp[idBess] << endl;
#endif

      }//m
    }//l
  }//p

#if DEBUG>=1
  ofs.close();
#endif
  

}//Transform


void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp, 
				vector< complex<r_8> >& FBalmk,
				vector<r_8>& fijk
				) {

  int NtotCoeffFB = Nalm_*Pmax_; 
  if(FBlmp.size() != (size_t)NtotCoeffFB) throw LagSHTError("Laguerre2Bessel::Synthesis size error");

  cout << "Dump  Laguerre2Bessel::Synthesis START "<<endl;
    

  int NtotCoeffFL = Nalm_*Nmax_; 
  FBalmk.resize(NtotCoeffFL); 
  
  /*
    This is the McEwen & Leistedh definition of the Fourier-Bessel Transform
    FBalmk = alm^{FB}(R_k) = Sum_{p=0}^{Pmax-1} Kapa_{lp} FBlmp  j_l(k_{lp} R_k)
    with
    kapa_{lp} = Sqrt(2pi) Rmax^{-3}/j_{l+1}^2(q_{lp})
    k_{lp} = q_{lp}/Rmax
    R_k = tau r_k
    r_k: k-th nodes of Laguerre of order Nmax as for the Spherical-Laguerre Transform.
    tau: Rmax/r_{N-1} from Laguerre-Transform
    => k_{lp} R_k = q_{lp} r_k/r_{N-1}

    Rdeminder: Fourier-Laguerre case where
    alm^{FL}(R_k) = Sum_{n=0}^{Nmax-1} FLlmn LaguerreFunc_n(r_k)
    
    in the Fourier-Bessel case the element of the sum, except the FBlmp, are depending
    on "l" also throw kapa_{lp}, k_{lp} and j_l(k_{lp}.) 
   
   */
  
  vector<r_8> vrk = lagTrans_->GetNodes();
  r_8 rNm1 = vrk[Nmax_-1];

  //pedestrian computation
  r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);

  vector< complex<r_8> > sumElmt(Pmax_);
  for(int k=0;k<Nmax_;k++){

    r_8 rk_red = vrk[k]/rNm1;

    for(int l=0;l<Lmax_;l++){
      for(int m=0;m<=l;m++){ 

	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
	int id0 = k*Nalm_ + lmOff7;

	for(int p=0;p<Pmax_;p++){ //Fill
 	  r_8 qlp = (*jnu_)(l,p);
 	  r_8 jlp1 = Bessel::jn(l+1,qlp);	  
 	  r_8 kapalp = fact/(jlp1*jlp1);
 	  int id1 = p*Nalm_ + lmOff7;

  	  sumElmt[p] = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
	}
	//kahan summation	FBalmk[id0] = LagSHT::accumulate(sumElmt.begin(),sumElmt.end());
	FBalmk[id0] = std::accumulate(sumElmt.begin(),sumElmt.end(),complex<r_8>(0.,0.));

      }//loop on m
    }//loop on l
  }//loop on k


  //from Almk to Fijk with libsharp

  int Ntot3DPix = Npix_*Nmax_;
  fijk.resize(Ntot3DPix); 

   int off7n, off7k;
   //perform the synthesis per shell
   for(int n =0; n<Nmax_; n++){
     off7n = n*Nalm_;
     off7k = n*Npix_; 

     sphere_->alm2map(reinterpret_cast<const r_8*>(&(FBalmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation

   }//loop on shell
  

}//Synthesis


void Laguerre2Bessel::Print(ostream& os, int what){
  if(what>0){

    //Debug the Jlnp coefficients
    int stride = Lmax_*Pmax_;
    for (int pc = 0; pc < Pmax_; pc++) {
      for (int lc = 0; lc < Lmax_; lc++) {
	int ploff7 = pc + lc*Pmax_;
	for (int nc = 0; nc < Nmax_; nc++) {
	  os << lc<<" "<<nc<<" "<<pc<<" " << setprecision(12)<< Jlnp_[ploff7 + stride*nc] << endl;
	}//loop nc
      }//loop pc
    }//loop lc
  }//what
}//Print

}//namespace