(JEC) 12/11/15 jlnp via DFBT before re-introduction of same numerical integration

src/laguerre2bessel.cc

@@ -6,6 +6,7 @@
 
 
 #include "walltimer.h"    //timing
+#include <fstream>       //dump
 
 namespace LagSHT {
 
@@ -51,6 +52,9 @@ void Laguerre2Bessel::ComputeJInteg(){
 
   r_8 tau = lagTrans_->GetTau();
   r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
+  
+  cout << "JEC: " << tau3sqrt*sqrt(2.0/M_PI) << endl;
+
   r_8 R = Rmax_ / tau; // R = r[N-1] voir acces direct  = (lagTrans_->GetNodes()).back()
   
   int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
@@ -70,13 +74,16 @@ void Laguerre2Bessel::ComputeJInteg(){
   }
   
     
-  int PPmax = max(128,Pmax_); // Sum maximum    (DOIT ETRE > PMax pour la Synthese car appel a qlp
+  int PPmax = max(512,Pmax_); // Sum maximum    (DOIT ETRE > PMax pour la Synthese car appel a qlp
   int PPmaxm1 = PPmax -1;
   
-  jnu_ = new Bessel(Lmax_,PPmax);
+  jnu_ = new Bessel(Lmax_,PPmax);  //To do put it in constructor
   
- 
   tstack_push("J Integ compute ");
+
+  std::ofstream ofs;
+  ofs.open ("jlnp.txt", std::ofstream::out);
+
    
   for(int p=0; p<Pmax_; p++){
     for(int l=0; l<Lmax_; l++){
@@ -89,7 +96,7 @@ void Laguerre2Bessel::ComputeJInteg(){
       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
+      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);
@@ -104,17 +111,22 @@ void Laguerre2Bessel::ComputeJInteg(){
 	  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]);
 		
+	ofs <<l<<" "
+	    <<n<<" "
+	    <<p<< " " << setprecision(12) << Jlnp_[ploff7 + n*Jstride] << endl;
+	
 	delete Ln;
 	
       }//n
     }//l
   }//p
   
+  ofs.close();
 
   tstack_pop("J Integ compute ");
   cout << "J Integ compute  END" << endl;
@@ -144,12 +156,13 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
       for (int m=0; m<=l; m++) {
 	int almoff7 = l+m*Lmax_-m*(m+1)/2;
 	int idBess = p*Nalm_ + almoff7;
-	FBlmp[idBess] = 0.0;
+	complex<r_8> sum = 0.0;
 
 	for (int n=0; n<Nmax_; n++){
 	  int idLag =  n*Nalm_ + almoff7;
-	  FBlmp[idBess] += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
+	  sum += FLlmn[idLag]*Jlnp_[ploff7 + n*Jstride];	  
 	}//n
+	FBlmp[idBess] = sum;
 
       }//m
     }//l
@@ -207,7 +220,8 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
 	int id0 = k*Nalm_ + lmOff7;
 
- 	FBalmk[id0]=0.;
+	
+	complex<r_8> sum = 0.;
 
 	for(int p=0;p<Pmax_;p++){
 
@@ -218,9 +232,12 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 	  r_8 kapalp = fact/(jlp1*jlp1);
 	  int id1 = p*Nalm_ + lmOff7;
 
-	  FBalmk[id0] += kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
+	  sum += (kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1]);
 
 	}//loop on p
+	FBalmk[id0] += sum;
+
+
       }//loop on m
     }//loop on l
   }//loop on k

test/Makefile

@@ -37,10 +37,10 @@ SHARPLIBN = -L$(SHARPLIB) -lsharp -lc_utils -lfftpack
 
 
 #======= LagSHT ======
-LAGSHTDIR = ..
+LAGSHTDIR = /Users/campagne/Travail/Software/Laguerre/LagSHT
 LAGSHTLIB = $(LAGSHTDIR)/lib
 LAGSHTINC = -I$(LAGSHTDIR)/include
-LAGSHTLIBN = $(LAGSHTLIB)/liblagsht.a
+LAGSHTLIBN = -L$(LAGSHTLIB) -llagsht
 
 
 SRC = ./src/
@@ -95,7 +95,7 @@ CXXSHOBJ =
 CXXHDR = 
 
 CPPFLAGS  += $(BLASYES) $(SHARPINC) $(BLASINC) $(LAGSHTINC)
-LDFLAGS   += $(SHARPLIBN) $(BLASLIBN)  $(LAGSHTLIBN)  -lm
+LDFLAGS   += $(LAGSHTLIBN) $(SHARPLIBN) $(BLASLIBN)    -lm
 
 
 ######################

test/test3D-t3.txt

-Max Memory size: 8589 MBytes
-Configuration parameters are set to:  Test number : 3
- Lmax, Nmax, Pmax, Rmax , alpha: 1024, 128, 0, 5, 2
-Geometry: Gauss Ntheta, Nphi: -1, -1
-Generate Cl(k) spectra and fill a 3D-ball
-Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell
-Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn
-Error analysis on Cl(k) input and output of the Laguerre Spherical Transf.
- >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<
- Err. Max. 4.65397e-08 [0, 298], Err. Rel. 3.97506e-05
- >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<
-
-Total wall clock time for 'Test 3': 51.8014s
-|
-+- Generate fijk 3D-ball : 54.67% (28.3221s)
-+- Lag Spherical Transf. : 42.01% (21.7609s)
-|  |
-|  +- SHT Analysis          : 90.37% (19.6658s)
-|  +- Laguerre MultiAnalysis:  4.23% ( 0.9205s)
-|  +- <unaccounted>         :  5.40% ( 1.1746s)
-|  
-+- Clk and Cln OUT       :  2.62% ( 1.3590s)
-+- error analysis        :  0.00% ( 0.0006s)
-+- <unaccounted>         :  0.69% ( 0.3588s)
-
-
-Accumulated timing overhead: approx. 0.0000s
----/ Fin bloc try ---- 
- ---- Programme test3D.cc -  FIN  (Rc=0) --- 

