(JEC) 3/11/15 new version of laguerre2bessel. The Jlnp intgral is computed...

(JEC) 3/11/15 new version of laguerre2bessel. The Jlnp intgral is computed with a Discret Fourier-Bessel Transformation.

Makefile

@@ -121,8 +121,8 @@ CXXHDR = $(SRC)lagsht_exceptions.h \
 	$(SRC)lagSphericTransform.h \
 	$(SRC)walltimer.h \
 	$(SRC)lagsht_bessel.h \
-	$(SRC)laguerre2bessel.h \
-	$(SRC)quadinteg.h
+	$(SRC)laguerre2bessel.h
+#	$(SRC)quadinteg.h
 
 CPPFLAGS  += $(BLASYES) $(SHARPINC) $(BLASINC)
 LDFLAGS   += $(SHARPLIBN) $(BLASLIBN)  -lm

src/lagsht_bessel.cc

@@ -7,6 +7,7 @@ namespace LagSHT {
 void Bessel::BesselRoots() {
   //  using namespace std;
 
+
   string msg;
   try {
     unsigned int n_zeros = nmax_;
@@ -16,6 +17,7 @@ void Bessel::BesselRoots() {
       r_8 order = (r_8)lcur + 0.5;
       vector<r_8> zeros;
       boost::math::cyl_bessel_j_zero(order,1,n_zeros, back_inserter(zeros));
+
       copy(zeros.begin(),zeros.end(),qln_.begin() + lcur*nmax_);
       
     }//l loop   

src/lagsht_bessel.h

@@ -26,6 +26,7 @@
 
 
 //STD
+#include <iostream>
 #include <vector>
 #include <algorithm>
 //Boost
@@ -69,10 +70,15 @@ class  Bessel {
   }//end Init
 
   // Sperical Bessel functions
-  static inline r_8 jn(int n, r_8 x) { return boost::math::sph_bessel(n, x); }
+  static inline r_8 jn(int n, r_8 x) { 
+    return boost::math::sph_bessel(n, x); 
+  }
 
   //Helper
-  r_8 operator()(int l, int n) const { return qln_[l*nmax_+n];} //qln(l,n)
+  r_8 operator()(int l, int n) const { 
+    if(l>lmax_ || n > nmax_) cout <<" ERROR Besel::Op(l,n) called with ("<<l<<","<<n<<") while "
+      << " lmax = " << lmax_ << " nmax = " << nmax_  <<endl;
+    return qln_[l*nmax_+n];} //qln(l,n)
   void GetVecRoots(std::vector<r_8>& roots, int l) const {
     roots.resize(nmax_);
     copy(qln_.begin()+l*nmax_,qln_.begin()+(l+1)*nmax_,roots.begin());

src/lagsht_testsuite.cc

@@ -383,51 +383,110 @@ void TestLaguerre2Bessel() {
   tstack_pop("Compute Alm(rk) from FB coeffs.");
   
   cout << "tstack 2 End" <<endl;
+
+  //verif F-Laguerre part
+  tstack_push("F-L Analysis verif");
+  vector< complex<r_8> > flmn(Ntot);
+  sphlagtrans.Analysis(fFLijk,flmn);
+  {
+    r_8 err_abs(0.);
+    r_8 err_rel(0.);
+    int imax = -1;
+    for(int i=0;i<Ntot;i++){
+      
+      complex<r_8> cdiff = (flmnOrig[i] -  flmn[i]) * conj(flmnOrig[i] -  flmn[i]);
+      r_16 diff = sqrt(cdiff.real());
+      if(diff>err_abs){ 
+	err_abs = diff;
+	imax = i;
+      }
+      complex<r_8> foriCAbs = flmnOrig[i]*conj(flmnOrig[i]);
+      r_8 foriAbs = sqrt(foriCAbs.real());
+      r_8 relatdiff = diff/foriAbs;
+	
+      if((relatdiff)>err_rel) err_rel = relatdiff;
+    }
+    cout << " >>>>>>>>>>>>>>>>>>>>> Fourrier-Laguerre part <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+    cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
+    cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+    
+  }
+  tstack_pop("F-L Analysis verif");
+
+
   tstack_pop("processing");
   
   
-  {//dump 1
+  {//Check 1
     cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
 
     std::ofstream ofs;
     ofs.open ("Almi.txt", std::ofstream::out);
     
     for(int n=0;n<Nmax;n++){
+      r_8 err_abs(0.);
+      r_8 err_rel(0.);
+
       for(int l=0;l<Lmax;l++){
 	for (int m=0;m<=l;m++) {
 	  int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
-	  cout << "Almk ("
-	       <<l<<","
-	       <<m<<","
-	       <<n<< "): FL : " << setprecision(12) << FLalmk[id] 
-	       << " FB : " << FBalmk[id]  << endl;
 	  ofs <<l<<" "
 	      <<m<<" "
 	      <<n<< " " << setprecision(12) << FLalmk[id].real() << " " <<   FLalmk[id].imag()
 	      << "  " << FBalmk[id].real() << " " << FBalmk[id].imag()  << endl;
 	  
+	  complex<r_8> cdiff = (FLalmk[id] -  FBalmk[id])*conj(FLalmk[id] -  FBalmk[id]);
+	  r_8 diff = sqrt(cdiff.real());
+	  if(diff>err_abs){
+	    err_abs = diff;
+	  }
+	  complex<r_8> cref= FLalmk[id]*conj(FLalmk[id]);
+	  r_8 ref = sqrt(cref.real());
+	  r_8 relatdiff = diff/ref;
+	  if(relatdiff>err_rel) err_rel = relatdiff;
 	}
       }
+      
+      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<n<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+      cout << "Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
+      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+
     }
 
     ofs.close();
-  }//dump 1
+  }//check 1
 
-  {//dump 2
+  {//check 2
     std::ofstream ofs;
     ofs.open ("fijk.txt", std::ofstream::out);
 
+    cout << "Dump FL or FB reconstructed fijk on each shell k" <<endl;
+
     for (int ish=0;ish<Nshell; ish++){
+      r_8 err_abs(0.);
+      r_8 err_rel(0.);
       for (int ip=0; ip<Npix; ip++) {
 	int id = ish*Npix+ip;
-	cout << "Shell("<<ish<<") Pix["<<ip<<"] FL fijk = " << fFLijk[id] << " FB fijk = " << fFBijk[id] << endl;
+	//	cout << "Shell("<<ish<<") Pix["<<ip<<"] FL fijk = " << fFLijk[id] << " FB fijk = " << fFBijk[id] << endl;
 	ofs << ish << " " << ip << " " << setprecision(12) <<  fFLijk[id] << " " << fFBijk[id] << endl;
-      }
-    }
+	r_8 diff = fabs(fFLijk[id]-fFBijk[id]);
+	if(diff>err_abs){
+	  err_abs = diff;
+	}
+	r_8 relatdiff = diff/ fFLijk[id];
+	if(relatdiff>err_rel) err_rel = relatdiff;
+      }//loop on px
+      
+            
+      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<ish<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+      cout << "Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
+      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+
+    }//loop on shell
 
 
     ofs.close();
-  }//dump 2
+  }//check 2
 
   
 
@@ -615,7 +674,7 @@ int main(int narg, char *arg[]) {
 
     case 5:
       {
-	if(Pmax<N){ param.Pmax = 16*N; cout << "(Warning): Pmax<Nmax => modify the value to "<< Pmax << endl;}
+	if(Pmax<N){ param.Pmax = 16*N; cout << "(Warning): Pmax<Nmax => modify the value to "<<  param.Pmax << endl;}
 	tstack_push("TestLaguerre2Bessel");
 	TestLaguerre2Bessel();
  	tstack_pop("TestLaguerre2Bessel");

src/laguerre2bessel.cc

@@ -7,8 +7,6 @@
 
 #include "walltimer.h"    //timing
 
-#include <omp.h> //TEST 
-
 namespace LagSHT {
 
   /* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre 
@@ -37,67 +35,92 @@ namespace LagSHT {
 
 }//Ctor
 
-  /*! compute J_{lnp} = J_{ln}(k_{lp} via Clenshaw-Curtis quadrature and local integration strategy
-   */
-  void Laguerre2Bessel::ComputeJInteg(){
-
-    typedef  ClenshawCurtisQuadrature<double> Integrator_t;
-    Integrator_t theQuad(20,"./data/ClenshawCurtisRuleData-20.txt",false); //false=do not recompute the weights and nodes of the quadrature
-
-    r_8 tau = lagTrans_->GetTau();
-    r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
-    //    cout << "ComputeJlpnInteg: tau " << tau << ", tau^(3/2): " << tau3sqrt << endl;
-
-
-    //normalisation factor 
-    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(2/pi)*tau^(3/2)*sqrt(n!/(n+alpha)!) en iteratif
-    facts[0] = invalphaFact*tau3sqrt*sqrt(2.0/M_PI); //JEC 7/10/15 add sqrt(2/pi)
-    for(int n=1;n<Nmax_;n++) {
-      facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) ); 
-    }
-
-
-    jnu_ = new Bessel(Lmax_,Pmax_);
-    int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
-    Jlnp_.resize(Nmax_*Jstride);
 
+/*!
+  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}
+ */
+void Laguerre2Bessel::ComputeJInteg(){
+
+  r_8 tau = lagTrans_->GetTau();
+  r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
+  r_8 R = Rmax_ / tau; // R = r[N-1] voir acces direct  = (lagTrans_->GetNodes()).back()
+  
+  int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
+  Jlnp_.resize(Nmax_*Jstride);
+  
+  
+  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(2/pi)*tau^(3/2)*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)) );
+  }
+  
     
-    r_8 integLowerBound = 0.;
-    r_8 integUpperBound = 10*Rmax_/tau; //for the integral 
-    r_8 tol = 1e-6;
-
-    tstack_push("J Integ compute");
-
-    Quadrature<r_8,Integrator_t>::values_t loc_val;
-    for(int p=0; p<Pmax_; p++){
-      for(int l=0; l<Lmax_; l++){
-
-	int ploff7 = p + l*Pmax_;
- 	r_8 klp = (*jnu_)(l,p)/Rmax_*tau; //tau-rescaling of the original klp definition
-
-	for (int n=0; n<Nmax_; n++){
-
-	  LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
-
-	  IntFunc1D f(l,klp,Ln);
-	  theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
-	  loc_val = theQuad.LocalAdapStrat(tol);
-	  Jlnp_[ploff7 + n*Jstride] = loc_val.first*facts[n];
-
-	  delete Ln;
-	}//n
-      }//l
-    }//p
+  int PPmax = max(128,Pmax_); // Sum maximum    (DOIT ETRE > PMax pour la Synthese car appel a qlp
+  int PPmaxm1 = PPmax -1;
+  
+  jnu_ = new Bessel(Lmax_,PPmax);
+  
+ 
+  tstack_push("J Integ compute ");
+   
+  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 qlpx = (*jnu_)(l,PPmaxm1);         //ql,PPmax-1
+      
+      r_8 Kapa_l = qlpx/R;                  // 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;  // tau^(3/2) Sqrt(2 Pi)/Kapa_l^3
+      
+      for (int n=0; n<Nmax_; n++){
+	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
+	
+	r_8 sum = 0.;
+
+	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 * R;          // rl,pp = (ql,pp/ql,PPmax-1) * R 
+	  r_8 jlqq = Bessel::jn(l,qlpp * qlp/qlpx);
+	  
+	  sum += jlqq/(jlp1*jlp1) * Ln->Value(rlpp);
+	}
+	
+
+	Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
+		
+	delete Ln;
+	
+      }//n
+    }//l
+  }//p
+  
 
