📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

Commit 51122e08 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 13/1/16 extend and fix bugs when alpha =/= 2

parent 72f5ae1d
......@@ -34,7 +34,7 @@ endif
# ===== Libsharp =====
#SHARPDIR = <yourdir>/libsharp-code/auto/
SHARPDIR = /exp/opera/kits/libsharp-code/auto
SHARPDIR = /Users/campagne/Travail/kits/libsharp-code/auto
SHARPLIB = $(SHARPDIR)/lib
SHARPINC = -I$(SHARPDIR)/include
SHARPLIBN = -L$(SHARPLIB) -lsharp -lc_utils -lfftpack
......
......@@ -21,13 +21,17 @@ using namespace LagSHT;
#define DEBUG 0
//JEC 12/1/16: alpha becomes double instead of int
//-------- Parameters set in the main and used in the different test functions
struct PARAM {
int Lmax;
int N;
r_8 R;
int Pmax;
int alpha;
// int alpha;
double alpha;
string geometry;
int ntheta;
int nphi;
......@@ -47,9 +51,6 @@ static double drand (double min, double max, int *state){
void TestBasic() {
// to access some general parameters... use param.
// int N = param.N;
// int alpha = param.alpha;
// int Lmax = param.Lmax;
......@@ -77,7 +78,8 @@ void TestBasic() {
void LaguerreQuadrature() {
int N = param.N;
int alpha = param.alpha;
// int alpha = param.alpha;
double alpha = param.alpha;
cout << " ___________ LaguerreQuadrature TEST _____________ " << endl;
......@@ -107,6 +109,7 @@ void LaguerreQuadrature() {
void MultiLaguerreTransformSdtAlone() {
int N = param.N;
r_8 R = param.R;
r_8 alpha = param.alpha; //JEC 12/1/16 to be tested with alpha =/= 2 (eg. alpha = 0)
cout << " ___________ MultiLaguerreTransformSdtAlone TEST _____________ " << endl;
......@@ -130,7 +133,7 @@ void MultiLaguerreTransformSdtAlone() {
std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
LaguerreTransform trans(N,R);
LaguerreTransform trans(N,R,alpha); //JEC 12/1/16 add alpha
vector< complex<r_8> > fi(Ntot);
#if DEBUG >= 1
......@@ -628,11 +631,15 @@ int main(int narg, char *arg[]) {
bool recomputeJlnp = false;
//JEC 12/1/16 alpha as input
double alpha = 2.0; //default
int ka=1;
while (ka<narg) {
if (strcmp(arg[ka],"-h")==0) {
cout << "usage: ./Objs/lagsht_testsuite -t <test number> [1]\n"
<< " [-n <Nmax>] [-l <Lmax>]\n"
<< " [-a <alpha> [2 by default] ] "
<< " [-n <Nmax> [128]] [-l <Lmax> [128]]\n"
<< " [-g <geometry> Gauss|Fejer1|Healpix [Gauss]] \n"
<< " [-ntheta <number of theta rings> [determined by the geometry]\n in case of Healpix gives the nside parameter]\n"
<< " [-nphi <number of pixel per rings> [determined by the geometry]]\n"
......@@ -641,6 +648,10 @@ int main(int narg, char *arg[]) {
<< endl;
return 0;
}
else if (strcmp(arg[ka],"-a")==0) {
alpha = atof(arg[ka+1]);
ka+=2;
}
else if (strcmp(arg[ka],"-n")==0) {
N=atoi(arg[ka+1]);
ka+=2;
......@@ -692,7 +703,7 @@ int main(int narg, char *arg[]) {
param.N = N;
param.Pmax = Pmax;
param.R = R;
param.alpha = 2;
param.alpha = alpha;
param.geometry = geometry;
param.ntheta = ntheta;
param.nphi = nphi;
......
......@@ -10,12 +10,12 @@ void LaguerreFuncQuad::QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights)
}//LaguerreFuncQuad::QuadWeightNodes
r_8 LaguerreFuncQuad::Value(r_8 x, r_8 scale) const{
return value(x,N_);
return value(x,N_,scale); //JEC 12/1/16
}//LaguerreFuncQuad::Value
void LaguerreFuncQuad::Values(r_8 x, vector<r_8>& vals, r_8 scale) {
if(vals.size() != (size_t)(N_+1) ) throw LagSHTError("Values with vals.size != (N+1)");
values(x,N_,vals);
values(x,N_,vals, scale); //JEC 12/1/16
}//LaguerreFuncQuad::Values
......
......@@ -29,6 +29,10 @@
#include <iostream>
#include <iomanip> // std::setprecision
//JEC 12/1/16: alpha becomes double instead of int
//JEC (N+alpha)! -> Gamma[N+alpha+1]. Use Boost::tgamma function
#include "boost/math/special_functions/gamma.hpp"
#include "lagsht_numbers.h"
using namespace std;
......@@ -49,7 +53,7 @@ namespace LagSHT {
class LaguerreFuncQuad {
public:
//! Creator
LaguerreFuncQuad(int N, int alpha=2): N_(N), alpha_(alpha),
LaguerreFuncQuad(int N, r_8 alpha=2): N_(N), alpha_(alpha),
ln2(0.693147180559945309417232121458176),
inv_ln2(1.4426950408889634073599246810018921),
fsmall(std::ldexp(1.,-large_exponent2)),
......@@ -70,7 +74,7 @@ public:
//! Getters
int GetOrder() const { return N_; }
int GetAlpha() const { return alpha_; }
r_8 GetAlpha() const { return alpha_; }
//! Compute Weights and Nodes of the Gauss-Laguerre quadrature
void QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) const;
......@@ -86,15 +90,15 @@ protected:
/*! returns the amount by which x must be shifted to go closer to a root of
L_N^alpha(x)
*/
double newtoncorr (double x, int N) const {
r_8 newtoncorr (r_8 x, int N) const {
using namespace std;
if(N<1) return 0.;
double pnm1=0., pn=1.;
double c1=alpha_-1., c2=alpha_-1.;
r_8 pnm1=0., pn=1.;
r_8 c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
r_8 pnm2 = pnm1;
pnm1 = pn;
c1+=2.;
c2+=1.;
......@@ -104,25 +108,26 @@ protected:
while (abs(pn)>=fbig)
{ pn*=fsmall; pnm1*=fsmall; }
}
double ppn=(N*pn-(N+alpha_)*pnm1)/x;
r_8 ppn=(N*pn-(N+alpha_)*pnm1)/x;
return -pn/ppn;
}//newtoncorr
/*! returns L_N^alpha(x)*exp(-x/2)
//JEC 12/1/16 add a normalization factor
/*! returns L_N^alpha(x)*exp(-x/2)*norm
*/
double value(double x, int N) const {
r_8 value(r_8 x, int N, r_8 norm = 1.0) const {
using namespace std;
if (N<0) return 0;
double log2expfac = -0.5*x*inv_ln2;
r_8 log2expfac = -0.5*x*inv_ln2;
int scale=int(log2expfac/large_exponent2);
log2expfac-=scale*large_exponent2;
double pnm1=0., pn=exp(ln2*log2expfac);
double c1=alpha_-1., c2=alpha_-1.;
r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
r_8 c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
r_8 pnm2 = pnm1;
pnm1 = pn;
c1+=2.;
c2+=1.;
......@@ -133,23 +138,24 @@ protected:
return (scale<minscale) ? 0 : pn*cf[scale-minscale];
}//value
/*! returns L_n^alpha(x)*exp(-x/2) [n=0..N] in res.
//JEC 12/1/16 add a normalization factor
/*! returns L_n^alpha(x)*exp(-x/2)*norm [n=0..N] in res.
*/
void values (double x, int N, std::vector<double> &res) {
void values (r_8 x, int N, std::vector<r_8> &res, r_8 norm=1.0) {
using namespace std;
if (N<0) { res.resize(0); return; }
res.resize(N+1);
res[0]=exp(-0.5*x);
res[0]=exp(-0.5*x)*norm;
double log2expfac = -0.5*x*inv_ln2;
r_8 log2expfac = -0.5*x*inv_ln2;
int scale=int(log2expfac/large_exponent2);
log2expfac-=scale*large_exponent2;
double pnm1=0., pn=exp(ln2*log2expfac);
double c1=alpha_-1., c2=alpha_-1.;
r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
r_8 c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
r_8 pnm2 = pnm1;
pnm1 = pn;
c1+=2.;
c2+=1.;
......@@ -163,18 +169,26 @@ protected:
/*! computes nodes and weights for Gauss-Laguerre integration of order N
NOTE: the weights are correct for modified Laguerre polynomials, i.e.
L_N^alpha(x)*exp(-x/2)
The nodes/weights are indexed i:0,...,N-1
*/
void nodes_and_weights(int N, std::vector<double> &nodes,
std::vector<double> &weights) const {
void nodes_and_weights(int N, std::vector<r_8> &nodes,
std::vector<r_8> &weights) const {
using namespace std;
nodes.resize(N);
weights.resize(N);
double tmp1 = N+alpha_;
for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
tmp1 /= (N+1.);
/* JEC 12/1/16 alpha is r_8. tmp1 = (N+alpha)!/N!/(N+1)^2 = Gamma(N+alpha+1)/Gamma(N+1)/(N+1)^2
r_8 tmp1 = N+alpha_;
for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
tmp1 /= (N+1.);
Use the ratio of Gamma functions as in Boost library: Todo see if one should tune the Policy
*/
r_8 tmp1 = boost::math::tgamma_ratio((r_8)N+alpha_+1.0, r_8(N)+1.0);
tmp1 /= (N+1.)*(N+1.);
double x=0.;
r_8 x=0.;
// Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
for(int k=0; k<N; k++) {
if (k==0)
......@@ -184,8 +198,8 @@ protected:
else
{
int km1 = k-1;
double c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
double c1 = 1.+0.3*alpha_;
r_8 c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
r_8 c1 = 1.+0.3*alpha_;
x += c0*(x-nodes[k-2])/c1;
}
......@@ -194,7 +208,7 @@ protected:
bool done=false;
while (true)
{
double xold=x;
r_8 xold=x;
x+=newtoncorr (x, N);
if (done) break;
......@@ -203,7 +217,7 @@ protected:
nodes[k] = x;
double tmp0 = value(x, N+1);
r_8 tmp0 = value(x, N+1);
weights[k] = x*tmp1/(tmp0*tmp0);
}//end for
......@@ -221,12 +235,17 @@ protected:
sumW += wip;
// sumW += weights[i]*exp(-nodes[i]/tau);
}
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
/* JEC 12/1/16 alpha ! = alphaFact = Gamma(alpha+1)
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
r_16 sumRTh = (r_16)(N*(N+alpha_));
r_16 sumWTh = (r_16)(alphaFact);
*/
r_16 sumRTh = (r_16)(N*(N+alpha_));
r_16 sumWTh = (r_16)(alphaFact);
r_16 sumWTh = boost::math::tgamma((r_16)(alpha_+1.0));
cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
......@@ -252,13 +271,13 @@ protected:
protected:
int N_; //!< Order of the polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
r_8 alpha_; //!< second parameter of the generalized Laguerre polynomial
private:
enum { large_exponent2=90, minscale=-14, maxscale=14 };
double ln2, inv_ln2, fsmall, fbig;
std::vector<double> xi;
std::vector<double> cf;
r_8 ln2, inv_ln2, fsmall, fbig;
std::vector<r_8> xi;
std::vector<r_8> cf;
}; //LaguerreFuncQuad
......
......@@ -21,7 +21,8 @@ using namespace std;
namespace LagSHT {
LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha, bool R2tau) :
// LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha, bool R2tau) :
LaguerreTransform::LaguerreTransform(int N, r_8 R, r_8 alpha, bool R2tau) :
BaseLaguerreTransform(N,alpha), R_(R) {
#if DEBUG >= 1
......@@ -29,16 +30,18 @@ namespace LagSHT {
#endif
if(alpha<0) throw LagSHTError("LagTransform call with alpha<0");
LaguerreFuncQuad lag(N_,alpha_);
//reminder: LaguerreFuncQuad is called with N to compute nodes & weioghts order N
LaguerreFuncQuad lag(N_, alpha_);
lag.QuadWeightNodes(nodes_,weights_);
int tmp = 1;
for(int i=1;i<=alpha_; i++) tmp *= i;
alphaFact_ = sqrt((r_8)tmp); //sqrt(alpha!)
//JEC 12/1/16 alpha is a r_8
// int tmp = 1;
// for(int i=1;i<=alpha_; i++) tmp *= i;
// alphaFact_ = sqrt((r_8)tmp); //sqrt(alpha!)
alphaFact_ = sqrt(boost::math::tgamma(alpha_+1.0)); //<---------- usage of alpha for normalieation
//JEC 22/9/15
if(R2tau) {
tau_ = R_/nodes_[N_-1];
r_8 pow3inv2 = 3./2.;
......@@ -74,28 +77,29 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
#if DEBUG >= 2
std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
r_8 invalphaFact = 1.0/alphaFact_; //1/sqrt(alpha)!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
//JEC 12/1/16 still valid for alpha r_8 (Todo: put facts[n] in Ctor)
r_8 invalphaFact = 1.0/alphaFact_; //1/sqrt(alpha)!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) the normalization
facts[0] = invalphaFact;
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1);
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16; N-1 => Values() computes with n in [0,N-1]
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) k,n:0...N-1
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
lag.Values(rk,LnAll,prk);
for (int n = 0; n<N_; n++ ){
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_; //JEC 22/9/15 add tau-fact
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_;
}
}
#ifdef CBLAS
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fn[0]), 2*stride);
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fn[0]), 2*stride);
#else
vector<complex<r_8> > vtmp(N_);
for (int l=0; l<stride; l++) {
......@@ -129,28 +133,31 @@ void LaguerreTransform::MultiAnalysisR(const vector<r_8>& fi,
#if DEBUG >= 2
std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
//JEC 12/1/16 still valid for alpha r_8
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) for normalization
facts[0] = invalphaFact;
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1);
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
lag.Values(rk,LnAll,prk);
for (int n = 0; n<N_; n++ ){
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_; //JEC 22/9/15 add tau-fact
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_;
}
}
//#ifdef CBLAS
// cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fn[0]), 2*stride);
// cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fn[0]), 2*stride);
// throw LagSHTError("CBLAS to be adapted for REAL vectors");
//#else
vector<r_8> vtmp(N_);
......@@ -187,26 +194,29 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
#if DEBUG >= 2
std::copy(fn.begin(), fn.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
//JEC 12/1/16 still valid for alpha r_8 (Todo: put facts[n] in Ctor)
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
facts[0] = invalphaFact;
for(int n=1; n<N_; n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1);
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) n,:0...N-1
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
lag.Values(rk,LnAll);
r_8 rk = nodes_[k];
r_8 prk = pow(rk, alphap_); //JEC 12/1/16
lag.Values(rk,LnAll,prk); //<---- use the scaling
for (int n = 0; n<N_; n++ ){
LnkMtx[n*N_+k] = LnAll[n]*facts[n]*tau3sqrtinv_; //JEC 22/9/15 add tau-fact
LnkMtx[n*N_+k] = LnAll[n]*facts[n]*tau3sqrtinv_; //JEC 12/1/16 add prk
}
}
#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);
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fi[0]), 2*stride);
#else
cout << "CBLAS desactivated................" << endl;
vector<complex<r_8> > vtmp(N_);
......
......@@ -30,20 +30,28 @@
#include "laguerreBuilder.h"
/* JEC 12/1/16
1) alpha is r_8 even if here it will be called by user as an integer 0, 2
*/
namespace LagSHT {
class BaseLaguerreTransform {
public:
//! Creator
BaseLaguerreTransform(int n, int alpha=2): N_(n), alpha_(alpha) {}
// BaseLaguerreTransform(int n, int alpha=2): N_(n), alpha_(alpha) {}
BaseLaguerreTransform(int n, r_8 alpha=2): N_(n), alpha_(alpha) { alphap_ = alpha/2.0-1.0; } //JEC 13/1/16
//! Destructor
virtual ~BaseLaguerreTransform() {}
//! Getters
int GetOrder() { return N_; }
int GetAlpha() { return alpha_;}
// int GetAlpha() { return alpha_;}
r_8 GetAlpha() { return alpha_;}
protected:
int N_; //!< Order of the Generalized Lagurerre polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
// int alpha_; //!< second parameter of the generalized Laguerre polynomial
r_8 alpha_; //!< second parameter of the generalized Laguerre polynomial
r_8 alphap_; // auxiliary parameter
};//BaseLaguerreTransform
......@@ -51,14 +59,15 @@ protected:
///////////////////////////////////////////////////////////////////////
//// ------------------- Class LaguerreTransform ----------------- //
//// Perform single or multiple 1D-Laguerre Transform and Inverse //
//// using generalized Laguerre Polynomials except for T:double //
//// using generalized Laguerre Polynomials except for T:r_8 //
//// where in this case the generalized Laguerre functions are used //
///////////////////////////////////////////////////////////////////////
//! Laguerre Transform Syntheis/Analysis
class LaguerreTransform : public BaseLaguerreTransform {
public:
//! Creator
LaguerreTransform(int N, r_8 R, int alpha=2, bool R2tau=true);
// LaguerreTransform(int N, r_8 R, int alpha=2, bool R2tau=true);
LaguerreTransform(int N, r_8 R, r_8 alpha=2, bool R2tau=true);
//! Destructor
virtual ~LaguerreTransform() {}
......
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