test/test3D-t5.txt

-Max Memory size: 8589 MBytes
-Configuration parameters are set to:  Test number : 5
- Lmax, Nmax, Pmax, Rmax , alpha: 1024, 128, 0, 5, 2
-Geometry: Gauss Ntheta, Nphi: -1, -1
-Generate a 3D-Cube 1024^3
-Fill the 3D-Ball from the Cube
- Fill a 3D-Ball (2,2,2)Mpc^3: Rmax is set to : 819.2
-Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell
-Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn
-
-Total wall clock time for 'Test 5': 274.0201s
-|
-+- Generate a 3D-Cube             : 42.19% (115.6206s)
-+- Fill the  3D-Ball from the Cube: 41.79% (114.4995s)
-+- Lag Spherical Transf.          : 14.48% ( 39.6700s)
-|  |
-|  +- SHT Analysis                   : 79.37% ( 31.4845s)
-|  +- Laguerre MultiAnalysis         :  9.26% (  3.6734s)
-|  +- <unaccounted>                  : 11.37% (  4.5120s)
-|  
-+- Clk and Cln OUT                :  0.44% (  1.2128s)
-+- <unaccounted>                  :  1.10% (  3.0172s)
-
-
-Accumulated timing overhead: approx. 0.0010s
----/ Fin bloc try ---- 
- ---- Programme test3D.cc -  FIN  (Rc=0) --- 

test/test3D.cc

@@ -8,7 +8,7 @@
 #include "LagSHT/walltimer.h"    
 #include "LagSHT/lagSphericTransform.h"
 
-
+#include "LagSHT/laguerre2bessel.h" 
 
 using namespace LagSHT;
 
@@ -173,6 +173,7 @@ void test3() {
   //initial values
   int Lmax = param.Lmax;
   int Nmax = param.N;
+  int Pmax = param.Pmax; //JEC new
   r_8 Rmax = param.R;
   string geom = param.geometry;
   int ntheta = param.ntheta;
@@ -237,15 +238,15 @@ void test3() {
   int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
   int Nshell = sphlagtrans.GetOrder();
   int Ntot = Nshell * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
-  //  int Npix = sphlagtrans.GetSphericalGeometry()->Npix();
-  //  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels
+  int Npix = sphlagtrans.GetSphericalGeometry()->Npix();
+  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels
 
-//   cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
-//        <<  Npix << " "
-//        <<  NpTot << " "
-//        <<  Nalm << " "
-//        <<  Nshell << " "
-//        <<  Ntot << endl;
+  cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
+       <<  Npix << " "
+       <<  NpTot << " "
+       <<  Nalm << " "
+       <<  Nshell << " "
+       <<  Ntot << endl;
 
   vector< complex<r_8> > flmn(Ntot);
   vector< complex<r_8> > almk(Ntot);
@@ -253,6 +254,84 @@ void test3() {
 
   tstack_pop("Lag Spherical Transf.");
 
+
+  //JEC 9/11/15 START
+  tstack_push("Laguerre 2 Bessel");
+  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
+  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
+
+  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax);
+  vector< complex<r_8> > FBlmp(Nalm*Pmax);
+  lag2bess.Lag2BessCoeff(flmn,FBlmp);
+  vector< complex<r_8> > FBalmk;
+  vector<r_8> fFBijk;
+  lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
+  tstack_pop("Laguerre 2 Bessel");
+
+  {//Check 1
+    cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
+  
+    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 << "n,lm: " << n << " ("<<l<<","<<m<<") : " << almk[id] << ", "<<FBalmk[id] << endl;
+	  
+	  complex<r_8> cdiff = (almk[id] - FBalmk[id])*conj(almk[id] -  FBalmk[id]);
+	  r_8 diff = sqrt(cdiff.real());
+	  if(diff>err_abs){
+	    err_abs = diff;
+	  }
+	  complex<r_8> cref= almk[id]*conj(almk[id]);
+	  r_8 ref = sqrt(cref.real());
+	  r_8 relatdiff = diff/ref;
+	  if(relatdiff>err_rel) err_rel = relatdiff;
+	}
+      }
+      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<n<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+      cout << "Almk: Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
+      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+    }
+  }//check 1
+
+
+
+  {//check 2
+
+    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;
+	r_8 diff = fabs(fijk[id]-fFBijk[id]);
+	if(diff>err_abs){
+	  err_abs = diff;
+	}
+	r_8 relatdiff = diff/ fijk[id];
+	if(relatdiff>err_rel) err_rel = relatdiff;
+      }//loop on px
+      
+            
+      cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<ish<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+      cout << "Fijk: Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
+      cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+
+    }//loop on shell
+
+
+  }//check 2
+
+
+  //JEC 9/11/15 END
+  
+
+
   tstack_push("Clk and Cln OUT");
   cout << "Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn" << endl;
 
@@ -528,8 +607,8 @@ int main(int narg, char *arg[]) {
 
   int N = 128; 
   r_8 R = 5.;
-  int Lmax = 1024;
-  int Pmax = 0;
+  int Lmax = 128;
+  int Pmax = 128;
 
   int test=3;