Skip to content

GitLab

Commit 3828e28f authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 23/11/15 computation of the Jlnp integral via a combined method: DFBT +... 

(JEC) 23/11/15 computation of the Jlnp integral via a combined method: DFBT + Clenshaw-Curtis quadrature
parent 7abb5300
Expand all Show whitespace changes
Inline Side-by-side
Showing with 2195 additions and 46 deletions
Makefile
View file @ 3828e28f
......@@ -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/lagSphericTransform.h
View file @ 3828e28f
......@@ -52,7 +52,7 @@ class LaguerreSphericalTransform {
LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2):
r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2, bool R2tau=true):
Lmax_(Lmax), Mmax_(Mmax), R_(R) {
//Factory
......@@ -68,7 +68,7 @@ class LaguerreSphericalTransform {
throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
}
lagTrans_ = new LaguerreTransform(Nmax, R, alpha);
lagTrans_ = new LaguerreTransform(Nmax, R, alpha, R2tau); //JEC 16/11/15 add R2tau
//do the map pixelization
......
src/lagsht_bessel.h
View file @ 3828e28f
......@@ -76,8 +76,8 @@ class Bessel {
//Helper
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;
/* 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_);
......
src/lagsht_testsuite.cc
View file @ 3828e28f
......@@ -18,7 +18,7 @@ using namespace std;
using namespace LagSHT;
#define DEBUG 1
#define DEBUG 2
//-------- Parameters set in the main and used in the different test functions
struct PARAM {
......@@ -233,11 +233,11 @@ void MultiSphericalLaguerreTransform() {
#endif
vector<r_8> fijk(NpTot);
sphlagtrans.Synthesis(flmnOrig,fijk);
#if DEBUG >= 2
//#if DEBUG >= 2
for (int i=0; i<NpTot; i++){
cout << "fijk("<<i<<"): " << fijk(i) << endl;
cout << "fijk("<<i<<"): " << fijk[i] << endl;
}
#endif
//#endif
tstack_pop("processing part Synthesis");
tstack_push("processing part Analysis");
......@@ -333,6 +333,12 @@ void TestLaguerre2Bessel() {
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++;
//verif
int idLag = n*Nalm + l+m*Lmax-m*(m+1)/2 ;
if(id != idLag)
throw LagSHTError("ERROR: FLag index error");
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);
......@@ -341,7 +347,23 @@ void TestLaguerre2Bessel() {
}//end loop on shell
#if DEBUG >=2
std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
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;
#endif
cout << "tstack 1 End" <<endl;
......@@ -361,6 +383,11 @@ void TestLaguerre2Bessel() {
vector< complex<r_8> > FBlmp(Nalm*Pmax);
lag2bess.Lag2BessCoeff(flmnOrig,FBlmp);
#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
cout << "tstack 2b End" <<endl;
tstack_pop("Compute Fourier-Bessel coeffs.");
......@@ -375,12 +402,18 @@ void TestLaguerre2Bessel() {
cout << "tstack 2d Start" <<endl;
tstack_push("Compute Alm(rk) from FB coeffs.");
tstack_push("Compute Alm(rk) from FL coeffs.");
vector< complex<r_8> > FLalmk;
vector<r_8> fFLijk;
sphlagtrans.Synthesis(flmnOrig,fFLijk,FLalmk);
#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
cout << "tstack 2d End" <<endl;
tstack_pop("Compute Alm(rk) from FB coeffs.");
tstack_pop("Compute Alm(rk) from FL coeffs.");
cout << "tstack 2 End" <<endl;
......
src/laguerre2bessel.cc
View file @ 3828e28f
......@@ -14,6 +14,7 @@ namespace LagSHT {
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.
*/
Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
LaguerreTransform* lagTrans,
......@@ -47,37 +48,54 @@ namespace LagSHT {
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(){
r_8 tau = lagTrans_->GetTau();
r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
cout << "JEC: " << tau3sqrt*sqrt(2.0/M_PI) << 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 = 3.0 * (rNm1-rNm2);
r_8 integLowerBound = rNm1;
r_8 integUpperBound = rNm1+deltaR; //for the integral
r_8 tol = 1e-10;
typedef ClenshawCurtisQuadrature<double> Integrator_t;
Integrator_t theQuad(20,"/Users/campagne/Travail/Software/Laguerre/LagSHT/data/ClenshawCurtisRuleData-20.txt",false); //false=do not recompute the weights and nodes of the quadrature
r_8 R = Rmax_ / tau; // R = r[N-1] voir acces direct = (lagTrans_->GetNodes()).back()
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(2/pi)*tau^(3/2)*sqrt(n!/(n+alpha)!) en iteratif
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 PPmax = max(512,Pmax_); // Sum maximum (DOIT ETRE > PMax pour la Synthese car appel a qlp
int PPmax = 1024; // Sum maximum To be put in parameter
int PPmaxm1 = PPmax -1;
jnu_ = new Bessel(Lmax_,PPmax); //To do put it in constructor
jnu_ = new Bessel(Lmax_,max(Pmax_,PPmax)); //To do put it in constructor
tstack_push("J Integ compute ");
......@@ -93,7 +111,7 @@ void Laguerre2Bessel::ComputeJInteg(){
r_8 qlpx = (*jnu_)(l,PPmaxm1); //ql,PPmax-1
r_8 Kapa_l = qlpx/R; // Kapa_l = (ql,PPmax-1)/R
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
......@@ -102,20 +120,34 @@ void Laguerre2Bessel::ComputeJInteg(){
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 * R; // rl,pp = (ql,pp/ql,PPmax-1) * R
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]);
//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));
//End
ofs <<l<<" "
<<n<<" "
<<p<< " " << setprecision(12) << Jlnp_[ploff7 + n*Jstride] << endl;
......@@ -135,8 +167,8 @@ void Laguerre2Bessel::ComputeJInteg(){
/*!
FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{ln}(k_{lp})
= Sqrt[2/Pi] Sum_n FL_{lmn} J_{lnp}
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
......@@ -147,6 +179,10 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
int Jstride = Lmax_*Pmax_;
std::ofstream ofs;
ofs.open ("FBlmp.txt", std::ofstream::out);
//try here a pedestrian computation
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
......@@ -164,10 +200,19 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
}//n
FBlmp[idBess] = sum;
ofs <<l<<" "
<<m<<" "
<<p<< " " << setprecision(12) << FBlmp[idBess] << endl;
}//m
}//l
}//p
ofs.close();
}//Transform
......@@ -221,22 +266,33 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
int id0 = k*Nalm_ + lmOff7;
complex<r_8> sum = 0.;
//JEC 18/11/15 simple summation algo
// complex<r_8> sum = 0.;
// for(int p=0;p<Pmax_;p++){
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);
r_8 qlp = (*jnu_)(l,p);
// int id1 = p*Nalm_ + lmOff7;
r_8 jlp1 = Bessel::jn(l+1,qlp);
// complex<r_8> tmp = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
// sum += tmp;
// }//loop on p
// FBalmk[id0] = sum;
//JEC 18/11/15 Kahan summation algo
vector< complex<r_8> > sumElmt(Pmax_);
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;
sum += (kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1]);
}//loop on p
FBalmk[id0] += sum;
sumElmt[p] = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
}
// 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
......
src/laguerre2bessel.h
View file @ 3828e28f
......@@ -27,12 +27,15 @@
#include <iostream>
#include <vector>
#include <complex>
#include <utility> //move
#include "lagsht_numbers.h"
#include "laguerreTransform.h"
#include "lagsht_spheregeom.h"
#include "lagsht_bessel.h"
#include "quadinteg.h" //JEC 12/11/15
using namespace std;
namespace LagSHT {
......@@ -41,6 +44,28 @@ namespace LagSHT {
//// --------- Class to transform Spherical Laguerre into Spherical Bessel -- //
////////////////////////////////////////////////////////////////////////////////
//implementation of the Kahan summation algorithm
template <class InIt>
typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
typedef typename std::iterator_traits<InIt>::value_type real;
real sum = real();
real running_error = real();
real temp;
real difference;
for (; begin != end; ++begin) {
difference = *begin;
difference -= running_error;
temp = sum;
temp += difference;
running_error = temp;
running_error -= sum;
running_error -= difference;
sum = std::move(temp);
}
return sum;
}
class Laguerre2Bessel {
public:
......@@ -87,6 +112,25 @@ private:
void ComputeJInteg();
/*!
function for numerical integration
*/
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)
......
src/laguerreBuilder.h
View file @ 3828e28f
......@@ -27,13 +27,14 @@
#include <cmath> //very important for abs definition
#include <vector>
#include <iostream>
#include <iomanip> // std::setprecision
#include "lagsht_numbers.h"
using namespace std;
#define DEBUG 0
#define DEBUG 1
namespace LagSHT {
......@@ -206,11 +207,15 @@ protected:
weights[k] = x*tmp1/(tmp0*tmp0);
}//end for
#if DEBUG >= 1
//verification
r_16 sumR = 0.;
r_16 sumW = 0.;
r_16 wip;
for (int i=0;i<N;i++){
cout << "root["<<i<<"]: " << setprecision(16) << nodes[i] << " LagN(ri) = " << Value(nodes[i]) << endl;
sumR += nodes[i];
wip = exp(log(abs(weights[i]))-nodes[i]); if(weights[i]<0) wip *= -1.0;
sumW += wip;
......@@ -220,7 +225,6 @@ protected:
for(int i=1;i<=alpha_; i++) alphaFact *= i;
#if DEBUG >= 1
r_16 sumRTh = (r_16)(N*(N+alpha_));
r_16 sumWTh = (r_16)(alphaFact);
......
src/quadinteg.h 0 → 100644
View file @ 3828e28f
This diff is collapsed.
test/test3D.cc
View file @ 3828e28f
......@@ -203,7 +203,8 @@ void test3() {
r_8 sigma = Lmax/20.; //fix
r_8 arg = (l-mean)/sigma;
arg *= arg;
Clkin(l,k) = exp(-arg/2.);
r_8 norm = exp(-((r_8)(k/Nmax)));
Clkin(l,k) = norm*exp(-arg/2.);
Clkin(l,k) /= (r_8)(2*l+1); //nomalisation
Clin[l] = Clkin(l,k);
}//l-loop
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment