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

(JEC) 25/4/15 decouple 2D-spherical-geometry management from the Spherical Laguerre Transform part

parent f23ef8b6
......@@ -63,11 +63,13 @@ clean :
#C++ common Objects
CXXOBJ = $(OBJ)laguerreBuilder.o \
$(OBJ)laguerreTransform.o \
$(OBJ)lagsht_spheregeom.o \
$(OBJ)lagSphericTransform.o \
$(OBJ)walltimer.o
CXXSHOBJ = laguerreBuilder.o \
laguerreTransform.o \
lagsht_spheregeom.o \
lagSphericTransform.o \
walltimer.o
......@@ -77,6 +79,7 @@ CXXHDR = lagsht_exceptions.h \
lagsht_numbers.h \
lagsht_utils.h \
lagsht_geom.h \
lagsht_spheregeom.h \
lagsht_healpixhelper.h \
laguerreBuilder.h \
laguerreTransform.h \
......@@ -130,8 +133,8 @@ fullcheck :
./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
./Objs/lagsht_testsuite -t 4 -l 1024 -n 128 -nphi 2048 -g Gauss
./Objs/lagsht_testsuite -t 4 -l 1024 -n 128 -ntheta 512 -g Healpix
......
......@@ -12,8 +12,8 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git
* lagsht_execptions.h: class derived from std::exceptions.
* walltimer.h (.cc): utilities to profile the program.
* lagsht_utils.h: code by David Robert Nadeau to get the avalable memory
* lagsht_annserver.h : access to the ANN tool (see below)
* lagsht_geom.h : define some typedef used in the code
* lagsht_spheregeom.h (.cc) : set of classes to deeal with the 2D sphere pixelization and libsharp initialization
* lagsht_geom.h : define some typedef used in the code for geometry
* lagsht_testsuite.c : a simple program to test different piece of the code
* Makefile that should be tuned to the local platform (*.inc files)
* doxydoc / footer.html : input file to doxygen tool to generate the class documentation
......@@ -34,12 +34,6 @@ If the variable `CBLAS` is defined then the code of laguerreTransform.cc use the
`git clone https://github.com/xianyi/OpenBLAS.git`
**WARNING**: on Darwin the Makefile use the `Accelerate Framework` features which are a bit more effective than OpenBLAS.
* ANN (under test) :
get the code here: `http://www.cs.umd.edu/~mount/ANN`
ANN is a library written in C++ which supports data structures and algorithms for both exact and approximate nearest neighbor searching in arbitrarily high dimensions. It is distributed under the terms of the GNU Lesser Public License. The code require an ANSII C++ and is tested under different platforms. I have build it using g++ compilor on Mac OS X which was not mentionned by the author. ANN is used to perform to find the closest pixel index given the (theta, phi) position of a point.
# Compilation
> edit Makefile and adapat to local platform :
......@@ -51,10 +45,15 @@ ANN is a library written in C++ which supports data structures and algorithms fo
> make MACH=Linux BLAS=1 : on Linux with BLAS ON
> make : (default) on Darwin and use Accelerate Framework
```
The result of "make" is binanry file under ./Objs directory as well as `lagsht.a` library
> make check
runs a computation of the Nodes & Weights for N=1024
and compare to the xxx.txt.TEST files to produce a
Numerical Error estimate.
> make fullcheck
runs a series of `lagsht_testsuite` which activates different parts of the code and shows how-to call this test program.
# Plateform tested
Mac OS X 10.9.5 + gcc 4.8.4
Linux SLC 6.6 + gcc 4.9.1 20140922
......@@ -70,22 +69,28 @@ ANN is a library written in C++ which supports data structures and algorithms fo
0: basic test for numerical application
1: to get the nodes & weights of the generalized Gauss-laguerre function quadrature.
2: to test 1D- Laguerre Transform (Synthesis & Analysis)
3: to test the full 3D Spherical Harmonic Laguerre trabsform (Synthesis & Analysis)
3: to test 2D- pixelization schemes from libsharp
4: to test the full 3D Spherical Harmonic Laguerre trabsform (Synthesis & Analysis)
> <geometry>
Gauss : use the Gauss quadrature by default Ntheta = Lmax and Nphi = 2Lmax - 1
Fejer1: use the Fejer's sum rule n°1 by default Ntheta = 2Lmax-1 and Nphi = 2Lmax-1
Healpix: use the Healpix map with Ntheta giving the Nside parameter.
Gauss : use the Gauss quadrature by default Ntheta = Lmax and
Nphi = 2Lmax - 1.
Fejer1: use the Fejer's sum rule n°1 by default
Ntheta = 2Lmax-1 and Nphi = 2Lmax-1.
Healpix: use the Healpix map with Ntheta giving the Nsides
parameter.
> nphi: libsharp works faster if this number is increased (compared to the minimal number required by the geometry) up to a power of 2 number. Eg. for Lmax = 1024, Nphi_min = 2L-1 = 2047, so Nphi = 2048 will run faster.
# Examples
* ./Objs/lagsht_testsuite -n 1024
gives 2 files `lagNodes-1024-Func.txt` and `lagWeights-1024-Func.txt` with the Nodes & Weights computed
* ./Objs/lagsht_testsuite -t 3 -l 1024 -n 128
* ./Objs/lagsht_testsuite -t 4 -l 1024 -n 128
perform a Laguerre SH transform with the default geometry (Gauss) using Lmax = 1024 and Nmax = 128. The program generates a set of flmn coefficients, then performs a Synthesis followed by an Analysis operations. Finally it computes the maximal absolute and relative errors comparing the original flmn complex values to the results of the Analysis. The timing of the different processes is shown in a sorted tree (here on Darwin Mac Os X with Accelerate Framework).
````
./Objs/lagsht_testsuite -t 3 -l 1024 -n 128
./Objs/lagsht_testsuite -t 4 -l 1024 -n 128
Max Memory size: 8589 MBytes
___________ MultiSphericalLaguerreTransform TEST _____________
LaguerreTransform start....
......
......@@ -5,7 +5,7 @@
#include "lagSphericTransform.h"
#include "lagsht_exceptions.h"
#include "lagsht_healpixhelper.h"
//#include "lagsht_healpixhelper.h"
#include "walltimer.h" //timing
......@@ -16,224 +16,6 @@
namespace LagSHT {
//------------------------- 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_;
}//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;
}
//dans sharp_testsuite.c Gauss
//nrings = Lmax+1
//npixper ring = 2Mmax+1
Nrings_ = (nrings == -1) ? L_ : nrings;
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;
}
//dans sharp_testsuite.c Healpix
//nside = Lmax/2;
Nrings_ = (nsides == -1 ) ? (L_ -1)/2: nsides; // this is the nside argument in fact
Nphi_ = -1; // is not used
sharpjob_.set_Healpix_geometry(Nrings_);
//do the common features
SetCommonMapFeatures();
geomDone_=true;
}//SetMap
//---------------------------------------
void BaseGeometry::SetCommonMapFeatures() {
sharpjob_.set_triangular_alm_info(L_-1,M_-1);
Npix_ = get_npix(sharpjob_.get_geom_info());
Nalm_ = get_nalms(sharpjob_.get_alm_info());
//fill the ANN-kdTree to prepare the filling of the 2D spheres
// geom_t vecPix;
SetPixelPositions();
cout << "SetThetaPhiMap: Geom, Nring, Nphi, Npix, Nalm: " << geomType_ << ", "
<< Nrings_ << ", "
<< Nphi_ << ", "
<< Npix_ << ", "
<< Nalm_
<< endl;
}//SetThetaPhiMap
//----------------------------------------------------------
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;
for(int j=0; j<nph1; j++) {
double phi1 = phi0 + j*(2.0*M_PI/nph1);
// cout << "vecPix_s ( "<< theta1 << ", " << phi1 << ")" << endl;
vecPix_.push_back(make_pair(theta1,phi1));
}
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++) {
double phi2 = phi0 + j*(2.0*M_PI/nph2);
// cout << "vecPix_s ( "<< theta2 << ", " << phi2 << ")" << endl;
vecPix_.push_back(make_pair(theta2,phi2));
}
}
}
sort(vecPix_.begin(),vecPix_.end()); //See if it is necessary...
sort(thetaPix_.begin(), thetaPix_.end()); //mandatory for PixIndexSph
}//SetPixelPositions
//------------------------------------------------
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk) {
int Ntot3DPix = Npix_*N_;
......@@ -270,10 +52,7 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
#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
sphere_->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
......@@ -313,10 +92,6 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
cout << "Start Harmonic Spheric... at n = "<< n << endl;
#endif
// 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
sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(flmk.data()[off7n])),false); //the false is to forbidd accumulation
......
#ifndef LAGSPHERICTRANSFORM_SEEN
#define LAGSPHERICTRANSFORM_SEEN
#include <utility> // std::pair, std::make_pair
using namespace std;
// #include <utility> // std::pair, std::make_pair
// using namespace std;
//Sharp
#include "sharp_lowlevel.h"
#include "sharp_geomhelpers.h"
#include "sharp_almhelpers.h"
#include "sharp_cxx.h"
// #include "sharp_lowlevel.h"
// #include "sharp_geomhelpers.h"
// #include "sharp_almhelpers.h"
// #include "sharp_cxx.h"
#include "lagsht_spheregeom.h"
#include "laguerreTransform.h"
#include "lagsht_geom.h"
//#include "lagsht_geom.h"
#include "lagsht_exceptions.h"
......@@ -23,242 +23,6 @@ namespace LagSHT {
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////
//Todo: See if Nmax and alpha can be transfered from BaseGeomerty.... to LaguerreSphericalTransform
class BaseGeometry {
public:
//! Creator
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 ~BaseGeometry() {}
//! Number of colatitude rings of the 2D-sphere representation
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() const { return Nrings(); }
//! Number of phi values (pixels) per colatitude rings
int NPhi() const { return Nphi_; }
//! Total number of 2D pixels
int Npix() const { return Npix_; }
//! Total number of Alm coefficients
int Nalm() const { return Nalm_; }
//! Type of geomery
string Geometry() const { return geomType_; }
//! Sorted List of pixel center (theta, phi) positions
geom_t GetPixelPositions() const { return vecPix_; }
/*! 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 = 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;
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;
for (int i=0; i<ginfo->npairs; ++i){
res += ginfo->pair[i].r1.nph;
if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
}
return res;
}
//from sharp test suite
static ptrdiff_t get_nalms(const sharp_alm_info *ainfo) {
ptrdiff_t res=0;
for (int i=0; i<ainfo->nm; ++i)
res += ainfo->lmax-ainfo->mval[i]+1;
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 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_; //!< 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
}; //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 ){}
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:
......@@ -269,12 +33,11 @@ class LaguerreSphericalTransform {