Commit 44eeb1f6 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

bug alm2stat

parent f5adceda
......@@ -97,16 +97,16 @@ class PowSpecExt : public PowSpec
}
// Map Object
Status FromMap(const MapExt<T> &map)
Status FromMap(MapExt<T> &map)
{
AlmExt<double> Alm;
Alm.Set(this->Lmax(), this->Lmax());
CD_iCheck( Alm.FromMap(map) );
CD_iCheck( Alm.FromMap(map) ) ;
CD_iCheck( FromAlm(Alm) );
CD_return(0);
}
Status FromMap(const MapExt<T> &map1,const MapExt<T> &map2)
Status FromMap(MapExt<T> &map1,MapExt<T> &map2)
{
AlmExt<double> Alm1;
AlmExt<double> Alm2;
......
......@@ -22,7 +22,7 @@
"""spherelib is an extension of the Healpix c++ library which provides tools for spherical harmonic and wavelet analysis on the sphere and filtering.
"""
import map
from map import *
from spherelib_map import regress, fastfill, cat2mask, apodize_mask
from spherelib_alm import _alm2cov
from fitsfunc2 import *
......
......@@ -38,12 +38,13 @@ def alm2cov (list_alm):
nin = len(list_alm)
nalm = len(list_alm[0])
almtab = zeros((nin, nalm))
almtab = zeros((nin, nalm), dtype="complex128")
lmax = int(round((sqrt(1+8*nalm)-3) / 2) )
for i in range(nin):
almtab[i,:] = list_alm[i]
lmax = int(round((sqrt(1+8*nalm)-3) / 2) )
ncross = nin*(nin+1)/2
covtab = zeros((ncross, lmax+1))
......
......@@ -19,8 +19,9 @@
from pylab import *
import bin
import map
def unbeamcov (cov, beams, max_inv=10000):
def beam_correction (cov, beams, max_inv=10000):
""" Correct covariance matrices from beam effect
Parameters
......@@ -28,6 +29,10 @@ def unbeamcov (cov, beams, max_inv=10000):
cov: 3D array, symetric covariance matrices
beams: vector, beam fwhm in arcmin
max_inv: upper limit for inversion
Returns
-------
2D array (nin, lmax+1) deconvolver
"""
nin = shape(cov)[0]
......@@ -41,9 +46,51 @@ def unbeamcov (cov, beams, max_inv=10000):
for ell in range(lmax+1):
cov[:,:,ell] = dot(dot(diag(ib[:,ell]),cov[:,:,ell]),diag(ib[:,ell]))
return ib
def mask_correction (cov, mask, corr="fsky", nmode=None):
"""
Parameters
----------
cov: 3D array, symetric covariance matrices
mask: array-like, healpix map
corr: string, correction to apply (fsky, mll)
Returns
-------
scalar, fsky value.
square matrix (lmax+1,lmax+1) mll'.
"""
if corr=="fsky":
fsky = mean(double(mask)**2)
cov /= fsky
if nmode is not None:
nmode *= fsky
return fsky
elif corr=="mll":
lmax = shape(cov)[2]-1
mll = map.mask2mll (mask, lmax)
fsky = sum(mll[0,:])
nmap = shape(cov)[0]
for i in range(nmap):
for j in range(i,nmap):
cov[i,j,:] = linalg.solve(mll , squeeze(cov[i,j,:]))
if nmode is not None:
nmode *= fsky
return fsky, mll
else:
error("incorrect correction keyword (should be 'fsky' or 'mll')")
def plotautocov (cov, figfile=None, leg=None):
"""
def plotautocov (cov, figfile=None, leg=None, tit=None, loc="upper right"):
""" Plot auto spectra of covariance matrices.
Parameters
----------
cov: 3D array, symetric covariance matrices
figfile: string, filename
leg: list of key string.
"""
nmap = shape(cov)[0]
figure()
......@@ -51,8 +98,34 @@ def plotautocov (cov, figfile=None, leg=None):
semilogy(squeeze(cov[m,m,:]))
if leg is not None:
legend(leg)
legend(leg, loc=loc)
if tit is not None:
title(tit)
if figfile is None:
show()
else:
savefig(figfile)
def wiener (RModel):
""" Compute multi-components wiener filter.
Parameters
----------
RModel : array-like, shape (m,m,N,C). 4D
array which contains covariance matrices for each bin (from 0 to
N) and each component (from 0 to C).
Returns
-------
array-like, shape (C,m,m,N).
"""
ndet = shape(RModel)[0]
ncomp = shape(RModel)[3]
nbin = shape(RModel)[2]
R = sum(RModel, axis=3)
W = zeros((ncomp, ndet, ndet, nbin))
for comp in range(ncomp):
for q in range(nbin):
W[comp,:,:,q] = dot(RModel[:,:,q,comp],inv(R[:,:,q]))
return W
......@@ -18,5 +18,32 @@ Healpix map computations.
# 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
def rien():
return 0
import spherelib_map
from numpy import *
def mask2mll(mask, lmax):
""" Compute the mode coupling matrix.
Parameters
----------
mask: array-like, healpix map
lmax: integer, maximum order of analysis
Returns
-------
square matrix (lmax+1, lmax+1).
"""
mll = zeros((lmax+1, lmax+1))
spherelib_map._mask2mll (mask, mll)
# for i in range(lmax+1):
# for j in range(0,i):
# mll[i,j] = mll[j,i]
# for i in range(lmax+1):
# for j in range(lmax+1):
# mll[i,j] = mll[i,j] * (2.0 * j +1.0)
return mll
......@@ -51,8 +51,11 @@ void _alm2cov(npy_cdouble *tabalm, int nbin, int nalm, double* cov, int ncross,
j++;}
Alm[m].Set(arralm, lmax, mmax);
}
Alm[0].PrintAlm(10, 10, "alm");
CovMat<double> CovMat;
CovMat.FromAlm (Alm, nbin);
......
......@@ -30,7 +30,7 @@
%}
// Needed (see the documentation of numpy.i)
%numpy_typemaps(npy_cdouble , NPY_CDOUBLE, int)
%numpy_typemaps(npy_cdouble , NPY_CDOUBLE, int);
// alm2stat
%apply (npy_cdouble* IN_ARRAY2, int DIM1, int DIM2) {(npy_cdouble *alm, int nin, int nalm)};
......
......@@ -19,6 +19,9 @@
#include "fastfiller.h"
#include "Python.h"
#include "apodizemask.cpp"
#include "MapExt.h"
#include "AlmExt.h"
#include "PowSpecExt.h"
void regress(double *tabmap, int npix, double *tabmask, int npix2, int order=1)
{
......@@ -98,4 +101,67 @@ void apodize_mask(double *tabmap, int npix,double radius=20,bool inside=false ,b
}
void _mask2mll (double* tabmap, int npix, double* mll, int nell1, int nell2)
{
if (nell2 != nell1)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",nell2, nell1);
arr<double> arrmap(tabmap, npix);
MapExt<double> map(arrmap, RING);
int lmax = nell1-1;
PowSpecExt<double> cl (1, lmax);
cl.FromMap (map);
arr<double> wl (2*lmax+1);
for (int ell=0; ell<lmax+1; ell++)
wl[ell] = cl.tt(ell);
//arr2<double> M (lmax+1, lmax+1);
lmax = 2*lmax;
double lrf[3*lmax+2];
lrf[0] = 0;
for(int i = 0; i <= 3*lmax; i++){
lrf[i+1] = lrf[i] + log(2. * (2. * i + 1.) / (i+1.) ) ;
}
double tmp ;
double pi_4 = 1. / (4 * 3.1415926535897931159979 );
int j=0;
for(int l1 = 0 ; l1 <= lmax/2 ; l1++){
for (int tl = 0; tl < l1 ; tl++){
mll[j] = mll[tl*(lmax/2 + 1)+l1] ; // set symetric
j++;
}
for(int l2 = l1 ; l2 <= lmax/2 ; l2++){
int minl = abs(l1-l2);
int maxl = l1+l2;
tmp = 0.0 ;
for(int l3 = minl ; l3 <= maxl ; l3++){
int L = l1+l2+l3;
int L_2 = L / 2 ;
if(2*L_2 == L){
double c = wl[l3] * (2 * l3 + 1);
tmp += exp( lrf[L_2 - l1]
+ lrf[L_2 - l2]
+ lrf[L_2 - l3]
- lrf[L_2]
- log(L + 1)) * c;
}
}
tmp = tmp * pi_4 ;
mll[j] = tmp;
j++;
}
}
j=0;
for(int l1 = 0; l1 <= lmax/2; l1++){
for(int l2 = 0; l2 <= lmax/2; l2++){
mll[j] = mll[j] * (2 * l2 +1) ; j++;
}
}
}
......@@ -31,5 +31,6 @@ void regress (double *map1, int npix1, double *map2, int npix2, int order=1);
void fastfill(double *map, int npix, double eps=1e-6, int niter=-1, int startscale=0);
void cat2mask(double *map, int npix, double* theta, int ntheta, double* phi, int nphi, double radius=20 );
void apodize_mask(double *map, int npix, double radius=20,bool inside=false ,bool pixbool=0,bool spline=false );
void _mask2mll(double *map, int npix, double* mll, int nell1, int nell2);
#endif /* !_SPHERELIB_MAP_H_ */
......@@ -48,6 +48,10 @@
%feature("docstring", "Convert a binary mask (0 or 1) of arbitrary shape into smooth mask. ") apodize_mask ;
%feature("autodoc",0) apodize_mask;
// mask2mll
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {(double* mll, int nell1, int nell2)};
%feature("docstring", " Compute the mode coupling matrix.") _mask2mll ;
%feature("autodoc",0) _mask2mll;
%include "spherelib_map.h"
......
......@@ -3,9 +3,12 @@ from pylab import *
import healpy
import spherelib
map = healpy.read_map('/home/lejeune/Bureau/mapilc512.fits')
map2 = healpy.read_map('/home/lejeune/Bureau/map512.fits')
#mask = healpy.read_map('/home/lejeune/Bureau/unmask.fits')
#map = healpy.read_map('/home/lejeune/Bureau/mapilc512.fits')
#map2 = healpy.read_map('/home/lejeune/Bureau/map512.fits')
mask = healpy.read_map('/home/lejeune/Bureau/mask.fits')
mll = spherelib.mask2mll(mask, 500)
mll[0:3, 0:3]
#healpy.mollview(squeeze(map))
......@@ -15,25 +18,25 @@ map2 = healpy.read_map('/home/lejeune/Bureau/map512.fits')
#map = ones((12*512*512, 1))
alm = healpy.map2alm(squeeze(map), 200)
alm2 = healpy.map2alm(squeeze(map2), 200)
# alm = healpy.map2alm(squeeze(map), 200)
# alm2 = healpy.map2alm(squeeze(map2), 200)
spherelib.write_alm("a.fits", alm)
b = healpy.read_alm ("a.fits")
# spherelib.write_alm("a.fits", alm)
# b = healpy.read_alm ("a.fits")
print alm
print b
# print alm
# print b
cov,nm = spherelib.alm2cov([alm, alm2])
# cov,nm = spherelib.alm2cov([alm, alm2])
figure(1)
semilogy(squeeze(cov[1,1,:]))
# figure(1)
# semilogy(squeeze(cov[1,1,:]))
# spherelib.write_alm ('alm.fits', alm)
# # spherelib.write_alm ('alm.fits', alm)
theta = array([0,0.2])
phi = array([0.2,0])
nside = 256
mask = ones((nside*nside*12,))
mask = spherelib.cat2mask(mask,theta, phi,50)
# theta = array([0,0.2])
# phi = array([0.2,0])
# nside = 256
# mask = ones((nside*nside*12,))
# mask = spherelib.cat2mask(mask,theta, phi,50)
......@@ -56,6 +56,7 @@ def build(bld):
libs.append('gomp')
target=bld.env["HEALPIX_TARGET"]
abspath = bld.root.abspath()
cxx_spherelib = bld(
features = ['cxx', 'cshlib'],
source = ['../lib/src/MapExt.cpp', '../lib/src/AlmExt.cpp', '../lib/src/Bin.cpp', '../lib/src/PowSpecExt.cpp', '../lib/src/CovMat.cpp'],
......@@ -76,6 +77,8 @@ def build(bld):
includes = ['spherelib', '../include', '../lib/src'],
source = ['spherelib/spherelib_map.i','spherelib/spherelib_map.cpp'],
target = '_spherelib_map',
#defines = ['HEALPIXDATA=\"'+abspath+'../healpix/data/\"'],
defines = ['HEALPIXDATA=\"/home/lejeune/Soft/Healpix_2.14a/data/\"'],
swig_flags = '-c++ -python -Wall',
lib = libs,
staticlib = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
......
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