lagsht_testsuite.cc 55 KB
Newer Older
1 2
#include <stdlib.h>
#include <iostream>
3
#include <iomanip>      // std::setprecision
4
#include <fstream>
5
#include <sstream>
6 7 8 9 10 11 12 13
#include <string.h>
#include <limits>       // std::numeric_limits
#include <string>
#include <typeinfo>
using namespace std;

#include "lagsht_exceptions.h"
#include "lagsht_utils.h" //getMemorySize
14
#include "walltimer.h"    //timing
15
#include <omp.h>          //OpenMP for sampling
16 17

#include "lagSphericTransform.h"
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
18 19
#include "laguerre2bessel.h"

20 21 22 23 24 25 26 27 28
//Test FFT
#include <vector>
#include <algorithm>
#include <functional>
#include <math.h>
#include <fftw3.h>



29 30
using namespace LagSHT;

31
#define DEBUG 10
32 33 34
#define DEBVEC 0
#define NUM_OMP_THREADS   4
#define CHECK_BY_CCQUAD 0
35 36 37

//JEC 12/1/16: alpha becomes double instead of int

38 39 40 41 42
//-------- Parameters set in the main and  used in the different test functions
struct PARAM {
  int Lmax;
  int N;
  r_8 R;
43
  int Pmax;
44
  int spin;
45 46
  //  int alpha;
  double alpha;
47 48 49
  string geometry;
  int ntheta;
  int nphi;
50 51 52
  string clenshawDir;
  string jlnpDir;
  bool recomputeJlnp;
53 54 55 56 57 58 59 60 61 62 63 64
} param ; 

//simple random generator (taken from sharp_testsuite
static double drand (double min, double max, int *state){
  *state = (((*state) * 1103515245) + 12345) & 0x7fffffff;
  return min + (max-min)*(*state)/(0x7fffffff+1.0);
}//rand


//----------------------------------------------------
void TestBasic() {

65
  // to access some general parameters... use param.
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92



  cout << " ___________  BASIC TEST _____________ " << endl;
  

  cout << " sizeof(double): " << sizeof(r_8) << " bytes" 
       << " min " << std::numeric_limits<r_8>::min() << "-ln(min): " << -log(std::numeric_limits<r_8>::min()) 
       << " max " << std::numeric_limits<r_8>::max()  
       << " min exp_10 " << std::numeric_limits<r_8>::min_exponent10
       << " max exp_10 " << std::numeric_limits<r_8>::max_exponent10  
       <<  endl;
  


  cout << " sizeof(long double): " << sizeof(r_16) << " bytes" 
       << " min " << std::numeric_limits<r_16>::min() <<  "-ln(min): " << -logl(std::numeric_limits<r_16>::min()) 
       << " max " << std::numeric_limits<r_16>::max()  
       << " min exp_10 " << std::numeric_limits<r_16>::min_exponent10
       << " max exp_10 " << std::numeric_limits<r_16>::max_exponent10  
       <<  endl;
}//TestBasic

//----------------
void LaguerreQuadrature() {

  int N = param.N;
93 94
  //  int alpha = param.alpha;
  double alpha = param.alpha;
95 96 97 98 99

  cout << " ___________  LaguerreQuadrature  TEST _____________ " << endl;

  LaguerreFuncQuad lagFunc(N,alpha);
  cout << "N = " << lagFunc.GetOrder() << " alpha = " << lagFunc.GetAlpha() << endl;
100 101
  vector<r_8> nodes;
  vector<r_8> weights;
102 103
  lagFunc.QuadWeightNodes(nodes,weights);

104 105 106 107 108 109 110 111
  {
    char buff1[80], buff2[80];
    ofstream ofs1, ofs2;
    sprintf(buff1,"lagNodes-%d-Func.txt",N);
    sprintf(buff2,"lagWeights-%d-Func.txt",N);
    ofs1.open(buff1,ofstream::out);
    ofs2.open(buff2,ofstream::out);
    for (int i=0;i<N;i++){
112 113
      ofs1 << setprecision(20) << nodes[i] << endl;
      ofs2 << setprecision(20) << weights[i] << endl;
114
    }
115 116
  }

117

118 119 120 121 122 123
}//LaguerreQuadrature

