Skip to content

GitLab

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

(JEC) 14/10/15 add new directory test with 3D analysis exemples

parent 1fef3817
Hide whitespace changes
Inline Side-by-side
Showing with 822 additions and 0 deletions
test/Makefile 0 → 100644
View file @ 75d74b7c
ifndef MACH
MACH := $(shell uname)
endif
ifeq ($(MACH), Linux)
include ../Linux_g++_make.inc
ifeq ($(BLAS), 1)
# ===== OpenBLAS ========
BLASYES = -DCBLAS
BLASDIR = /exp/opera/kits/OpenBLAS
BLASINC = -I$(BLASDIR)/include
BLASLIB = -L$(BLASDIR)/lib
BLASLIBN = $(BLASLIB) -lopenblas
endif
else
include ../Darwin_g++_omp_make.inc
# ===== nNative BLAS (Darwin)
BLASYES = -DCBLAS -flax-vector-conversions
BLASLIBN = -framework Accelerate
ifeq ($(BLAS), 0)
BLASYES =
BLASLIBN =
endif
endif
# ===== Libsharp =====
#SHARPDIR = <yourdir>/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
#======= LagSHT ======
LAGSHTDIR = ..
LAGSHTLIB = $(LAGSHTDIR)/lib
LAGSHTINC = -I$(LAGSHTDIR)/include
LAGSHTLIBN = $(LAGSHTLIB)/liblagsht.a
SRC = ./src/
LIB = ./lib/
OBJ = ./objs/
EXE = ./bin/
DATA = ./data/
INCL = ./include/LagSHT/
# Define our target list
.PHONY: default
default: makedir test3D
.PHONY: all
all : makedir test3D
.PHONY: makedir
makedir :
@mkdir -p $(OBJ)
@mkdir -p $(EXE)
confcheck :
@echo "Config check mach=<$(MACH)> sharp=<$(SHARPDIR)> blas=<$(BLASDIR)>"
ifeq ($(BLAS), 1)
@if [ ! -d "$(BLASDIR)" ]; then \
echo "BLAS directory not defined properly!"; \
exit 1; \
fi ;
endif
@if [ ! -d "$(SHARPDIR)" ]; then \
echo "SHARP directory not defined properly!"; \
exit 1; \
fi ;
.PHONY: tidy
tidy :
find . -name "*~" | xargs -I {} rm {}
.PHONY: clean
clean :
rm -f $(OBJ)/*
rm -f $(EXE)/*
#C++ common Objects
CXXOBJ =
CXXSHOBJ =
#C++ common Headers
CXXHDR =
CPPFLAGS += $(BLASYES) $(SHARPINC) $(BLASINC) $(LAGSHTINC)
LDFLAGS += $(SHARPLIBN) $(BLASLIBN) $(LAGSHTLIBN) -lm
######################
.PHONY: test3D
test3D: $(EXE)test3D
echo '---- test3D made'
$(EXE)test3D : $(OBJ)test3D.o $(CXXOBJ)
echo "Link..."
$(CXXLINK) -o $@ $(OBJ)test3D.o $(CXXOBJ) $(LDFLAGS)
$(OBJ)test3D.o: test3D.cc $(CXXHDR)
echo "compile... $<"
$(CXXCOMPILE) $< -o $@
test/test3D-t3.txt 0 → 100644
View file @ 75d74b7c
Max Memory size: 8589 MBytes
Configuration parameters are set to: Test number : 3
Lmax, Nmax, Pmax, Rmax , alpha: 1024, 128, 0, 5, 2
Geometry: Gauss Ntheta, Nphi: -1, -1
Generate Cl(k) spectra and fill a 3D-ball
Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell
Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn
Error analysis on Cl(k) input and output of the Laguerre Spherical Transf.
>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<
Err. Max. 4.65397e-08 [0, 298], Err. Rel. 3.97506e-05
>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<
Total wall clock time for 'Test 3': 51.8014s
|
+- Generate fijk 3D-ball : 54.67% (28.3221s)
+- Lag Spherical Transf. : 42.01% (21.7609s)
| |
| +- SHT Analysis : 90.37% (19.6658s)
| +- Laguerre MultiAnalysis: 4.23% ( 0.9205s)
| +- <unaccounted> : 5.40% ( 1.1746s)
|
+- Clk and Cln OUT : 2.62% ( 1.3590s)
+- error analysis : 0.00% ( 0.0006s)
+- <unaccounted> : 0.69% ( 0.3588s)
Accumulated timing overhead: approx. 0.0000s
---/ Fin bloc try ----
---- Programme test3D.cc - FIN (Rc=0) ---
test/test3D-t5.txt 0 → 100644
View file @ 75d74b7c
Max Memory size: 8589 MBytes
Configuration parameters are set to: Test number : 5
Lmax, Nmax, Pmax, Rmax , alpha: 1024, 128, 0, 5, 2
Geometry: Gauss Ntheta, Nphi: -1, -1
Generate a 3D-Cube 1024^3
Fill the 3D-Ball from the Cube
Fill a 3D-Ball (2,2,2)Mpc^3: Rmax is set to : 819.2
Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell
Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn
Total wall clock time for 'Test 5': 274.0201s
|
+- Generate a 3D-Cube : 42.19% (115.6206s)
+- Fill the 3D-Ball from the Cube: 41.79% (114.4995s)
+- Lag Spherical Transf. : 14.48% ( 39.6700s)
| |
| +- SHT Analysis : 79.37% ( 31.4845s)
| +- Laguerre MultiAnalysis : 9.26% ( 3.6734s)
| +- <unaccounted> : 11.37% ( 4.5120s)
|
+- Clk and Cln OUT : 0.44% ( 1.2128s)
+- <unaccounted> : 1.10% ( 3.0172s)
Accumulated timing overhead: approx. 0.0010s
---/ Fin bloc try ----
---- Programme test3D.cc - FIN (Rc=0) ---
test/test3D.cc 0 → 100644
View file @ 75d74b7c
#include <string.h> //strcmp
#include <utility> // std::pair
//LagSHT
#include "LagSHT/lagsht_numbers.h"
#include "LagSHT/lagsht_exceptions.h"
#include "LagSHT/lagsht_utils.h" //getMemorySize
#include "LagSHT/walltimer.h"
#include "LagSHT/lagSphericTransform.h"
using namespace LagSHT;
template <typename T>
class Vector2d {
public:
Vector2d(size_t d1=0, size_t d2=0, T const & t=T()) :
d1(d1), d2(d2), data(d1*d2, t){ totsize = d1*d2; }
T & operator()(size_t i, size_t j) {
return data[i + j*d1];
}
T const & operator()(size_t i, size_t j, size_t k) const {
return data[i + j*d1];
}
size_t size() const { return totsize; }
private:
size_t d1,d2;
size_t totsize;
vector<T> data;
};
template <typename T>
class Vector3d {
public:
Vector3d(size_t d1=0, size_t d2=0, size_t d3=0, T const & t=T()) :
d1(d1), d2(d2), d3(d3), data(d1*d2*d3, t){ totsize = d1*d2*d3; }
T & operator()(size_t i, size_t j, size_t k) {
return data[i + j*d1 + k*d1*d2];
}
T const & operator()(size_t i, size_t j, size_t k) const {
return data[i + j*d1 + k*d1*d2];
}
size_t size() const { return totsize; }
private:
size_t d1,d2,d3;
size_t totsize;
vector<T> data;
};
//-------- 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;
string geometry;
int ntheta;
int nphi;
} param ;
//---------- Utility functions and class for this test program
namespace LagSHT {
//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
static double gaussrand(int* state,double mean=0., double sigma=1.0) {
r_8 rd1 = drand(0.,1.,state);
while (rd1 == 0.) rd1=drand(0.,1.,state);
r_8 rd2 = drand(0.,1.,state);
return (sqrt(-2.*log(rd1))*cos(2.*M_PI*rd2))*sigma + mean;
}
class TestTransform : public GaussGeometry {
public:
//! Constructor
TestTransform(int Lmax, int Mmax, int Nrings, int Nphi):
GaussGeometry(Lmax, Mmax, Nrings, Nphi) {
}//Ctor
//! Destructor
virtual ~TestTransform() {}
//! Forward map -> coeff
void Forward(const vector<r_8>& map, vector<complex<r_8> >& coeff){
// cout << "Forward map: " << map.size() << endl;
// cout << "coeff: " << coeff.size() << endl;
map2alm(&(map.data()[0]),reinterpret_cast<r_8*>(&coeff.data()[0]),false);
}
//! Backward coeff -> map
void Backward(const vector< complex<r_8> >& coeff, vector<r_8>& map) {
// cout << "Backward map: " << map.size() << endl;
// cout << "coeff: " << coeff.size() << endl;
alm2map(reinterpret_cast<const r_8*>(&coeff.data()[0]),&(map.data()[0]),false);
}
//
void ComputeCl(const vector<complex<r_8> >& alm, vector<r_8>& cl){
if(cl.size() != (size_t)L_) {
cout << "WARNING: resising Cl" << endl;
cl.resize(L_,0.);
}
for(int l=0;l<L_;l++) {
for (int m=0;m<=l;m++){
//sharp LambdaLM
size_t i = l+m*L_-m*(m+1)/2;
// cout << "(l,m): " << l << "," << m << " a-> " << alm[i] << endl;
r_8 alm2 = alm[i].real()*alm[i].real() + alm[i].imag()*alm[i].imag();
cl[l] += (m==0) ? alm2 : 2.0*alm2;
}
cl[l] /= (r_8)(2*l+1);
}
}//ComputeCl
void GenAlmFromCl(const vector<r_8>& cl, vector< complex<r_8> >& alm){
int state=1234567 + 8912 ; //random seed
for(int l=0;l<L_;l++) {
r_8 rms = sqrt(cl[l]);
for (int m=0;m<=l;m++){
//sharp LambdaLM
size_t i = l+m*L_-m*(m+1)/2;
if(m==0){
alm[i] = complex<r_8>(gaussrand(&state,0.,rms));
} else {
complex<r_8>tmp = complex<r_8>(gaussrand(&state,0.,1.),
gaussrand(&state,0.,1.));
alm[i] = ((r_8)(rms*M_SQRT1_2))*tmp;
}
}
}
}//GenAlmFromCl
};
}//end namespace
//------------------------------------
//test3:
// generate Clin(k) k here index of shells "a la" CMB Cl gaussian profile
// compute Alm(k) from Clin(k)
// compute a map(k) a from Inverse Spherical Harmonic Transform and fill fijk 3D-vector
// proceed to LagSHT Analysis transform to get flmn and AlmOut(k) coeff.
// compute Clout(k) from AlmOut(k) which should be identical to Clin(k)
// compute Clout(n) a la CMB by reducing the "m" index
//------------------------------------
void test3() {
//initial values
int Lmax = param.Lmax;
int Nmax = param.N;
r_8 Rmax = param.R;
string geom = param.geometry;
int ntheta = param.ntheta;
int nphi = param.nphi;
tstack_push("Generate fijk 3D-ball");
cout << "Generate Cl(k) spectra and fill a 3D-ball" << endl;
//Class to perform 2D Ylm transform with the convention in lagSHT
TestTransform blst(Lmax,-1,-1,-1);
int npix = blst.Npix();
int nalms = blst.Nalm();
vector<r_8> fijk(npix*Nmax,0); // the 3D pixelization
Vector2d<r_8> Clkin(Lmax,Nmax,0); // Generate Cl(k) with k index of subshell
for(int k=0;k<Nmax;k++){ //shell index
// cout << "Shell[" << k << "]" << endl;
// Generate true Cl Profile as in 2D case but evolving in k-direction
vector<r_8> Clin(Lmax); // Cl on the subshell
for(int l=0;l<Lmax;l++){
r_8 mean = Lmax/3.*(1. + k/(Nmax-1.)); // a k dependance Lmax/3 -> 2Lmax/3
r_8 sigma = Lmax/20.; //fix
r_8 arg = (l-mean)/sigma;
arg *= arg;
Clkin(l,k) = exp(-arg/2.);
Clkin(l,k) /= (r_8)(2*l+1); //nomalisation
Clin[l] = Clkin(l,k);
}//l-loop
//From the Clin compute the Alm and generate a 2D map
vector< complex<r_8> > alm(nalms);
vector<r_8> shmap(npix,0);
blst.GenAlmFromCl(Clin,alm);
blst.Backward(alm,shmap);
//transfert the 2D Map into the 3D pixelization
// cout << "Transfert Start at fijk[" << k*npix << "]" << endl;
copy(shmap.begin(),shmap.end(),fijk.begin()+k*npix);
}//shell-loop
tstack_pop("Generate fijk 3D-ball");
tstack_push("Lag Spherical Transf.");
cout << "Perform the Laguerre Spherical Transf. to get both the flmn coef. and the alm on each shell" << endl;
//Perform 3D fijk -> flmn & flmk
LaguerreSphericalTransform
sphlagtrans(geom, //type of geometry
Lmax, //l=0,...,Lmaw-1
-1, //Mmax default = Lmax set if <0
Nmax, //Nmax : n:0...Nmax-1
Rmax, //Rmax : R in [0, Rmas]
ntheta, //number of collatitude rings (or nsides ) default computed if <0
nphi //number of phi per rings default computed if <0
);
int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
int Nshell = sphlagtrans.GetOrder();
int Ntot = Nshell * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
// int Npix = sphlagtrans.GetSphericalGeometry()->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> > flmn(Ntot);
vector< complex<r_8> > almk(Ntot);
sphlagtrans.Analysis(fijk,flmn, almk); // <---------- the job
tstack_pop("Lag Spherical Transf.");
tstack_push("Clk and Cln OUT");
cout << "Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn" << endl;
//Compute the Clkout with k the index of subshell and should be comapred to Clkin from the almk
//And the Cl(n) n the 3rd index of the flmn
//Compute back the Cl(l,n)
//and the Cl(l,k) from the alm(k)
Vector2d<r_8> Clnout(Lmax,Nmax,0);
Vector2d<r_8> Clkout(Lmax,Nmax,0);
for(int n=0;n<Nmax;n++){
for(int l=0;l<Lmax;l++){
for (int m=0;m<=l;m++) {
int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering
r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
Clnout(l,n) += (m==0) ? flmn2 : 2.0*flmn2;
Clkout(l,n) += (m==0) ? almk2 : 2.0*almk2;
}
Clnout(l,n) /= (r_8)(2*l+1); //nomalisation "a la Cl 2D"
Clkout(l,n) /= (r_8)(2*l+1);
}
}
tstack_pop("Clk and Cln OUT");
tstack_push("error analysis");
cout << "Error analysis on Cl(k) input and output of the Laguerre Spherical Transf." << endl;
r_8 err_abs(0.);
r_8 err_rel(0.);
pair<int,int> imax = make_pair(-1,-1);
for(int k=0;k<Nmax;k++){ //shell index
for(int l=0;l<Lmax;l++){
r_8 diff = Clkin(l,k)-Clkout(l,k); diff *=diff;
if (diff > err_abs) {
err_abs = diff;
imax = make_pair(k,l);
}
r_8 relatdiff = diff/fabs((r_8)Clkin(l,k));
if((relatdiff)>err_rel) err_rel = relatdiff;
}
}
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
cout << " Err. Max. " << err_abs << " [" << imax.first << ", "
<< imax.second << "], Err. Rel. " << err_rel << endl;
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
tstack_pop("error analysis");
// {//save
// #ifdef SOPHYALIB
// char buff[80];
// sprintf(buff,"test3-L%d-N%d.ppf",Lmax,Nmax);
// POutPersist po(buff);
// po << PPFNameTag("clkin") << Clkin;
// po << PPFNameTag("clkout") << Clkout;
// po << PPFNameTag("clnout") << Clnout;
// }//save
}//TEST 3
//-------------------------------------
// test5:
// o generated a 3D Cube with Nc^3 pixels filled with random gaussian number N(0,sigma) : fixed sigma
// o fill 3D Ball embedded into the Cube
// o perform a LagSHT Analysis to get the flmn and almk coefficients
// o compute the Cl(n) and the Cl(k) by reducing the "m" index. Notice that the Cl(k) are CMB like Cl on each
// subshells
void test5() {
using namespace LagSHT;
//input values
int Lmax = param.Lmax;
int Nmax = param.N;
r_8 Rmax = param.R; // <----------------- will be changed below
string geom = param.geometry;
int ntheta = param.ntheta;
int nphi = param.nphi;
//the 3D Cube of Nc^3 pixels
int state=1234567 + 8912;
int Nc = 256*4;
tstack_push("Generate a 3D-Cube");
cout << "Generate a 3D-Cube " << Nc << "^3" << endl;
Vector3d<r_8> cube(Nc,Nc,Nc);
r_8 rms = 1.0; // constant RMS for the moment
for(size_t ix=0; ix<(size_t)Nc; ix++) {
for(size_t iy=0; iy<(size_t)Nc; iy++) {
for(size_t iz=0; iz<(size_t)Nc; iz++) {
cube(ix,iy,iz) = (r_8)gaussrand(&state,0.,rms);
}
}
}
tstack_pop("Generate a 3D-Cube");
tstack_push("Fill the 3D-Ball from the Cube");
cout << "Fill the 3D-Ball from the Cube" << endl;
r_8 Lo = 256 * 8.0 ; //8Mpc per cell
r_8 Lx = Lo; //Total dimension of the Cube
r_8 Ly = Lo;
r_8 Lz = Lo;
//Gives physical dimensions of the Cube
r_8 Dx = Lx/Nc; //Mpc a cell X dimension
r_8 Dy = Ly/Nc; //Mpc a cell Y dimension
r_8 Dz = Lz/Nc; //Mpc a cell Z dimension
//define Rmax to include a 3D Ball in the 3D Cube
Rmax = min(min(Lx,Ly),Lz)/2.*0.8; //80% of the minimal half length of the cube
cout << " Fill a 3D-Ball (" << Dx << "," << Dy << "," << Dz << ")Mpc^3: "
<< "Rmax is set to : " << Rmax << endl;
LaguerreSphericalTransform
sphlagtrans(geom, //type of geometry
Lmax, //l=0,...,Lmaw-1
-1, //Mmax default = Lmax set if <0
Nmax, //Nmax : n:0...Nmax-1
Rmax, //Rmax : R in [0, Rmas]
ntheta, //number of collatitude rings (or nsides ) default computed if <0
nphi //number of phi per rings default computed if <0
);
BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
if(sphere == NULL) throw LagSHTError("test5 sphere == NULL");
int Nalm = sphere->Nalm();
int Nshell = sphlagtrans.GetOrder(); //number of shells
int Ntot = Nshell * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
int Npix = sphere->Npix(); //2D sphere #of pixels
int NpTot= Nshell * Npix; //totoal number of 3D-pixels
// cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
// << Npix << " "
// << NpTot << " "
// << Nalm << " "
// << Nshell << " "
// << Ntot << endl;
LaguerreTransform* lagtrans = sphlagtrans.GetLagTransform(); //the pure Laguerre Part
if(lagtrans == NULL) throw LagSHTError("test5 lagtrans == NULL");
vector<r_8> fijk(NpTot,0); // the 3D-ball
//scan the pixels of the different shells of the 3D Ball
for(int ks=0; ks<Nshell; ks++){ //shell k
r_8 rShell = lagtrans->Index2R(ks);
// cout << "Shell["<<ks<<"] radius = " << rShell << endl;
if(rShell > Rmax) throw LagSHTError("test5 a shell radius too large...");
for(int ipix=0; ipix<Npix; ipix++){ //pixels of the shell
r_8 theta; //here by convention theta = 0 means North Pole
r_8 phi;
sphere->PixThetaPhi(ipix,theta,phi);
//3D xyz position of the pixel
r_8 xpix = rShell * sin(theta) * cos(phi);
r_8 ypix = rShell * sin(theta) * sin(phi);
r_8 zpix = rShell * cos(theta);
//find the corresponding Cube cell ijk indices
int i = int((xpix + Lx/2.0 - Dx/2.0)/Dx);
if(i<0||i>Nc)throw LagSHTError("test5 i out of range");
int j = int((ypix + Ly/2.0 - Dy/2.0)/Dy);
if(j<0||j>Nc)throw LagSHTError("test5 j out of range");
int k = int((zpix + Lz/2.0 - Dz/2.0)/Dz);
if(k<0||k>Nc)throw LagSHTError("test5 k out of range");
//affect the 3D-cube ijk cell value to the shell ipix
int iter = ks*Npix+ipix; //location of the ipix of k-th shell
if(iter<0 || iter>NpTot)throw LagSHTError("test5 iter out of range");
fijk[iter] = cube(i,j,k);
}//end loop pixels of a shell
}