-    tstack_pop("J Integ compute");
+  tstack_pop("J Integ compute ");
+  cout << "J Integ compute  END" << endl;
+  
+}//ComputeJlpnInteg
 
-  }//ComputeJlpnInteg
 
 /*!
   FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{ln}(k_{lp})
@@ -158,10 +181,7 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
     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
-
-    23/9/15 in fact Leistedt & McEwen  et al in IEEE Vol 60  Num 12 (Dec 2012) do not 
-    use the same convention on Bessel Transform as Lanusse et al. A&A 578, A10 (2015) so
-    in the Sum_{p=0}^{Pmax-1}... the k_{lp} should be squared
+    => 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)
@@ -172,7 +192,6 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
    */
   
   vector<r_8> vrk = lagTrans_->GetNodes();
-  //OK  cout << "Verif: tau " << lagTrans_->GetTau() << ", Rmax/r[N-1]: " << Rmax_/vrk[Nmax_-1] << endl;
   r_8 rNm1 = vrk[Nmax_-1];
 
   //pedestrian computation
@@ -183,7 +202,7 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
     r_8 rk_red = vrk[k]/rNm1;
 
     for(int l=0;l<Lmax_;l++){
-      for(int m=0;m<=l;m++){ //TODO use Mmax also to limit "m"
+      for(int m=0;m<=l;m++){ 
 
 	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
 	int id0 = k*Nalm_ + lmOff7;
@@ -193,7 +212,9 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 	for(int p=0;p<Pmax_;p++){
 
 	  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;
 
@@ -205,8 +226,10 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
   }//loop on k
 
 
+  //from Almk to Fijk with libsharp
+
   int Ntot3DPix = Npix_*Nmax_;
-  fijk.resize(Ntot3DPix); //Not used for the moment
+  fijk.resize(Ntot3DPix); 
 
   int off7n, off7k;
   //perform the synthesis per shell

src/laguerre2bessel.h

@@ -33,8 +33,6 @@
 #include "lagsht_spheregeom.h"
 #include "lagsht_bessel.h"
 
-#include "quadinteg.h"
-
 using namespace std;
 
 namespace LagSHT {
@@ -88,19 +86,6 @@ private:
   */
   void ComputeJInteg();
 
-  class IntFunc1D : public ClassFunc1D {
-  public:
-    IntFunc1D(int l, r_8 klp,LaguerreFuncQuad* Ln):
-      l_(l),klp_(klp),Ln_(Ln) {}
-    virtual double operator()(double x) const {
-      return (x*x)*(Bessel::jn(l_,x*klp_))*Ln_->Value(x);
-    }
-    virtual ~IntFunc1D() {}
-  private:
-    int l_;
-    r_8 klp_;
-    LaguerreFuncQuad* Ln_; //no ownership
-  };//class IntFunc1D
  
 protected:
   BaseGeometry* sphere_; //<! the 2D pixelization of the shells (not the owner)
@@ -118,8 +103,6 @@ protected:
   r_8 Rmax_; //!< Radial maximal value  
 
   
-  //JEC 28/9/15  int K0_; //!< Order of the LaguerreTransform to compute the Bessel-Laguerre scalar products Jlpn
-  
   Bessel* jnu_; //!< Spherical Bessel fnt (owner of the ptr)
   vector<r_8> Jlnp_; //!< the  Bessel-Laguerre scalar product
   

src/laguerreTransform.cc

@@ -208,6 +208,7 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
 #if CBLAS
     cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fi[0]), 2*stride);
 #else
+    cout << "CBLAS desactivated................" << endl;
     vector<complex<r_8> > vtmp(N_);
     for (int l=0; l<stride; l++) {
       vtmp.assign(N_,0.);