//--------------------------------------
void MultiLaguerreTransformSdtAlone() {
  int N = param.N;
  r_8 R  = param.R;
124
  r_8 alpha = param.alpha; //JEC 12/1/16 to be tested with alpha =/= 2 (eg. alpha = 0)
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144

  cout << " ___________   MultiLaguerreTransformSdtAlone  TEST _____________ " << endl;

  //    Timer tm("Test 3 (r_16)");
      
  int mapSize=2; //this is the stride
      
  int Ntot = N*mapSize;

#if DEBUG >= 1
  cout << "Start random generation...." << endl;
#endif
  vector< complex<r_8> > fnOrig(Ntot);
  int state=1234567 + 8912 ; //random seed
  for(int i=0;i<Ntot;i++){
    r_8 rv = drand(-1,1,&state);
    r_8 iv = drand(-1,1,&state);
    fnOrig[i] = complex<r_8>(rv,iv);
  }
#if DEBUG >=2 
145
  std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
146 147
#endif
    
148
  LaguerreTransform trans(N,R,alpha);  //JEC 12/1/16 add alpha
149 150 151 152 153 154 155 156

  vector< complex<r_8> > fi(Ntot);
#if DEBUG >= 1
  cout << "Start MultiSynthesis...." << endl;
#endif
  trans.MultiSynthesis(fnOrig,fi,mapSize);
    
#if DEBUG >=2 
157
  std::copy(fi.begin(), fi.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
158 159 160 161 162 163 164
#endif

  vector< complex<r_8> > fn(Ntot);
  trans.MultiAnalysis(fi,fn,mapSize);

  r_8 err_abs(0.);
  r_8 err_rel(0.);
165
  int imax=-1;
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
  for(int i=0;i<Ntot;i++){
    if(i<10)cout << "("<<i<<"): " << fnOrig[i] << " <-> " << fi[i] << " <-> "  
		 << fn[i] << endl;
    complex<r_8> cdiff = (fnOrig[i] -  fn[i]) * conj(fnOrig[i] -  fn[i]);
    r_16 diff = sqrt(cdiff.real());
    if(diff>err_abs){ 
      err_abs = diff;
      imax = i;
    }
    complex<r_8> foriCAbs = fnOrig[i]*conj(fnOrig[i]);
    r_8 foriAbs = sqrt(foriCAbs.real());
    r_8 relatdiff = diff/foriAbs;
	
    if((relatdiff)>err_rel) err_rel = relatdiff;
  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << "r_16 Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  

}//MultiLaguerreTransformSdtAlone

//---------------------------------------
void MultiSphericalLaguerreTransform() {

191
  r_8 alpha = param.alpha;  //JEC 13/1/16
192
  int Nmax = param.N;
193
  int Lmax = param.Lmax;
194
  r_8 Rmax = param.R;
195 196 197
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
198 199 200 201


  cout << " ___________  MultiSphericalLaguerreTransform   TEST _____________ " << endl;

202 203
  tstack_push("data input");

204 205
  int state = 1234567 + 8912 ; //random seed

206 207 208 209 210 211
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
212 213
					 nphi,
					 alpha
214
					 );
215

216 217 218

  int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
  int Nshell = sphlagtrans.GetOrder();
219
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
220
  int Npix =  sphlagtrans.GetSphericalGeometry()->Npix();
221
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels
222

223 224 225
  cout << "Verif: Ntheta, Nphi, Npix, Nptot, Nalm, Nshell, Ntot: "
       << sphlagtrans.GetSphericalGeometry()->NTheta() << " "
       << sphlagtrans.GetSphericalGeometry()->NPhi() << " "
226 227 228 229 230
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;
231

232 233 234 235 236 237 238
  unsigned int maxmemsize =  getMemorySize()/1e6;
  unsigned int estimateMem = 8*((uint_8)(NpTot+2*Ntot))/1e6;
  cout << "Estimate usage memory: " << estimateMem << " MBytes" << endl;
  if ( estimateMem > (0.9*maxmemsize) ) {
    cout << ">>>>> Warning: estimate Mem " << estimateMem <<" MB > 90\% "<< maxmemsize << " MB" << endl;
  }//test memory

239 240
  vector< complex<r_8> > flmnOrig(Ntot);
  int id=-1;
241
  for(int n=0;n<Nmax;n++){ //on each shell
242 243 244 245 246 247 248 249 250 251 252
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
	r_8 rv = drand(-1,1,&state);
	r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
	flmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
253
  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
254 255
#endif
      
256 257 258 259
  
  tstack_pop("data input");
  tstack_push("processing");
  tstack_push("processing part Synthesis");
260 261 262 263 264 265

#if DEBUG >= 1
  cout << "Main:Synthesis r_8 function..." << endl;
#endif
  vector<r_8> fijk(NpTot);
  sphlagtrans.Synthesis(flmnOrig,fijk);
266
#if DEBUG >= 2
267
  for (int i=0; i<NpTot; i++){
268
    cout << "fijk("<<i<<"): " << fijk[i] << endl;
269
  }
270
#endif
271 272 273 274
  
  tstack_pop("processing part Synthesis");
  tstack_push("processing part Analysis");
  
275 276 277 278 279 280
#if DEBUG >= 1
  cout << "Main:Analysis r_8 function..." << endl;
#endif
  vector< complex<r_8> > flmn(Ntot);
  sphlagtrans.Analysis(fijk,flmn);
    
281 282 283
  tstack_pop("processing part Analysis");
  tstack_pop("processing");  

284 285 286 287 288 289


  cout << "Error analysis..." << endl;

  r_8 err_abs(0.);
  r_8 err_rel(0.);
290
  int imax = -1;
291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307
  for(int i=0;i<Ntot;i++){
    if(i>(Ntot-9)) cout << "(" << i << ") : flnm Orig " <<  flmnOrig[i] << " <-> flmn Rec " <<  flmn[i] << endl;

    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 << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
308
  cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
309 310 311 312 313
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 


}//MultiSphericalLaguerreTransform

314 315 316
//---------------------------------------
void SpinMultiSphericalLaguerreTransform() {

317
  r_8 alpha = param.alpha;  
318 319 320 321 322 323
  int Nmax = param.N;
  int Lmax = param.Lmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
324
  int spin = param.spin; //JEC 18/1/16
325

326 327 328
  if(spin == 0)  
    throw LagSHTError("WARNING: spin =0 : use test MultiSphericalLaguerreTransform"); 

329
  cout << " ___________  SpinMultiSphericalLaguerreTransform   TEST _____________ " << endl;
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351

  tstack_push("data input");

  int state = 1234567 + 8912 ; //random seed

  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi,
					 alpha
					 );


  int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
  int Npix =  sphlagtrans.GetSphericalGeometry()->Npix();
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels

352 353 354 355

  cout << "Verif: Ntheta, Nphi, Npix, Nptot, Nalm, Nshell, Ntot: "
       << sphlagtrans.GetSphericalGeometry()->NTheta() << " "
       << sphlagtrans.GetSphericalGeometry()->NPhi() << " "
356 357 358 359 360 361
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;

362 363 364 365 366 367 368 369 370

  unsigned int maxmemsize =  getMemorySize()/1e6;
  unsigned int estimateMem = 2*8*((uint_8)(NpTot+2*Ntot))/1e6;
  cout << "Estimate usage memory: " << estimateMem << " MBytes" << endl;
  if ( estimateMem > (0.9*maxmemsize) ) {
    cout << ">>>>> Warning: estimate Mem " << estimateMem <<" MB > 90\% "<< maxmemsize << " MB" << endl;
  }//test memory


371 372
  vector< complex<r_8> > ElmnOrig(Ntot);
  vector< complex<r_8> > BlmnOrig(Ntot);
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
373 374
  int id
=-1;
375

376 377 378 379
  for(int n=0;n<Nmax;n++){ //on each shell
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
380 381
	r_8 rv = (l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	r_8 iv = (m==0 || l<abs(spin)) ? 0.0 : drand(-1,1,&state);
382
	ElmnOrig[id] = complex<r_8>(rv,iv);	    
383 384
	rv = (l<abs(spin)) ? 0.0 : drand(-1,1,&state);
	iv = (m==0 || l<abs(spin)) ? 0.0 : drand(-1,1,&state);
385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
	BlmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
  std::copy(ElmnOrig.begin(), ElmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
  std::copy(BlmnOrig.begin(), BlmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
      
  
  tstack_pop("data input");
  tstack_push("processing");
  tstack_push("processing part Synthesis");

#if DEBUG >= 1
  cout << "Main:Synthesis complex function..." << endl;
#endif
  vector<r_8> fijk_re; //NpTot
  vector<r_8> fijk_im;
405 406
//   vector< complex<r_8> > Elmk; //calcul intermediaire
//   vector< complex<r_8> > Blmk;
407

408 409
//  sphlagtrans.Synthesis(ElmnOrig, BlmnOrig, spin, fijk_re, fijk_im, Elmk, Blmk);
  sphlagtrans.Synthesis(ElmnOrig, BlmnOrig, spin, fijk_re, fijk_im);
410
#if DEBUG >= 2
411 412 413
  for (int i=0; i<NpTot; i++){
    cout << "fijk("<<i<<"): " << fijk_re[i] << ", " << fijk_im[i] << endl;
  }
414
#endif
415 416 417 418 419 420 421
  
  tstack_pop("processing part Synthesis");
  tstack_push("processing part Analysis");
  
#if DEBUG >= 1
  cout << "Main:Analysis complex function..." << endl;
#endif
422 423
  vector< complex<r_8> > Elmn(Ntot);
  vector< complex<r_8> > Blmn(Ntot);
424

425 426
//   vector< complex<r_8> > Elmk_ana;
//   vector< complex<r_8> > Blmk_ana;
427

428 429
//  sphlagtrans.Analysis(fijk_re, fijk_im, spin, Elmn, Blmn, Elmk_ana, Blmk_ana);
  sphlagtrans.Analysis(fijk_re, fijk_im, spin, Elmn, Blmn);
430 431 432 433 434 435 436 437 438 439 440
    
  tstack_pop("processing part Analysis");
  tstack_pop("processing");  


  cout << "Error analysis..." << endl;

  r_8 err_abs(0.);
  r_8 err_rel(0.);
  int imax = -1;
  for(int i=0;i<Ntot;i++){
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
441
    if(i>(Ntot-9)) {
442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
      cout << "(" << i << ") : Elnm Orig " <<  ElmnOrig[i] << " <-> Elmn Rec " <<  Elmn[i] << endl;
      cout << "............  Blnm Orig " <<  BlmnOrig[i] << " <-> Blmn Rec " <<  Blmn[i] << endl;
    }

    complex<r_8> cdiff = ((ElmnOrig[i] -  Elmn[i]) * conj(ElmnOrig[i] -  Elmn[i])
			  + (BlmnOrig[i] -  Blmn[i]) * conj(BlmnOrig[i] -  Blmn[i]))/2. ;
    r_16 diff = sqrt(cdiff.real());
    if(diff>err_abs){ 
      err_abs = diff;
      imax = i;
    }
    complex<r_8> foriCAbs = (ElmnOrig[i]*conj(ElmnOrig[i]) + BlmnOrig[i]*conj(BlmnOrig[i]))/2. ;
    r_8 foriAbs = sqrt(foriCAbs.real());
    r_8 relatdiff = diff/foriAbs;
	
    if((relatdiff)>err_rel) err_rel = relatdiff;

  }
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
  cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 


}//SpinMultiSphericalLaguerreTransform

467 468 469 470 471 472 473 474 475 476
//-------------------------------------------------------------
void TestJlnpComputation() {

  int Nmax = param.N;
  int Lmax = param.Lmax;
  int Pmax = param.Pmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
477 478 479
  string clenshawDir = param.clenshawDir;
  string jlnpDir = param.jlnpDir;
  bool recomputeJlnp = param.recomputeJlnp;
480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496

  cout << " ______________ TestJlnpComputation  START ___________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

  
  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();


  tstack_push("Ctor lag2bess");
497
  recomputeJlnp = true;
498 499 500 501
  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
			   clenshawDir,jlnpDir,recomputeJlnp);
  tstack_pop("Ctor lag2bess");

502 503 504



505 506 507
  cout << " ______________ TestJlnpComputation  End  ___________ " << endl;

}
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243
//---------------------------------------------------------------------------------
class JBess : public ClassFunc1D {
public:
  JBess(int l, r_8 kt):l_(l),kt_(kt) {norm_ = sqrt(2./M_PI);}
  virtual double operator()(double x) const {
    return norm_*Bessel::jn(l_,x*kt_);
  }
  virtual ~JBess() {}
private:
  int l_;
  r_8 kt_;
  r_8 norm_;
};
//------------------
class KLag : public ClassFunc1D {
public:
  KLag(LaguerreFuncQuad* Ln, r_8 norm):Ln_(Ln),norm_(norm){}
  virtual double operator()(double x) const {
    return norm_*x*x*Ln_->Value(x);
  }
  virtual ~KLag() {}
private:
  LaguerreFuncQuad* Ln_; //no ownership
  r_8 norm_;
};
//------------------
class IntFunc1D : public ClassFunc1D {
public:
  IntFunc1D(int l, r_8 klp,LaguerreFuncQuad* Ln, r_8 norm=1.0):
    l_(l),klp_(klp),Ln_(Ln),norm_(norm) {}
  virtual double operator()(double x) const {
    return norm_*(x*x)*(Bessel::jn(l_,x*klp_))*Ln_->Value(x);
  }
  virtual ~IntFunc1D() {}
private:
  int l_;
  r_8 klp_;
  LaguerreFuncQuad* Ln_; //no ownership
  r_8 norm_;
};//class IntFunc1D
//------------------
void ChebyshevSampling(int n, vector<r_8>& val, const ClassFunc1D& f, r_8 a, r_8 b){
  r_8 bma = 0.5*(b-a); r_8 bpa = 0.5*(b+a);
  r_8 cte = M_PI/((r_8)n);

  int k;
  r_8 x;
#pragma omp parallel for private(x) num_threads(NUM_OMP_THREADS)
  for(k=0;k<=n;k++){
    x = cos(k*cte)*bma+bpa;
    val[k] = f(x);
  }
}
//------------------
class FFTPlanning {
public:
  FFTPlanning(int n, vector<r_8>& vec) {
    createPlan(n,vec);
  }
  ~FFTPlanning() { DestroyPlan(); }
  void DestroyPlan() {
    if(plan_) {
      cout << "Destroy Plan ["<< this << "]" << endl;
      fftw_destroy_plan(plan_); 
      plan_=NULL;
    }
  }
  void Execute() const { fftw_execute(plan_); }
private:
  void createPlan(int n, vector<r_8>& vec) {
    cout << "Create Plan ["<< this << "]" << endl;
    plan_ = fftw_plan_r2r_1d(n + 1,&vec.data()[0], &vec.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
  }
  fftw_plan plan_;
};
//------------------
//size of val = n+1 
//void ChebyshevCoeffFFT(int n, const fftw_plan plan, vector<r_8>& val){
void ChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){


  tstack_push("FFTW plan exec...");
  //  fftw_execute ( plan );
  plan.Execute();
  tstack_pop("FFTW plan exec...");
  /*
    Chebyshev Coeff. the noramization of FFTW is propto  1/(val.size()-1) = 1/n
  */
  tstack_push("FFTW norm...");
  r_8 norm = 1./((r_8)n); 
  transform(val.begin(), val.end(), val.begin(),
		std::bind1st(multiplies<r_8>(),norm));

  val[0] *= 0.5;
  val[n] *= 0.5;
  tstack_pop("FFTW norm...");

}
//------------------
void InverseChebyshevCoeffFFT(int n, const FFTPlanning& plan, vector<r_8>& val){
  val[0] *= 2.0;
  val[n] *= 2.0;
  //  fftw_execute ( plan );
  plan.Execute();
  
  transform(val.begin(), val.end(), val.begin(),
		std::bind1st(multiplies<r_8>(),0.5));

}
//------------------
void VectorDebug(const string& msg, const vector<r_8>&vec, int beg, int end){
  cout << "(JEC) debug " << msg << endl;
  for(int i=beg;i<end;i++) {
    cout << vec[i] << ", ";
  }
  cout << endl;
}
//------------------
//extrait de quadinteg.h dans un premier temps.
//rappel l'original calcule les poids (et nodes) pour une quadrature sur [-1, 1]
void ClenshawCurtisWeights(int  order, vector<r_8>& w) throw(string) {
  double b;
  int i;
  int j;
  double pi = M_PI; //JEC use cmath definition of PI
  double theta;

  //Todo this kind of error should be tracked before
  if ( order < 1 ) {
    cerr << "\n";
    cerr << "CLENSHAW_CURTIS_COMPUTE - Fatal error!\n";
    cerr << "  Illegal value of ORDER = " << order << "\n";
    exit (EXIT_FAILURE);
  }
  
 
  if ( order == 1 ) {
    w[0] = 2.0;
  } else {
    
    for ( i = 0; i < order; i++ ) {
      theta = ( double ) ( i ) * pi / ( double ) ( order - 1 );  
      w[i] = 1.0;
      
      for ( j = 1; j <= ( order - 1 ) / 2; j++ ) {
	if ( 2 * j == ( order - 1 ) ) {
	  b = 1.0;
	} else {
	  b = 2.0;
	}
	
	w[i] = w[i] - b *  cos ( 2.0 * ( double ) ( j ) * theta )
	  / ( double ) ( 4 * j * j - 1 );
      }
    }
    
    w[0] = w[0] / ( double ) ( order - 1 );
    for ( i = 1; i < order - 1; i++ ) {
      w[i] = 2.0 * w[i] / ( double ) ( order - 1 );
    }
    w[order-1] = w[order-1] / ( double ) ( order - 1 );
  }
}

//weights defined for [-1 1]
//FFTW Y_k = 2 Sum''_j=0^(N-1) X_j Cos[Pi k j/(N-1)]
//CC weights
//W_k = 4/(N-1) a_k  Sum''_{j=0, j even}^(N-1) 1/(1-j^2) Cos[Pi k j/(N-1)]
void ClenshawCurtisWeightsFast(int n, const FFTPlanning&  plan, vector<r_8>& w) {

  //dim of w is n
  fill(w.begin(), w.end(), (r_8)0.0);

  for(int k=0;k<n; k +=2){
    w[k] = 1./(1.-(r_8)(k*k));
  }
  
  //  fftw_execute(plan);
  plan.Execute();
  

  r_8 norm = 2.0/(r_8)(n-1); //2 * FFTW DCT-1 
  transform(w.begin(), w.end(), w.begin(),
		std::bind1st(multiplies<r_8>(),norm));
  w[0] /= (r_8)2; w[n-1] /= (r_8)2;
  
}

/*
  estimates the order of the Chebyshev expansion approx using the number
  of roots and min/max of the cosine function asymptotic apprx of Spherical Bessel
  in the range [lowBnd, uppBnd]
 */
int EstimChebyshevOrdSpherBessel(int el, r_8 kt, r_8 lowBnd, r_8 uppBnd){
  int n=0;
  int qq = - ( el/2 +1 );
  r_8 val = M_PI/(2.0* kt) * (2*qq + el + 2 );
  while(val<lowBnd) {
    qq++;
    val = M_PI/(2.0* kt) * (2*qq + el + 2);
  }
  while(val<=uppBnd){
    qq++;
    n++;
    val = M_PI/(2.0* kt) * (2*qq + el + 2);
  }

  int nNodesAndMinMax = 2*n; //factor 2 to take into account min/max
  
  return floor(log2((r_8)nNodesAndMinMax))+1;
}
/*
  estimates the order of the Chebyshev expansion approx using the number
  of roots of Laguerre polynomial in the range [lowBnd, uppBnd]
 */
int EstimChebyshevOrdLagFunc(const vector<r_8>& nodes, r_8 lowBnd, r_8 uppBnd) {

  //If C++11 not available find an alternative with Class fucntion
  int n= count_if(nodes.begin(), nodes.end(), [&](r_8 x){ return (x>=lowBnd)&&(x<=uppBnd);} );

  if(n==0)n=1; //no root for Laguerre[0,x]

  int nNodesAndMinMax = 2*n; //factor 2 to take into account min/max
  return floor(log2((r_8)nNodesAndMinMax))+1;

}


//---------------------------------------------------------------------------------
void TestJnlpOrder() {
  int Nmax = param.N;
  int Lmax = param.Lmax;
  int Pmax = param.Pmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
  double alpha = param.alpha;
  cout << "_______________ TestJnlpOrder Start _____________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
  Bessel jnu(Lmax,std::max(Pmax,Nmax));  //the max is here just in case Nmax =/= Pmax

  r_8 tau = lagTrans->GetTau();
  r_8 rNm1 = Rmax / tau; // R = r[N-1] with N the original transform value. Todo (lagTrans_->GetNodes()).back()

  vector<r_8> integUpperBound(Nmax);
  integUpperBound[0] = min(100.,rNm1); //this is the case of L_0^{\alpha}(x) =1 with no root (L
  integUpperBound[1] = min(120.,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 << "FFT-order"<<Lmax<<"-N"<<Nmax<<"-P"<<Pmax<<".txt";
  std::ofstream ofs;
  string fname("./"+ss.str());
  ofs.open (fname.c_str(), std::ofstream::out);


  for(int p=0;p<Pmax;p++){
    for(int l=0; l<Lmax; l++){

      r_8 qlp = jnu(l,p);
      r_8 klptau = qlp/rNm1; 
      JBess jFunc(l,klptau);
  
      r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau;
      r_8 ql2scaled = jnu(l,1)/klptau;   //  the 2nd BesselRoot rescaled
  
      for (int n=0; n<Nmax; n++) {
  
	r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);
	
	int nJ = EstimChebyshevOrdSpherBessel(l,klptau,intLowBound,intUppBound);

	LaguerreFuncQuad lag(n,alpha);
	vector<r_8> nodes;
	vector<r_8> weights;
	lag.QuadWeightNodes(nodes,weights);

	int nK = EstimChebyshevOrdLagFunc(nodes,intLowBound,intUppBound);
	
	ofs << l << " " 
	     << n << " "
	     << p << " "
	     << nJ << " " 
	     << nK << endl;

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

  ofs.close();

  cout << "_______________ TestJnlpOrder End _____________ " << endl;

}//TestJnlpOrder
//---------------------------------------------------------------------------------
void TestJlnpComputationByFFT() {
  int Nmax = param.N;
  int Lmax = param.Lmax;
  int Pmax = param.Pmax;
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
  double alpha = param.alpha;

  string jlnpDir = param.jlnpDir;
  string clenshawDir = param.clenshawDir;
 
  cout << " ______________ TestJlnpComputation FFT START ___________ " << endl;
  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );


  tstack_push("Ctors....");  
  //  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
  Bessel jnu(Lmax,std::max(Pmax,Nmax));  //the max is here just in case Nmax =/= Pmax


  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()

  cout << "(JEC) ComputeJInteg Facts START " << endl;

  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)) );

  cout << "(JEC) ComputeJInteg Facts END " << endl;

  cout << "(JEC) ComputeJInteg  integUpperBound START " << endl;


  tstack_pop("Ctors....");  
  

  tstack_push("Bounds init...");  

  vector<r_8> integUpperBound(Nmax);
  integUpperBound[0] = min(100.,rNm1); //this is the case of L_0^{\alpha}(x) =1 with no root
  integUpperBound[1] = min(120.,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
  }
  tstack_pop("Bounds init...");  

  cout << "(JEC) ComputeJInteg  integUpperBound END " << endl;


  tstack_push("Common FFTW Plan init...");  

  cout << "(JEC) Plan Init START" << endl;

  //Try a common order for K & J function this make the FFTW plan for once
  int iOrdGen = 10; 
  int nOrdGen=pow(2.,(r_8)iOrdGen)-1; //the -1 is just to use FFTW with 2^iOrdGen 
  vector<r_8> VecDCT1(nOrdGen+1,0.);  //but for  FFTW_REDFT00 it is requiered n+1... see FFTW doc
  FFTPlanning planJK(nOrdGen,VecDCT1);

  cout << "(JEC) Plan Init END" << endl;

  cout << "(JEC) Plan Inverse Init START" << endl;

  //Unique Product plan
  int nOrdProd= 2*nOrdGen+1; //order max of product is 2*nOrdGen but to get a power of 2 for the FFT I add 1.
  vector<r_8> VecDCT1Inv(nOrdProd+1,0.);
  FFTPlanning planInv(nOrdProd,VecDCT1Inv);

  cout << "(JEC) Plan Inverse Init END" << endl;

  tstack_pop("Common FFTW Plan init...");  
  
  tstack_push("Common CC weights init...");  

  cout << "(JEC) CC weights START" << endl;

  //Clenshow-Curtis single quadratuer weight
  vector<r_8> wCC(nOrdProd+1,0);
  FFTPlanning planCC(nOrdProd, wCC);

  ClenshawCurtisWeightsFast(nOrdProd+1,planCC,wCC);
#if DEBVEC>0
  VectorDebug("DCT CCurtis W [0,10]: ",wCC,0,10);
  VectorDebug("DCT CCurtis W [last-10,last]: ",wCC,nOrdProd-11,nOrdProd+1);
#endif



  cout << "(JEC) CC weights END" << endl;

  tstack_pop("Common CC weights init...");  

  stringstream ss; ss << "OMP-FFTClass-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 p=0;p<10;p++){
    for(int l=0; l<Lmax; l++){
      //for(int l=0; l<10; l++){

      r_8 qlp = jnu(l,p);

      r_8 klptau = qlp/rNm1; 
      JBess jFunc(l,klptau);
  
      //---------> r_8 klptau =1.; //JUST FOR TEST wrt M10
  

      r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau;
      r_8 ql2scaled = jnu(l,1)/klptau;   //  the 2nd BesselRoot rescaled
  


      for (int n=0; n<Nmax; n++) {
	//for (int n=0; n<10; n++) {
  
	r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);
	
#if DEBUG>10

	cout << "lowBnd, uppBnd: " << setprecision(30) 
	     << intLowBound << ", " << intUppBound 
	     << setprecision(6)
	     << endl;
	cout << "l,n,p,klptau: " 
	     << l << ", " 
	     << n << ", " 
	     << p << ", " 
	     << klptau << endl;
#endif
  
	tstack_push("Integral Computation...");
  
	tstack_push("Chebyshev J...");  


	//Get Chebyshev coeff of Bessel
	//  int iOrdJ = 8; int nOrdJ=pow(2.,(r_8)iOrdJ);
	int nOrdJ = nOrdGen; //USE COMMON dimension
	vector<r_8> coefJ(nOrdJ+1,0.);
	{
	  tstack_push("Chebyshev J SAMPLING...");  
	  fill(VecDCT1.begin(), VecDCT1.end(), 0.0); //may be unusefull here
	  ChebyshevSampling(nOrdJ,VecDCT1,jFunc,intLowBound,intUppBound);

#if DEBVEC>0
	  VectorDebug("Sampling J-func",VecDCT1,0,10);
#endif
	  tstack_pop("Chebyshev J SAMPLING...");  
	  tstack_push("Chebyshev J FFT...");  

	  ChebyshevCoeffFFT(nOrdJ,planJK, VecDCT1);
	  copy(VecDCT1.begin(),VecDCT1.end(),coefJ.begin());

	  tstack_pop("Chebyshev J FFT...");  

#if DEBVEC>0
	  VectorDebug("Che J-func",coefJ,0,10);
#endif
	}

	tstack_pop("Chebyshev J...");  


	tstack_push("Chebyshev K...");  
  
	//Get Chebyshev coeff of Laguerre Func
	//  int iOrdK = 6; int nOrdK = pow(2.,(r_8)iOrdK);
	int nOrdK =  nOrdGen; //USE COMMON dimension
	vector<r_8> coefK(nOrdK+1,0.);
	{
	  tstack_push("Chebyshev K SAMPLING..."); 
	  LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
	  KLag kLagFunc(Ln, facts[n]);
	  fill(VecDCT1.begin(), VecDCT1.end(), 0.0); //may be unusefull here
	  ChebyshevSampling(nOrdK,VecDCT1,kLagFunc,intLowBound,intUppBound);

#if DEBVEC>0
	  VectorDebug("Sampling K-func",VecDCT1,0,10);
#endif
	  tstack_pop("Chebyshev K SAMPLING..."); 
	  tstack_push("Chebyshev K FFT..."); 
	  ChebyshevCoeffFFT(nOrdK, planJK, VecDCT1);
	  copy(VecDCT1.begin(),VecDCT1.end(),coefK.begin());

	  tstack_pop("Chebyshev K FFT..."); 
#if DEBVEC>0
	  VectorDebug("Che K-func",coefK,0,10);
#endif

	  delete Ln;
	}
	tstack_pop("Chebyshev K...");  
  

	tstack_push("Chebyshev Product...");  

	//Product size
	//   int nOrdProd = nOrdJ*nOrdK;
	//perform Inverse Chebyshev transform in the new basis

	tstack_push("Inverse Chebyshev J...");  

	vector<r_8>coefJExt(nOrdProd+1,0.); 
	fill(VecDCT1Inv.begin(), VecDCT1Inv.end(), 0.0);
	copy(coefJ.begin(),coefJ.end(),VecDCT1Inv.begin());  //make the VecDCT1Inv[0:nOrdJ-1] = coefJ; then VecDCT1Inv is comleted by 0
	InverseChebyshevCoeffFFT(nOrdProd,planInv,VecDCT1Inv);
	copy(VecDCT1Inv.begin(), VecDCT1Inv.end(), coefJExt.begin());

#if DEBVEC>0
	VectorDebug("Inv Che J",coefJExt,0,10);
#endif

	tstack_pop("Inverse Chebyshev J...");  
	tstack_push("Inverse Chebyshev K...");  

	vector<r_8>coefKExt(nOrdProd+1,0.); 
	fill(VecDCT1Inv.begin(), VecDCT1Inv.end(), 0.0);
	copy(coefK.begin(),coefK.end(),VecDCT1Inv.begin());
	InverseChebyshevCoeffFFT(nOrdProd,planInv,VecDCT1Inv);
	copy(VecDCT1Inv.begin(), VecDCT1Inv.end(), coefKExt.begin());

#if DEBVEC>0
	VectorDebug("Inv Che K",coefKExt,0,10);
#endif

	tstack_pop("Inverse Chebyshev K...");  


	tstack_push("Inverse Chebyshev J*K...");  

	vector<r_8>invCoefProd(nOrdProd+1,0.); //is there an in-place product?
	transform(coefJExt.begin(),coefJExt.end(),
		  coefKExt.begin(),
		  invCoefProd.begin(),
		  multiplies<r_8>());

#if DEBVEC>0
	VectorDebug("Inv Che Prod",invCoefProd,0,10);
#endif
  
	tstack_pop("Inverse Chebyshev J*K...");  
	tstack_pop("Chebyshev Product...");  
  
	//Compute integral
	tstack_push("CC quadrature...");
	r_8 integral = inner_product(invCoefProd.begin(),invCoefProd.end(),wCC.begin(),0.);
	integral *= (intUppBound - intLowBound)/2.0; //the 2 division comme from CC quadrature weights definition in [-1,1];
	tstack_pop("CC quadrature...");

#if CHECK_BY_CCQUAD>0
	cout << "Integral via FFT               = " << setprecision(20) << integral <<  setprecision(6) << endl;
#endif
	ofs <<l<<" "
	    <<n<<" "
 	    <<p<< " " << setprecision(30) << integral 
	    << endl;

	tstack_pop("Integral Computation...");


#if CHECK_BY_CCQUAD>0
	tstack_push("Recursive CC quadrature");
 
	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;
	r_8 tol4Integral = 1e-10;
  
	LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
	IntFunc1D f(l,klptau,Ln);
	theQuad.SetFuncBounded(f,intLowBound,intUppBound);
	integ_val = theQuad.GlobalAdapStrat(tol4Integral);
	delete Ln;
	tstack_pop("Recursive CC quadrature");

	cout << "Integral via Recursive CC Quad = " 
	     << setprecision(20) 
	     << integ_val.first * sqrt(2.0/M_PI) * facts[n]
	     << setprecision(6) 
	     << endl;
#endif

      }//n-loop
    }//l-loop
  }//p-loop

  //file to save computations
  ofs.close();

  cout << " ______________ TestJlnpComputation FFT END ___________ " << endl;


}


//---------------------------------------------------------------------------------
// 
//Version 0 du calcul des Jlnp par FFTW mais non optimisee pour les plans.
// 
//
// void ChebyshevCoeffFFT(int n, vector<r_8>& val){
//   /* Planning In-place
//    */
//   tstack_push("FFTW plan init...");
//   fftw_plan plan;
//   plan = fftw_plan_r2r_1d(n + 1,&val.data()[0], &val.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
//   tstack_pop("FFTW plan init...");

//   tstack_push("FFTW plan exec...");
//   fftw_execute ( plan );
//   tstack_pop("FFTW plan exec...");
//   /*
//     Chebyshev Coeff.
//   */
//   tstack_push("FFTW norm...");
//   r_8 norm = 1./((r_8)n);
//   transform(val.begin(), val.end(), val.begin(),
// 		std::bind1st(multiplies<r_8>(),norm));
//   //  for(int k=0;k<=n;k++) val[k] *= norm;
//   val[0] *= 0.5;
//   val[n] *= 0.5;
//   tstack_pop("FFTW norm...");

//   fftw_destroy_plan(plan);
// }
// //------------------
// void InverseChebyshevCoeffFFT(int n, vector<r_8>& val){
//   val[0] *= 2.0;
//   val[n] *= 2.0;
//   //  for(int k=0;k<=n;k++) val[k] *= n; //not necessary
//   fftw_plan plan;
//   plan = fftw_plan_r2r_1d(n + 1,&val.data()[0], &val.data()[0], FFTW_REDFT00, FFTW_ESTIMATE);
//   fftw_execute ( plan );
//   //  for(int k=0;k<=n;k++) val[k] /= (2.0*n); //not necessary
//   for(int k=0;k<=n;k++) val[k] *= 0.5;
//   fftw_destroy_plan(plan);
// }
// //------------------
// void VectorDebug(const string& msg, const vector<r_8>&vec, int beg, int end){
//   cout << "(JEC) debug " << msg << endl;
//   for(int i=beg;i<end;i++) {
//     cout << vec[i] << ", ";
//   }
//   cout << endl;
// }
// void TestJlnpComputationByFFT() {
//   int Nmax = param.N;
//   int Lmax = param.Lmax;
//   int Pmax = param.Pmax;
//   r_8 Rmax = param.R;
//   string geometry = param.geometry;
//   int ntheta = param.ntheta;
//   int nphi = param.nphi;
//   double alpha = param.alpha;
//   string clenshawDir = param.clenshawDir;
 
//   cout << " ______________ TestJlnpComputation FFT START ___________ " << endl;
//   LaguerreSphericalTransform sphlagtrans(geometry,
// 					 Lmax,
// 					 -1,
// 					 Nmax,
// 					 Rmax,
// 					 ntheta,
// 					 nphi
// 					 );


//   tstack_push("Ctors....");  
//   //  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
//   LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
//   Bessel jnu(Lmax,std::max(Pmax,Nmax));  //the max is here just in case Nmax =/= Pmax


//   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()

//   cout << "(JEC) ComputeJInteg Facts START " << endl;

//   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)) );

//   cout << "(JEC) ComputeJInteg Facts END " << endl;

//   cout << "(JEC) ComputeJInteg  integUpperBound START " << endl;


//   tstack_pop("Ctors....");  
  

//   tstack_push("Bounds init...");  

//   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
//   }
//   tstack_pop("Bounds init...");  

//   cout << "(JEC) ComputeJInteg  integUpperBound END " << endl;
1244

1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402
//   tstack_push("Integral Computation...");

//   int p=256;
//   int l=64;
//   int n=32;
//   r_8 qlp = jnu(l,p);

//   //----------->   r_8 klptau = qlp/rNm1; 
  
//   r_8 klptau =1.; //JUST FOR TEST wrt M10
  
//   cout << "n,l,p,klptau: " 
//        << n << ", " 
//        << l << ", " 
//        << p << ", " 
//        << klptau << endl;

//   r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau;
//   r_8 ql2scaled = jnu(l,1)/klptau;   //  the 2nd BesselRoot rescaled
  
//   r_8 intUppBound = std::max(integUpperBound[n], ql2scaled);

//   cout << "lowBnd, uppBnd: " << setprecision(30) 
//        << intLowBound << ", " << intUppBound 
//        << setprecision(6)
//        << endl;

//   tstack_push("Chebyshev J...");  


//   //Get Chebyshev coeff of Bessel
//   int iOrdJ = 8; int nOrdJ=pow(2.,(r_8)iOrdJ);
//   vector<r_8> coefJ(nOrdJ+1,0.);
//   {
//     JBess jFunc(l,klptau);
//     tstack_push("Chebyshev J SAMPLING...");  
//     ChebyshevSampling(nOrdJ,coefJ,jFunc,intLowBound,intUppBound);

// #if DEBVEC>0
//     VectorDebug("Sampling J-func",coefJ,0,10);
// #endif
//     tstack_pop("Chebyshev J SAMPLING...");  
//     tstack_push("Chebyshev J FFT...");  

//     ChebyshevCoeffFFT(nOrdJ,coefJ);

//     tstack_pop("Chebyshev J FFT...");  

// #if DEBVEC>0
//     VectorDebug("Che J-func",coefJ,0,10);
// #endif
//   }

//   tstack_pop("Chebyshev J...");  
//   tstack_push("Chebyshev K...");  
  
//   //Get Chebyshev coeff of Laguerre Func
//   int iOrdK = 6; int nOrdK = pow(2.,(r_8)iOrdK);
//   vector<r_8> coefK(nOrdK+1,0.);
//   {  
//     LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
//     KLag kLagFunc(Ln, facts[n]);
//     ChebyshevSampling(nOrdK,coefK,kLagFunc,intLowBound,intUppBound);

// #if DEBVEC>0
//     VectorDebug("Sampling K-func",coefK,0,10);
// #endif

//     ChebyshevCoeffFFT(nOrdK,coefK);

// #if DEBVEC>0
//     VectorDebug("Che K-func",coefK,0,10);
// #endif

//     delete Ln;
//   }
//    tstack_pop("Chebyshev K...");  
  

//    tstack_push("Chebyshev Product...");  

//   //Product size
//    int nOrdProd = nOrdJ*nOrdK;
//   //extend the Chebyshev coeff vector with 0. Is there a simpler way?
//   vector<r_8>coefJExt(nOrdProd+1,0.); copy(coefJ.begin(),coefJ.end(),coefJExt.begin());
//   vector<r_8>coefKExt(nOrdProd+1,0.); copy(coefK.begin(),coefK.end(),coefKExt.begin());
//   //perform Inverse Chebyshev transform in the new basis
//   InverseChebyshevCoeffFFT(nOrdProd,coefJExt);

// #if DEBVEC>0
//   VectorDebug("Inv Che J",coefJExt,0,10);
// #endif

//   InverseChebyshevCoeffFFT(nOrdProd,coefKExt);

// #if DEBVEC>0
//   VectorDebug("Inv Che K",coefKExt,0,10);
// #endif

//   vector<r_8>invCoefProd(nOrdProd); //is there an in-place product?
//   transform(coefJExt.begin(),coefJExt.end(),
// 	    coefKExt.begin(),
// 	    invCoefProd.begin(),
// 	    multiplies<r_8>());

// #if DEBVEC>0
//   VectorDebug("Inv Che Prod",invCoefProd,0,10);
// #endif
  
//   tstack_pop("Chebyshev Product...");  
  
//   tstack_push("CCurtis weights...");
//   int iOrdCC = iOrdJ+iOrdK; int nOrdCC = pow(2.0,(r_8)iOrdCC)+1;
//   vector<r_8> wCC(nOrdCC,0);
//   ClenshawCurtisWeights(nOrdCC,wCC);

// #if DEBVEC>0
//   VectorDebug("CCurtis W [0,10]: ",wCC,0,10);
//   VectorDebug("CCurtis W [last-10,last]: ",wCC,nOrdCC-11,nOrdCC);
// #endif
//   tstack_pop("CCurtis weights...");
  
//   //Compute integral
//   tstack_push("CC quadrature...");
//   r_8 integral = inner_product(invCoefProd.begin(),invCoefProd.end(),wCC.begin(),0.);
//   integral *= (intUppBound - intLowBound)/2.0; //the 2 division comme from CC quadrature weights definition in [-1,1];
//   tstack_pop("CC quadrature...");
//   cout << "Integral via FFT               = " << setprecision(20) << integral <<  setprecision(6) << endl;

//   tstack_pop("Integral Computation...");


//   tstack_push("Recursive CC quadrature");
 
//   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;
//   r_8 tol4Integral = 1e-10;
  
//   LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n, alpha);
//   IntFunc1D f(l,klptau,Ln);
//   theQuad.SetFuncBounded(f,intLowBound,intUppBound);
//   integ_val = theQuad.GlobalAdapStrat(tol4Integral);
//   delete Ln;
//   tstack_pop("Recursive CC quadrature");

//   cout << "Integral via Recursive CC Quad = " 
//        << setprecision(20) 
//        << integ_val.first * sqrt(2.0/M_PI) * facts[n]
//        << setprecision(6) 
//        << endl;
  
  
//   cout << " ______________ TestJlnpComputation FFT END ___________ " << endl;


// }
1403

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1404 1405 1406 1407 1408 1409 1410
//-------------------------------------------------------------
void TestLaguerre2Bessel() {

  //Todo introduce Pmax !

  int Nmax = param.N;
  int Lmax = param.Lmax;
1411
  int Pmax = param.Pmax;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1412 1413 1414 1415
  r_8 Rmax = param.R;
  string geometry = param.geometry;
  int ntheta = param.ntheta;
  int nphi = param.nphi;
1416 1417 1418 1419
  string clenshawDir = param.clenshawDir;   //JEC 23/11/15
  string jlnpDir = param.jlnpDir;           //       "
  bool recomputeJlnp = param.recomputeJlnp; //       "

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1420 1421 1422

  cout << " ______________ TestLaguerre2Bessel TEST ___________ " << endl;
  
1423 1424
   cout << "tstack 1 Start" <<endl;
   tstack_push("data input");
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459

  int state = 1234567 + 8912 ; //random seed

  LaguerreSphericalTransform sphlagtrans(geometry,
					 Lmax,
					 -1,
					 Nmax,
					 Rmax,
					 ntheta,
					 nphi
					 );

  
  BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
  LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();

  int Nalm = sphere->Nalm();
  int Nshell = sphlagtrans.GetOrder();
  int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
  int Npix =  sphere->Npix();
  int NpTot=  Nshell * Npix; //totoal number of 3D-pixels

  cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
       <<  Npix << " "
       <<  NpTot << " "
       <<  Nalm << " "
       <<  Nshell << " "
       <<  Ntot << endl;

  vector< complex<r_8> > flmnOrig(Ntot);
  int id=-1;
  for(int n=0;n<Nmax;n++){ //on each shell
    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
      for(int l=m;l<Lmax;l++){
	id++;
1460 1461 1462 1463 1464 1465
	
	//verif
	int idLag =  n*Nalm + l+m*Lmax-m*(m+1)/2 ;
	if(id != idLag)
	  throw LagSHTError("ERROR: FLag index error"); 

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1466 1467 1468 1469 1470 1471 1472 1473
	r_8 rv = drand(-1,1,&state);
	r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
	flmnOrig[id] = complex<r_8>(rv,iv);	    
      }//end l-loop
    }//end m-loop
  }//end loop on shell

#if DEBUG >=2 
1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490
  cout << "F-Lag coeff orig ......... START " << endl;
  //  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  
  for(int l=0; l<Lmax; l++){
    for (int m=0; m<=l; m++) {
      int almoff7 = l+m*Lmax-m*(m+1)/2;
      for (int n=0; n<Nmax; n++){
	int idLag =  n*Nalm + almoff7;
	cout << l << " " 
	     << m << " "
	     << n << " "
	     << setprecision(12) << flmnOrig[idLag] << endl;
      }
    }
  }
  
  cout << "F-Lag coeff orig ......... END " << endl;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1491 1492 1493 1494 1495 1496 1497 1498 1499 1500
#endif
        
  cout << "tstack 1 End" <<endl;
  tstack_pop("data input");


  cout << "tstack 2 Start" <<endl;
  tstack_push("processing");
  cout << "tstack 2a Start" <<endl;
  tstack_push("Init Laguerre 2 Bessel");
1501 1502
  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
			   clenshawDir,jlnpDir,recomputeJlnp);
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1503 1504
  cout << "tstack 2a End" <<endl;
  tstack_pop("Init Laguerre 2 Bessel");
1505

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1506 1507 1508 1509
  cout << "tstack 2b Start" <<endl;
  tstack_push("Compute Fourier-Bessel coeffs.");

  vector< complex<r_8> > FBlmp(Nalm*Pmax);
1510
  lag2bess.Lag2BessCoeff(flmnOrig,FBlmp);
1511 1512 1513 1514 1515
#if DEBUG >=2 
  cout << "F-Bessel coeff ......... START " << endl;
  std::copy(FBlmp.begin(), FBlmp.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  cout << "F-Bessel coeff ......... END " << endl;
#endif
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1516 1517 1518 1519 1520 1521 1522
  cout << "tstack 2b End" <<endl;
  tstack_pop("Compute Fourier-Bessel coeffs.");
  

  cout << "tstack 2c Start" <<endl;
  tstack_push("Compute Alm(rk) from FB coeffs.");
  vector< complex<r_8> > FBalmk;
1523
  vector<r_8> fFBijk;
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1524 1525 1526 1527 1528 1529
  lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
  cout << "tstack 2c End" <<endl;
  tstack_pop("Compute Alm(rk) from FB coeffs.");
  

  cout << "tstack 2d Start" <<endl;
1530
  tstack_push("Compute Alm(rk) from FL coeffs.");
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1531 1532 1533
  vector< complex<r_8> > FLalmk;
  vector<r_8> fFLijk;
  sphlagtrans.Synthesis(flmnOrig,fFLijk,FLalmk);
1534 1535 1536 1537 1538 1539
#if DEBUG >=2 
  cout << "F-Lag fijk coeff ......... START " << endl;
  std::copy(fFLijk.begin(), fFLijk.end(), std::ostream_iterator< complex<r_8> >(std::cout, "\n"));
  cout << "F-Lag fijk coeff ......... END " << endl;
#endif

Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1540
  cout << "tstack 2d End" <<endl;
1541
  tstack_pop("Compute Alm(rk) from FL coeffs.");
Jean-Eric Campagne's avatar
Jean-Eric Campagne committed
1542
  
1543
  cout << "tstack 2 End" <<endl;
1544

1545 1546 1547
  tstack_pop("processing");


1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577
  //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");


1578 1579