Commit b48a74d8 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

add lib

parent 7241c879
This diff is collapsed.
CC= gcc -fPIC
CCFLAGS_NO_C= -W -Wall -I$(INCDIR) -O3 -march=native -g0 -s -fno-strict-aliasing -fomit-frame-pointer -ffast-math -fopenmp
CCFLAGS= $(CCFLAGS_NO_C) -c
CXX= g++ -fPIC
CXXL= g++ -fPIC
CXXWARNFLAGS= -W -Wall -Wshadow -Wwrite-strings -Wredundant-decls -Woverloaded-virtual -Wcast-qual -Wcast-align -Wpointer-arith -Wconversion -Wold-style-cast -Wno-unknown-pragmas
CXXCFLAGS_NO_C= $(CXXWARNFLAGS) -ansi -I$(INCDIR) -O3 -march=native -fopenmp -g0 -s -ffast-math -fomit-frame-pointer
CXXCFLAGS= $(CXXCFLAGS_NO_C) -c
CXXLFLAGS= -L. -L$(LIBDIR) -s -ffast-math -fopenmp -O3 -march=native
ifeq ($(shell uname),SunOS)
CXX_EXTRALIBS= -lnsl -lsocket
endif
ARCREATE= ar crv
include $(LEVELS_SRC)/config/rules.common
This diff is collapsed.
This diff is collapsed.
/*
* Project: Spherelib
* Date: 08/10
* Authors: J.M. Colley - M. Le Jeune
*
*/
/* Copyright (C) 2008 APC CNRS Université Paris Diderot <lejeune@apc.univ-paris7.fr>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/gpl.html
*/
#define ALM_DBG
#include <stdlib.h>
#include <string.h>
#include "CD_CommonDefinition.h"
#include "xcomplex.h"
#include "alm.h"
#include "healpix_base.h"
#include "fitshandle.h"
#include "alm_map_tools.h"
#include "healpix_data_io.h"
#include "AlmExt.h"
using namespace std;
/******************************************************************************
* Print alm values
*
*/
template <typename T>
void AlmExt<T>::PrintAlm(long lmax,long mmax,const string &title)
{
cout << "\n"<<title; // print title of the alm
// check on Alm consistency
if ((this->Lmax() ==0) && (this->Mmax() ==0)) {
CD_WarningStatus(1);
return;
}
if ( (lmax < 0) || (mmax < 0))
{ CD_WarningStatus(2);
return;
}
if ( (lmax > this->Lmax() ) ||(mmax > this->Mmax()) )
{
char Mes[500];
sprintf(Mes,"\nLmax = %d - %ld , Mmax = %d - %ld \n", this->Lmax(),lmax, this->Mmax(),mmax);
CD_WarningStatusMes(3, Mes);
if (lmax > this->Lmax())
lmax = this->Lmax();
if (mmax > this->Mmax())
mmax = this->Mmax();
}
//print alm value for each l and m up to lmax and mmax
for(long Myl=0; Myl <= lmax; Myl++)
{
for(long Mym=0; Mym <= Myl; Mym++)
cout <<"["<<Myl<<","<<Mym<<"]="<<(*this)(Myl,Mym)<<" ";
cout<<"\n";
}
return;
}
template void AlmExt<double>::PrintAlm(long lmax,long mmax,const string &title) ;
template void AlmExt<float>::PrintAlm(long lmax,long mmax,const string &title) ;
/******************************************************************************
* Print attribut of class Amm
*
*/
template <typename T>
void AlmExt<T>::PrintAttribut(const string &Name)
{
cout << "\nAlmExt object Name :'" << Name << "'";
cout << "\n Lmax : " << this->Lmax();
cout << " Mmax : " << this->Mmax();
cout << "\n PathWeight: '" << m_PathWeight << "'";
cout << "\n Status : " << m_Status;
//cout << "\n & arr : " << this->mstart();
PrintAlm(3, 3, " First Alm values: \n");
}
template void AlmExt<double>::PrintAttribut(const string &Name) ;
template void AlmExt<float>::PrintAttribut(const string &Name) ;
/******************************************************************************
*
* ScaleL exist also in base class with arr<T>
* //CD_return( (*this).ScaleL( (double *)&(Beam.tt()[0]) ) );
* Sept 2007, Jmc
*/
/******************************************************************************
*
* Compute the linear combination AlmIntern <- alpha.AlmIntern + beta.AlmB
*
*/
template <typename T>
Status AlmExt<T>::LinComb(const AlmExt<T> &AlmB, T alpha, T beta )
{
long Lmax = this->Lmax();
if (Lmax > AlmB.Lmax())
{
char Mes[500];
sprintf(Mes,"\nLmax = %ld - %ld\n", Lmax,AlmB.Lmax());
CD_WarningStatusMes(3, Mes);
Lmax = AlmB.Lmax();
}
this->Scale(alpha);
for (long ell=0; ell<=Lmax; ell++)
for (long emm=0; emm <=ell; emm++)
(*this)(ell, emm) += beta*AlmB(ell,emm);
CD_returnOK;
}
template Status AlmExt<double>::LinComb(const AlmExt<double> &AlmB, double alpha, double beta );
template Status AlmExt<float>::LinComb(const AlmExt<float> &AlmB, float alpha, float beta );
/******************************************************************************
*
* Compute the linear combination AlmIntern <- AlmIntern + beta.AlmB
*
*/
template <typename T>
Status AlmExt<T>::LinComb(const AlmExt<T> &AlmB,T beta )
{
(*this).LinComb(AlmB, 1, beta);
CD_returnOK;
}
template Status AlmExt<double>::LinComb(const AlmExt<double> &AlmB,double beta );
template Status AlmExt<float>::LinComb(const AlmExt<float> &AlmB,float beta );
/******************************************************************************
*
*/
template <typename T>
AlmExt<T>& AlmExt<T>::operator=(const AlmExt<T>& OtherAlm)
{
if (this != &OtherAlm)
{
this->Set(OtherAlm.Lmax(), OtherAlm.Mmax());
for (long Lidx=0; Lidx <= this->Lmax() ; Lidx++)
for (long Midx=0 ; Midx <= Lidx ; Midx++)
(*this)(Lidx, Midx) = OtherAlm(Lidx, Midx);
this->m_PathWeight = OtherAlm.m_PathWeight;
}
return *this;
}
template AlmExt<double>& AlmExt<double>::operator=(const AlmExt<double>& OtherAlm);
template AlmExt<float>& AlmExt<float>::operator=(const AlmExt<float>& OtherAlm);
/*
* Project: Spherelib
* Date: 08/10
* Authors: J.M. Colley - M. Le Jeune
*
*/
/* Copyright (C) 2008 APC CNRS Université Paris Diderot <lejeune@apc.univ-paris7.fr>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/gpl.html
*/
#ifndef ALMEXT_H
#define ALMEXT_H
#include "healpix_data_io.h"
#include "alm.h"
#include "trafos.h"
#include "alm_powspec_tools.h"
using namespace std;
// circular dependency between AlmExt and HealpixMapExt
// anticipated declaration of AlmExt
template <typename T> class MapExt;
template <typename T> class PowSpecExt;
/////////////////////////////////////////////////
/// AlmExt is an extension of the Healpix class alm
/////////////////////////////////////////////////
template<typename T>
class AlmExt : public Alm< xcomplex<T> >
{
public:
/////////////////////////////////////////////////
/// Heavy constructor that load the alm in memory
/////////////////////////////////////////////////
AlmExt(int lmax_=0, int mmax_=0) : Alm<xcomplex<T> >(lmax_,mmax_)
{InitAttribut();}
/////////////////////////////////////////////////
/// Destructor
/////////////////////////////////////////////////
virtual ~AlmExt() {}
// ====
// FUNC
// ====
// Getter
/////////////////////////////////////////////////
/// Return status of last call method
Status GetStatus() {return(m_Status);}
// Setter
/////////////////////////////////////////////////
/// Define path to Healpix ring weight files
/// \param PathWeight : path to the Healpix ring weight files
/////////////////////////////////////////////////
void SetPathWeight(const string &PathWeight)
{m_PathWeight = PathWeight;}
/////////////////////////////////////////////////
/// Copy alm values
/// \param OtherAlm : AlmExt object
/////////////////////////////////////////////////
AlmExt& operator=(const AlmExt& OtherAlm);
// Pretty print
/////////////////////////////////////////////////
/// Print alm values
/// \param lmax : maximum order of l to print
/// \param mmax : maximum order of m to print
/// \param title : title of the alm
/////////////////////////////////////////////////
void PrintAlm(long lmax, long mmax, const string &title="");
/////////////////////////////////////////////////
/// Print alm values
/// \param title : title of the alm
/////////////////////////////////////////////////
void PrintAlm(const string &title="")
{PrintAlm(this->Lmax(), this->Mmax(), title);}
/////////////////////////////////////////////////
/// Print attribut of class
/// \param Name : title of the alm
/////////////////////////////////////////////////
void PrintAttribut(const string &Name="");
// Check
/////////////////////////////////////////////////
/// Check compatibility between two Alm object
/// return 0 if OK
/// \param ObjAlm : other Alm
Status Conformable(const AlmExt<T> &ObjAlm) {
if (!Alm<xcomplex<T> >::conformable(ObjAlm)){ CD_return(2); }
CD_returnOK;}
// Alm <- SHT -> Map : SHT : Spherical Harmonic Transformation
/////////////////////////////////////////////////
/// Create a Map object computed from alm using iSHT
/// \param Nside : resolution of the map
/// \param scheme : ordering of the map
/////////////////////////////////////////////////
MapExt<T> * ToMap(int Nside,
Healpix_Ordering_Scheme scheme=RING)
{
MapExt<T> *pRetMap = new MapExt<T>;
if (pRetMap == NULL) {
CD_WarningStatus(7);
return(NULL);
}
// Call Map method to compute iSHT of Alm
Status status;
status = pRetMap->FromAlm(*this,Nside,scheme);
if ( status != 0) {
delete pRetMap;
CD_WarningStatus(status);
return(NULL);
}
return(pRetMap);
}
/// Initialize poalrized alm values computing an SHT from file map with current Lmax
/// and Mmax of AlmExt object
/// \param mapT : temperature map
/// \param mapQ : Q polar map
/// \param mapU : U polar map
/// \param AlmE : E modes
/// \param AlmB : B modes
/// \param N_iter : number of iterations (default is 0)
/////////////////////////////////////////////////
Status FromPolarMap(MapExt<T> &mapT,MapExt<T> &mapQ,MapExt<T> &mapU, AlmExt<T> &AlmE, AlmExt<T> &AlmB, int N_iter = 0)
{
// check that path is set
if (m_PathWeight=="") {CD_return(6);}
arr<double> weight_T, weight_Q, weight_U;
read_weight_ring (m_PathWeight, mapT.Nside(), weight_T);
// Add 1 like anafast
for (int m=0; m<weight_T.size(); ++m) {
weight_T[m] +=1 ; }
// replace undef with 0
mapT.SetUndefTo0();
mapQ.SetUndefTo0();
mapU.SetUndefTo0();
double avg=mapT.average();
mapT.add(-avg);
if (mapT.Scheme()==NEST) mapT.swap_scheme();
if (mapQ.Scheme()==NEST) mapQ.swap_scheme();
if (mapU.Scheme()==NEST) mapU.swap_scheme();
map2alm_pol_iter (mapT,mapQ,mapU,(*this),AlmE,AlmB,N_iter,weight_T);
(*this)(0,0) += avg*sqrt(fourpi);
mapT.add(avg);
CD_returnOK;
}
Trafo maketrafo (int num)
{
switch (num)
{
case 1: return Trafo(2000,2000,Equatorial,Galactic);
case 2: return Trafo(2000,2000,Galactic,Equatorial);
case 3: return Trafo(2000,2000,Equatorial,Ecliptic);
case 4: return Trafo(2000,2000,Ecliptic,Equatorial);
case 5: return Trafo(2000,2000,Ecliptic,Galactic);
case 6: return Trafo(2000,2000,Galactic,Ecliptic);
case 7: return Trafo(1950,1950,Equatorial,Galactic);
case 8: return Trafo(1950,1950,Galactic,Equatorial);
case 9: return Trafo(1950,1950,Equatorial,Ecliptic);
case 10: return Trafo(1950,1950,Ecliptic,Equatorial);
case 11: return Trafo(1950,1950,Ecliptic,Galactic);
case 12: return Trafo(1950,1950,Galactic,Ecliptic);
default: throw Message_error ("Unsupported transformation");
}
}
/////////////////////////////////////////////////
/// Rotate alm
/// \param rot : rotation parameter
/////////////////////////////////////////////////
Status Rotate (int trafo)
{
Trafo tr(maketrafo(trafo));
rotate_alm((*this),tr.Matrix());
}
/////////////////////////////////////////////////
/// Initialize alm values computing an SHT from file map
/// \param mapT : temperature map
/// \param mapQ : Q polar map
/// \param mapU : U polar map
/// \param AlmE : E modes
/// \param AlmB : B modes
/// \param Lmax : maximum order of l
/// \param Mmax : maximum order of m
/// \param N_iter : number of iterations (default is 0)
/////////////////////////////////////////////////
Status FromPolarMap(MapExt<T> &mapT,MapExt<T> &mapQ,MapExt<T> &mapU , AlmExt<T> &AlmE, AlmExt<T> &AlmB,int Lmax, int Mmax, int N_iter = 0)
{this->Set(Lmax, Mmax); AlmE.Set(Lmax, Mmax);AlmB.Set(Lmax, Mmax);CD_return(FromPolarMap(mapT,mapQ,mapU, AlmE, AlmB, N_iter)); }
/////////////////////////////////////////////////
/// Initialize alm values computing an SHT from object map with current Lmax
/// and Mmax of AlmExt object
/// \param Map : map to analyze
/// \param N_iter : number of iterations (default is 0)
/////////////////////////////////////////////////
Status FromMap(MapExt<T> &Map, int N_iter = 0)
{
arr<double> weight_T; // where to load ring weight coefficients
if (this->Lmax() == 0) {CD_WarningStatus(1);}
// check that path is set
if (m_PathWeight=="") {CD_return(6);}
// read ring weight
read_weight_ring (m_PathWeight, Map.Nside(), weight_T);
// Add 1 like anafast
for (int m=0; m<weight_T.size(); ++m)
weight_T[m] +=1 ;
// replace undef with 0
Map.SetUndefTo0();
// map minus mean like anafast
double avg = Map.average();
Map.add(-avg);
Healpix_Ordering_Scheme s = Map.Scheme();
if (s==NEST)
Map.swap_scheme();
map2alm_iter(Map, *this, N_iter,weight_T);
// like anafast
(*this)(0, 0) += avg*sqrt(fourpi);
// Add avg don't change map
Map.add(avg);
CD_returnOK;
}
/////////////////////////////////////////////////
/// Initialize alm values computing an SHT from object map
/// \param Map : map to analyze
/// \param Lmax : maximum order of l
/// \param Mmax : maximum order of m
/// \param N_iter : number of iterations (default is 0)
/////////////////////////////////////////////////
Status FromMap(MapExt<T> &Map, int Lmax, int Mmax, int N_iter = 0)
{
this->Set(Lmax, Mmax);
// return(FromMap(Map, N_iter));
CD_return(FromMap(Map, N_iter));
}
// Algèbre
/////////////////////////////////////////////////
/// Product all Alm for l=i with Vect[i]
/// \param Vect : Vect with Lmax size
/////////////////////////////////////////////////
template<typename T2> Status ScaleL(const T2 Vect[])
{
for(long Ll=0; Ll <= this->Lmax(); Ll++) {
for(long Lm=0; Lm <= Ll; Lm++) { (*this)(Ll, Lm) *= Vect[Ll]; }
}
CD_returnOK;
}
/////////////////////////////////////////////////
/// Product all Alm for l=i with Vect[i]
/// \param Vect : arr with Lmax size
/////////////////////////////////////////////////
template<typename T2> Status ScaleL (const arr<T2> Vect)
{
long Lmax = Vect.size()-1;
char Mes[500];
if (Lmax < this->Lmax())
{
sprintf(Mes,"\nLmax = %ld - %d\n", Lmax,this->Lmax());
CD_WarningStatusMes(3, Mes);
}
else
Lmax = this->Lmax();
for(long Ll=0; Ll <= Lmax; Ll++) {
for(long Lm=0; Lm <= Ll; Lm++) { (*this)(Ll, Lm) *= Vect[Ll]; }
}
CD_returnOK;
}
/////////////////////////////////////////////////
/// Perform a reconvolution on alms using beam from Beam and RefBeam up to Thresold
/// \param Beam : PowSpecExt object which contains beam values
/// \param Thresold : maximum value of deconvolver (default is 1000)
/// \param RefBeam : PowSpecExt object which contains reference beam values
/////////////////////////////////////////////////
Status Deconvol(arr<T> Beam, arr<T> RefBeam, double Thresold = 1000)
{
T inv;
arr<T> filter((int)(this->Lmax())+1);
// Inverse Beam with thresold
for(long Ll=0; Ll <= this->Lmax(); Ll++)
{
if (Ll > Beam.size())
inv = 0;//Thresold;
else if (Ll > RefBeam.size())
inv = 0;
else {//if ((Beam.tt(Ll)) > 0 && (RefBeam.tt(Ll)>0)) {
inv = RefBeam[Ll]/Beam[Ll];
if (inv > Thresold) inv = Thresold;}
filter[Ll] = inv;
}
// Scale by inverse beam
CD_return( (*this).ScaleL( filter));
}
// Filtering
/////////////////////////////////////////////////
/// Initialize alm values computing an L filtering of alm
/// \param ListAlm : collection of Alm to filter, dim NbAlm
/// \param Filter : collection of Lmax coefficient to apply, dim (NbAlm x Lmax)
/// \param NbAlm : number of Alm in ListAlm
/////////////////////////////////////////////////
template<typename T2> Status multiScaleL (AlmExt<T> **ListAlm,
arr2<T2> &Filter,
long NbAlm)
{
char Mes[500];
long F_lmax = 0;
long A_lmax = 0;
long Lmax = this->Lmax();
if (NbAlm > 0) {
// Check dimension between this and ListAlm
A_lmax = ListAlm[0]->Lmax();
if ((A_lmax < Lmax) || (ListAlm[0]->Mmax() < Lmax))
{
sprintf(Mes,"\nLmax = %ld - %ld\n", Lmax,A_lmax);
CD_WarningStatusMes(3, Mes);
Lmax = A_lmax;
}
// Check dimension between Alm Lmax and Filter Lmax
F_lmax = Filter.size2()-1;
if (F_lmax < Lmax)
{
sprintf(Mes,"\nLmax = %ld - %ld\n", Lmax,F_lmax);
CD_WarningStatusMes(3, Mes);
Lmax = F_lmax;
}
if (Filter.size1() < NbAlm)
{
sprintf(Mes,"\nInsufficient number of columns in filter (%ld - %ld)\n", Filter.size2(),NbAlm);
CD_returnMes(10, Mes);
}
// Perform filtering
(*this).SetToZero(); // set alm value to zeros first
for (long Idx=0;Idx < NbAlm; Idx++){
for (long Lidx=0; Lidx <= Lmax ; Lidx++)
for (long Midx=0 ; Midx <= Lidx ; Midx++)
(*this)(Lidx, Midx) += (Filter[Idx])[Lidx]*(*ListAlm[Idx])(Lidx, Midx) ;
}
CD_returnOK;
}
else
CD_WarningStatus(1);
}
/////////////////////////////////////////////////
/// Initialize alm values computing an L filtering of alm
/// \param ListAlm : collection of Alm to filter, dim NbAlm
/// \param Filter : collection of Lmax coefficient to apply, dim (NbAlm x Lmax)
/// \param NbAlm : number of Alm in ListAlm
/////////////////////////////////////////////////
template<typename T2> Status multiScaleL (AlmExt<T> ListAlm[],arr2<T2> &Filter,long NbAlm)
{
AlmExt<T> *pAlm[NbAlm];
for(int Li=0; Li < NbAlm; Li++) {
//ListAlm[Li].PrintAttribut();
pAlm[Li] = &(ListAlm[Li]);
}
CD_return( multiScaleL(pAlm, Filter, NbAlm) );
}
// Linear combination
/////////////////////////////////////////////////
/// Compute the linear combination AlmIntern <- alpha.AlmIntern + beta.AlmB
/// \param AlmB : alm to combine with
/// \param alpha : coefficient applied to the intern alm
/// \param beta : coefficient applied to AlmB
/////////////////////////////////////////////////
Status LinComb (const AlmExt<T> &AlmB, T alpha , T beta);
/////////////////////////////////////////////////
/// Compute the linear combination AlmIntern <- AlmIntern + beta.AlmB
/// \param AlmB : alm to combine with
/// \param beta : coefficient applied to AlmB
/////////////////////////////////////////////////
Status LinComb (const AlmExt<T> &AlmB, T beta);