Commit 1b065393 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC/MR) 24/4/15 restructuration of the 2D Geometry interface; introduction of...

(JEC/MR) 24/4/15 restructuration of the 2D Geometry interface; introduction of direct method to access 2D pixel location and index which remove dependance from ANN
parent e8a57ff7
......@@ -120,6 +120,18 @@ $(OBJ)lagsht_testsuite.o: lagsht_testsuite.cc $(CXXHDR)
echo "compile... $<"
$(CXXCOMPILE) $< -o $@
######################
.PHONY: fullcheck
fullcheck :
./Objs/lagsht_testsuite -t 0
./Objs/lagsht_testsuite -t 1 -n 1024
./Objs/lagsht_testsuite -t 2 -n 128
./Objs/lagsht_testsuite -t 3 -l 512 -g ECP
./Objs/lagsht_testsuite -t 3 -l 512 -g Gauss
./Objs/lagsht_testsuite -t 3 -l 512 -g Healpix
./Objs/lagsht_testsuite -t 4 -l 1024 -n 128 -nphi 2048 -g ECP
./Objs/lagsht_testsuite -t 3 -l 1024 -n 128 -nphi 2048 -g Gauss
./Objs/lagsht_testsuite -t 3 -l 1024 -n 128 -ntheta 512 -g Healpix
......
......@@ -5,40 +5,111 @@
#include "lagSphericTransform.h"
#include "lagsht_exceptions.h"
#include "lagsht_healpixhelper.h"
#include "walltimer.h" //timing
#define DEBUG 0
namespace LagSHT {
void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
//libsharp defini Lmax = L-1 avec une somme l:0...Lmax inclus
//------------------------- ECP/Fejer1 Geom ---------------------
int ECPGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
int stride_phi=1; //this migth be also a memeber data
int stride_theta=Nphi_; //"
r_8 phi0 = 0.; //for the time beeing
int itheta=int(Nrings_*theta/M_PI);
int iphi = int(Nphi_+ Nphi_*(phi-phi0)/(2*M_PI)+0.5)%Nphi_;
return stride_phi*iphi+stride_theta*itheta;
}//PixIndexSph
void ECPGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
//JEC if is not the case change the code int stride_phi=1; //this migth be also a memeber data
//JEC if is not the case change the code int stride_theta=Nphi_; //"
r_8 phi0 = 0.;
int iphi=k%Nphi_;
int itheta=k/Nphi_;
phi=phi0 + iphi*2*M_PI/Nphi_;
theta=(itheta+0.5)*M_PI/Nrings_;
//et pour les Alm
//set_trinagular_alm (Lmax,Mmax).
if(geomDone_){
cout << "-------------------------------------------------" <<endl;
cout << "----------------- New Geometry ------------- " <<endl;
cout << "-------------------------------------------------" <<endl;
geomDone_ = false;
}
if (choice == "ECP" || choice == "Fejer1") {
}//PixThetaPhi
void ECPGeometry::SetMap(int nrings, int nphi) {
geomType_ = "ECP";
if(geomDone_){
cout << "-------------------------------------------------" <<endl;
cout << "--------------- New ECP Geometry ------------- " <<endl;
cout << "-------------------------------------------------" <<endl;
geomDone_ = false;
}
//dans sharp_testsuite.c Fejer's 1st rule
//nrings = 2Lmax+1
//npixper ring = 2Mmax+1
Nrings_ = (nrings == -1) ? 2*L_-1 : nrings;
Nphi_ = (nphi == -1) ? 2*M_-1 : nphi;
sharpjob_.set_ECP_geometry(Nrings_,Nphi_);
//do the common features
SetCommonMapFeatures();
geomDone_=true;
}//SetMap
//------------------------- Gauss Geom ---------------------
int GaussGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
int stride_phi=1; //this migth be also a memeber data
int stride_theta=Nphi_; //"
r_8 phi0 = 0.; //for the time beeing
int itheta = std::min_element(thetaPix_.begin(),thetaPix_.end(), closer_to(theta)) - thetaPix_.begin();
int iphi = int(Nphi_+ Nphi_*(phi-phi0)/(2*M_PI)+0.5)%Nphi_;
return stride_phi*iphi+stride_theta*itheta;
}//PixIndexSph
void GaussGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
//Nb this is the same code as ECPGeometry
//JEC if is not the case change the code int stride_phi=1; //this migth be also a memeber data
//JEC if is not the case change the code int stride_theta=Nphi_; //"
r_8 phi0 = 0.;
int iphi=k%Nphi_;
int itheta=k/Nphi_;
phi=phi0 + iphi*2*M_PI/Nphi_;
theta=(itheta+0.5)*M_PI/Nrings_;
}//PixThetaPhi
void GaussGeometry::SetMap(int nrings, int nphi) {
geomType_ = "Gauss";
if(geomDone_){
cout << "-------------------------------------------------" <<endl;
cout << "--------------- New Gauss Geometry ------------- " <<endl;
cout << "-------------------------------------------------" <<endl;
geomDone_ = false;
}
} else if (choice == "Gauss") {
//dans sharp_testsuite.c Gauss
//nrings = Lmax+1
//npixper ring = 2Mmax+1
......@@ -47,18 +118,63 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
Nphi_ = (nphi == -1) ? 2*M_-1 : nphi;
sharpjob_.set_Gauss_geometry(Nrings_,Nphi_);
//do the common features
SetCommonMapFeatures();
geomDone_=true;
}//SetMap
//------------------------- HEALPIX Geom ---------------------
int HealpixGeometry::PixIndexSph(r_8 theta, r_8 phi) const {
if(theta<=0 || theta>=M_PI) LagSHTError("BaseLagSphTransform::PixIndexSph: Healpix theta out of range");
return ang2pix_ring_z_phi(Nrings_,cos(theta),phi);
}//PixIndexSph
void HealpixGeometry::PixThetaPhi(int k, r_8& theta, r_8& phi) const {
r_8 z;
pix2ang_ring_z_phi (Nrings_,k,&z,&phi);
theta = acos(z);
}//PixThetaPhi
void HealpixGeometry::SetMap(int nsides, int DUMMY) { //second argument not used
geomType_ = "Healpix";
if(geomDone_){
cout << "-------------------------------------------------" <<endl;
cout << "--------------- New Healpix Geometry --------- " <<endl;
cout << "-------------------------------------------------" <<endl;
geomDone_ = false;
}
} else if (choice == "Healpix") {
//dans sharp_testsuite.c Healpix
//nside = Lmax/2;
Nrings_ = (nrings == -1 ) ? (L_ -1)/2: nrings; // this is the nside argument in fact
Nrings_ = (nsides == -1 ) ? (L_ -1)/2: nsides; // this is the nside argument in fact
Nphi_ = -1; // is not used
sharpjob_.set_Healpix_geometry(Nrings_);
} else {
throw LagSHTError("BaseLagSphTransform::SetThetaPhiMap: Unknown geometry <"+choice+">");
}
//do the common features
SetCommonMapFeatures();
geomDone_=true;
}//SetMap
//---------------------------------------
void BaseGeometry::SetCommonMapFeatures() {
sharpjob_.set_triangular_alm_info(L_-1,M_-1);
......@@ -69,11 +185,9 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
//fill the ANN-kdTree to prepare the filling of the 2D spheres
// geom_t vecPix;
SetPixelPositions();
ann.SetkdTree(vecPix_);
geomDone_ = true;
cout << "SetThetaPhiMap: Geom, Nring, Nphi, Npix, Nalm: " << choice << ", "
cout << "SetThetaPhiMap: Geom, Nring, Nphi, Npix, Nalm: " << geomType_ << ", "
<< Nrings_ << ", "
<< Nphi_ << ", "
<< Npix_ << ", "
......@@ -82,11 +196,16 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
}//SetThetaPhiMap
void BaseLagSphTransform::SetPixelPositions() {
//----------------------------------------------------------
void BaseGeometry::SetPixelPositions() {
const sharp_geom_info *ginfo;
ginfo = sharpjob_.get_geom_info();
for (int i=0; i<ginfo->npairs; ++i) {
double theta1 = ginfo->pair[i].r1.theta;
thetaPix_.push_back(theta1);
// and so on for nph, phi0
int nph1 = ginfo->pair[i].r1.nph;
double phi0 = ginfo->pair[i].r1.phi0;
......@@ -97,6 +216,7 @@ void BaseLagSphTransform::SetPixelPositions() {
}
if (ginfo->pair[i].r2.nph>0) {// the second ring in this pair exists
double theta2 = ginfo->pair[i].r2.theta;
thetaPix_.push_back(theta2);
int nph2 = ginfo->pair[i].r2.nph;
double phi0 = ginfo->pair[i].r2.phi0;
for(int j=0; j<nph2; j++) {
......@@ -107,12 +227,12 @@ void BaseLagSphTransform::SetPixelPositions() {
}
}
sort(vecPix_.begin(),vecPix_.end()); //See if it is necessary...
sort(thetaPix_.begin(), thetaPix_.end()); //mandatory for PixIndexSph
sort(vecPix_.begin(),vecPix_.end());
}//SetPixelPositions
//------------------------------------------------
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk) {
......@@ -152,7 +272,8 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
#endif
// sharpjob_.alm2map(&(flmk.data()[off7n].real()),&fijk.data()[off7k],false); //the false is to forbidd accumulation
sharpjob_.alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
// sharpjob_.alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
sphere_->alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
}//loop on shell
......@@ -194,7 +315,9 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
// sharpjob_.map2alm(&(fijk.data()[off7k]),&(flmk.data()[off7n].real()),false); //the false is to forbidd accumulation
sharpjob_.map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(flmk.data()[off7n])),false); //the false is to forbidd accumulation
// sharpjob_.map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(flmk.data()[off7n])),false); //the false is to forbidd accumulation
sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(flmk.data()[off7n])),false); //the false is to forbidd accumulation
}//loop on shell
......
......@@ -12,7 +12,7 @@ using namespace std;
#include "laguerreTransform.h"
#include "lagsht_annserver.h"
#include "lagsht_geom.h"
#include "lagsht_exceptions.h"
......@@ -24,63 +24,80 @@ namespace LagSHT {
////////////////////////////////////////////////////////////////
class BaseLagSphTransform {
public:
//Todo: See if Nmax and alpha can be transfered from BaseGeomerty.... to LaguerreSphericalTransform
class BaseGeometry {
public:
//! Creator
BaseLagSphTransform(int Lmax, int Mmax, int Nmax, int alpha=2):
L_(Lmax), M_(Mmax), N_(Nmax), alpha_(alpha),
Nrings_(0), Nphi_(0), Npix_(0), Nalm_(0),geomDone_(false) {
BaseGeometry(int Lmax, int Mmax, int Nmax, int Nrings, int Nphi, int alpha=2):
L_(Lmax), M_(Mmax), N_(Nmax),
Nrings_(Nrings), Nphi_(Nphi), alpha_(alpha),
Npix_(0), Nalm_(0), geomType_(""), geomDone_(false) {
if (Mmax<0) M_=Lmax;
}
//! Destructor
virtual ~BaseLagSphTransform() {}
virtual ~BaseGeometry() {}
//! Number of colatitude rings of the 2D-sphere representation
int Nrings() { return Nrings_; }
int Nrings() const { return Nrings_; }
//! Side parameter for Healpix
int Nsides() const { return Nrings_; } //this is a convention here
//! Number of colatitude rings of the 2D-sphere representation
int NTheta() { return Nrings(); }
int NTheta() const { return Nrings(); }
//! Number of phi values (pixels) per colatitude rings
int NPhi() { return Nphi_; }
int NPhi() const { return Nphi_; }
//! Total number of 2D pixels
int Npix() { return Npix_; }
int Npix() const { return Npix_; }
//! Total number of Alm coefficients
int Nalm() { return Nalm_; }
int Nalm() const { return Nalm_; }
//! Type of geomery
string Geometry() const { return geomType_; }
//! Sorted List of pixel center (theta, phi) positions
geom_t GetPixelPositions() { return vecPix_; }
/*! Choose the 2D map representation & initialization of libsharp
\input ECP/Fejer1, Gauss, and Healpix
\input nrings
\input nphi
*/
void SetThetaPhiMap(string choice="Fejer1", int nrings=-1, int nphi=-1);
geom_t GetPixelPositions() const { return vecPix_; }
/*! Get Pixel index in the 2D map
\input theta in [0, pi]
\input phi in [0, 2pi]
*/
int PixIndexSph(r_8 theta, r_8 phi){ return ann.GetClosestId(theta, phi); }
virtual int PixIndexSph(r_8 theta, r_8 phi) const = 0;
/*! Get Theta Phi position of a given pixel index
\input k in [0, Npix-1] no boundchckeing done
*/
virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const = 0;
/*! Choose the 2D map representation & initialization of libsharp
\input ECP/Fejer1, Gauss, and Healpix
\input nrings
\input nphi
*/
virtual void SetMap(int nrings=-1, int nphi=-1) = 0;
void PixThetaPhi(int k, r_8& theta, r_8& phi) {
pixpos_t pix = vecPix_[k];
theta = pix.first;
phi = pix.second;
virtual void alm2map(const r_8 *alm, r_8 *map, bool add=false){
sharpjob_.alm2map(alm,map,add);
}
virtual void map2alm(const r_8 *map, r_8 *alm, bool add=false){
sharpjob_.map2alm(map,alm,add);
}
protected:
/*! Set the vector of pixel positions
*/
void SetPixelPositions();
/*! Manage the commonality of the different geometries
*/
void SetCommonMapFeatures();
//from sharp test suite
static ptrdiff_t get_npix(const sharp_geom_info *ginfo) {
ptrdiff_t res=0;
......@@ -99,39 +116,199 @@ protected:
return res;
}
protected:
int L_; //!< Spherical harmonic limit such that l:0,...,L-1 m:0,...,l
int M_; //!< maximum value of m
int N_; //!< Radial band limit such that n:0,...,N-1
int alpha_; //!< parameter of the generalized Laguerre Polynomials
int Nrings_; //!< Number of colatitude rings of the 2D-sphere representation
int Nphi_; //!< Number of phi values (pixels) per colatitude rings
int alpha_; //!< parameter of the generalized Laguerre Polynomials
int Npix_; //!< Total number of 2D pixels
int Nalm_; //!< Total number of Alm coefficients
geom_t vecPix_; //!< vector of pixel positions
geom_t vecPix_; //!< sorted vector of pixel positions
vector<r_8> thetaPix_; //!< sorted vector of theta rings
string geomType_; //!<< Type of pixelization
bool geomDone_; //!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 extend the validity of the computation
ANNServer ann; //!< utility to fill the spheric shells
}; //BaseLagSphTransform
}; //BaseGeometry
//-----------------------------------------------------
class ECPGeometry : public BaseGeometry {
public:
//! Constructor
ECPGeometry(int Lmax, int Mmax, int Nmax, int Nrings, int Nphi, int alpha=2):
BaseGeometry(Lmax, Mmax, Nmax, Nrings, Nphi, alpha) {}
//! Destructor
virtual ~ECPGeometry() {}
/*! Get Pixel index in the 2D map
\input theta in [0, pi]
\input phi in [0, 2pi]
*/
virtual int PixIndexSph(r_8 theta, r_8 phi) const;
/*! Get Theta Phi position of a given pixel index
\input k in [0, Npix-1] no boundchckeing done
*/
virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
/*! Choose the 2D map representation & initialization of libsharp
\input ECP/Fejer1
\input nrings
\input nphi
*/
virtual void SetMap(int nrings=-1, int nphi=-1);
}; //ECPGeometry
//-----------------------------------------------------
//alias of different geometry
#define Fejer1Geometry ECPGeometry
//-----------------------------------------------------
class GaussGeometry : public BaseGeometry {
public:
//! Constructor
GaussGeometry(int Lmax, int Mmax, int Nmax, int Nrings, int Nphi, int alpha=2):
BaseGeometry(Lmax, Mmax, Nmax, Nrings, Nphi, alpha) {}
//! Destructor
virtual ~GaussGeometry() {}
/*! Get Pixel index in the 2D map
\input theta in [0, pi]
\input phi in [0, 2pi]
*/
virtual int PixIndexSph(r_8 theta, r_8 phi) const;
/*! Get Theta Phi position of a given pixel index
\input k in [0, Npix-1] no boundchckeing done
*/
virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
/*! Choose the 2D map representation & initialization of libsharp
\input Gauss
\input nrings
\input nphi
*/
virtual void SetMap(int nrings=-1, int nphi=-1);
protected:
/*! Utility function object to perform PixIndexSph theta search
*/
struct closer_to {
r_8 target_;
explicit closer_to( r_8 target ) : target_( target ){}
class LaguerreSphericalTransform : public BaseLagSphTransform {
bool operator ()( r_8 left, r_8 right ) const
{
return std::abs( target_ - left ) < std::abs( target_ - right );
}
};
}; //GaussGeometry
//-----------------------------------------------------
class HealpixGeometry : public BaseGeometry {
public:
/*! Constructor
Nrigns = Nsides here
*/
HealpixGeometry(int Lmax, int Mmax, int Nmax, int Nsides, int alpha=2):
BaseGeometry(Lmax, Mmax, Nmax, Nsides,-1,alpha) {}
//! Destructor
virtual ~HealpixGeometry() {}
/*! Get Pixel index in the 2D map
\input theta in [0, pi]
\input phi in [0, 2pi]
*/
virtual int PixIndexSph(r_8 theta, r_8 phi) const;
/*! Get Theta Phi position of a given pixel index
\input k in [0, Npix-1] no boundchckeing done
*/
virtual void PixThetaPhi(int k, r_8& theta, r_8& phi) const;
/*! Choose the 2D map representation & initialization of libsharp
\input Healpix
\input nsides
*/
virtual void SetMap(int nsides=-1, int DUMMY=-1); //2nd arg not used but here for compatibility
}; //HealpixGeometry
//-----------------------------------------------------
class LaguerreSphericalTransform {
public:
//! Constructor
LaguerreSphericalTransform(int Lmax, int Mmax, int Nmax, r_8 R, int alpha=2):
BaseLagSphTransform(Lmax,Mmax,Nmax,alpha), R_(R), lagTrans_(Nmax,(r_8)R,alpha) {
LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2):
R_(R), lagTrans_(Nmax, R, alpha) {
//Factory
if (spheregeom == "ECP" || spheregeom == "Fejer1") {
sphere_ = new ECPGeometry(Lmax, Mmax, Nmax, Nrings, Nphi, alpha);
} else if (spheregeom == "Gauss") {
sphere_ = new GaussGeometry(Lmax, Mmax, Nmax, Nrings, Nphi, alpha);
} else if (spheregeom == "Healpix") {
sphere_ = new HealpixGeometry(Lmax, Mmax, Nmax, Nrings, alpha); //here Nsides = Nrings
} else {
cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << endl;
sphere_ = 0;
throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
}
//do the map pixelization
if(sphere_) {
sphere_->SetMap(Nrings, Nphi);
//store localy usefull parameters for the algorithms Analysis/Synthesis
Npix_ = sphere_->Npix();
Nalm_ = sphere_->Nalm();
N_ = lagTrans_.GetOrder();
alpha_= lagTrans_.GetAlpha();
}
}//Ctor
//! Destructor
virtual ~LaguerreSphericalTransform() {}
virtual ~LaguerreSphericalTransform() {
if(sphere_) delete sphere_;
sphere_ = 0;
}
//! Geometry helper
BaseGeometry* GetSphericalGeometry() const { return sphere_;}
//! Getters
int GetOrder() { return N_; }
int GetAlpha() { return alpha_;}
//! Laguerre Transform helper
//! Synthesis
/*! \brief Coeffs -> Pixels with Input/Output using T floating representation
......@@ -150,8 +327,17 @@ class LaguerreSphericalTransform : public BaseLagSphTransform {
protected:
r_8 R_; //!< Radial maximal value (not yet used)
BaseGeometry* sphere_; //<! the 2D pixelization of the shells
LaguerreTransform lagTrans_; //!< the Laguerre Transform
int Npix_; //!< Total number of 2D pixels
int Nalm_; //!< Total number of Alm coefficients