Commit 1e454347 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

+nilc

parent f8422b39
......@@ -22,7 +22,7 @@ from pylab import *
import utils
import spherelib_alm
import bin
import healpy as hp
## from http://www.phys.uu.nl/~haque/computing/WPark_recipes_in_python.html
## using numerical recipes
......@@ -51,32 +51,32 @@ def lmax2nside(lmax, level=2):
## nside >= (lmax+1)/sqrt(12)
## npix >= (lmax+1)^2
if level==1:
nside = 2 ** nextpow2( (lmax+1)/sqrt(12.0) )
nside = nextpow2( (lmax+1)/sqrt(12.0) )
## Healpix's advice
## nside >= lmax/2
## i.e. npix >= 3 * lmax^2
elif level==2:
nside = 2 ** nextpow2( lmax/2.0 )
nside = nextpow2( lmax/2.0 )
## More conservative
## nside >= lmax
elif level==3:
nside = 2 ** nextpow2( lmax )
nside = nextpow2( lmax )
else:
raise Exception ("level is 1, 2 or 3")
nside = max( nside, 4 ) ## because map2alm not happy with nside = 1 or 2
return nside
def alm2cov (list_alm):
def alm2cov (almtab):
""" Compute second order statistics of spherical harmonic coefficients.
Parameters
----------
list_alm : list of alm array of same size.
almtab : array of alm (nalm, nin)
Returns
-------
......@@ -84,13 +84,15 @@ def alm2cov (list_alm):
vector of number of modes used to compute each matrix.
"""
nin = len(list_alm)
nalm = len(list_alm[0])
almtab = zeros((nin, nalm), dtype="complex128")
nin = shape(almtab)[1]
nalm = shape(almtab)[0]
#nin = len(list_alm)
#nalm = len(list_alm[0])
#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]
#for i in range(nin):
# almtab[i,:] = list_alm[i]
ncross = nin*(nin+1)/2
......@@ -106,3 +108,5 @@ def alm2cov (list_alm):
nmode = bin.get_nmode(bin.uniformbin (0, lmax, 1))
return cov, nmode
......@@ -220,8 +220,10 @@ def smooth_cl (cl, dloverl=0.1):
ell = arange(0,lmax+1)
clsmooth = zeros(shape(cl));
for l in range(lmax+1):
llo = max(int(round(ell[l]*(1-dloverl/2))), 0)
lhi = min(int(round(ell[l]*(1+dloverl/2))), lmax)
llo = max(int(round(ell[l]*(1-dloverl/2.0))), 0)
lhi = min(int(round(ell[l]*(1+dloverl/2.0))), lmax)
if llo==lhi:
lhi=llo+1
clsmooth[l] = mean(cl[llo:lhi])
return clsmooth
......
......@@ -20,6 +20,7 @@
from pylab import *
import bin
import map
import healpy
def beam_correction (cov, beams, max_inv=10000):
""" Correct covariance matrices from beam effect
......@@ -129,3 +130,27 @@ def wiener (RModel):
for q in range(nbin):
W[comp,:,:,q] = dot(RModel[:,:,q,comp],inv(R[:,:,q]))
return W
def plotlocalizedfield (F, zone, figfile=None):
""" Plot an Healpix map of a localized field.
Parameters
----------
F : array-like, shape (1, N), with N number of zones.
zone : array-like, shape (npix, 1). An HEALPix map with pixel equal to zone number.
figfile : string. File name where to save the plot.
"""
nside = healpy.get_nside(zone)
npix = healpy.nside2npix(nside)
map = zeros((npix, 1))
if shape(F):
for p in range(npix):
map[p] = F[zone[p]]
else:
map = ones((npix, 1))*F
healpy.mollview(squeeze(map))
if figfile is None:
show()
else:
savefig(figfile)
......@@ -20,6 +20,7 @@ Healpix map computations.
import spherelib_map
from numpy import *
import utils
def mask2mll(mask, lmax):
""" Compute the mode coupling matrix.
......@@ -38,3 +39,81 @@ def mask2mll(mask, lmax):
return mll
def make_kp_zone (covers, refmap, nside=128):
""" Build kp-like zone.
Parameters
----------
covers: array-like, sky coverage of each zone.
refmap: reference map (like WMAP K band)
nside: integer, map resolution
"""
npix = 12*nside*nside
map = zeros((npix))
spherelib_map._make_kp_zone(map, covers, refmap)
return squeeze(map)
def make_grid_zone (nside_grid, nside=128):
""" Build healpix grid zone.
Parameters
----------
nside_grid: integer, grid resolution
nside: integer, map resolution
"""
npix = 12*nside*nside
map = zeros((npix))
spherelib_map._make_grid_zone(map, nside_grid)
return squeeze(map)
def map2cov (maptab, zone):
""" Compute second order statistics of healpix maps.
Parameters
----------
maptab : array of map (nin, npix)
zone: map which indices pixel sets from 0 to nz
Returns
-------
3D array representing symetric covariance matrices
vector of number of modes used to compute each matrix. (nin, nin, nz)
"""
nin = shape(maptab)[0]
npix = max(zone)+1
ncross = nin*(nin+1)/2
covtab = zeros((ncross, npix))
nmode = zeros((npix))
spherelib_map._map2cov (maptab, zone, covtab, nmode)
cov = zeros((nin, nin, npix))
for p in range(npix):
cov[:,:,p] = utils.vec2mat(covtab[:,p])
return cov, nmode
def combine_zone (maptab, zone, w):
""" Combine healpy maps by pixel sets.
Parameters
----------
maptab : array of map (nin, npix)
zone: map which indices pixel sets from 0 to nz
w : weight for combination (nin, nz)
Returns
-------
an healpix map (npix)
"""
nin = shape(maptab)[0]
npix = len(zone)
cmap = zeros((npix))
spherelib_map._combine_zone(maptab, zone, cmap, w)
return cmap
......@@ -183,7 +183,7 @@ def read_alm (file):
alm[j] = complex(vec[i], vec[i+1])
j = j+1
else:
alm = sp.read_alm(file)
alm = hp.read_alm(file)
return np.squeeze(alm)
......
......@@ -22,6 +22,7 @@
#include "MapExt.h"
#include "AlmExt.h"
#include "PowSpecExt.h"
#include "CovMat.h"
void regress(double *tabmap, int npix, double *tabmask, int npix2, int order=1)
{
......@@ -166,3 +167,115 @@ void _mask2mll (double* tabmap, int npix, double* mll, int nell1, int nell2)
}
void _make_kp_zone(double *tabmap, int npix,double* tabcovers, int nz,double* tabmapref, int npixref ){
arr<double> arrmapref(tabmapref, npixref);
MapExt<double> KMap(arrmapref, RING);
arr<double> arrmap(tabmap, npix);
MapExt<double> fMask(arrmap, RING);
arr<double> covers(tabcovers, nz);
KPMask (covers, fMask, KMap);
arrmap = fMask.Map();
}
void _make_grid_zone(double *tabmap, int npix , int nside_in){
arr<double> arrmap(tabmap, npix);
MapExt<double> Out(arrmap, RING);
if (nside_in==0)
Out.fill(0);
else
{
Healpix_Map<double> Mask ;
Mask.SetNside(nside_in, RING);
for (int pix=0; pix<Mask.Npix(); pix++)
Mask[pix] = pix;
Out.Import_upgrade(Mask);
}
arrmap = Out.Map();
}
void _map2cov(double* tabmap, int nbin, int npix, double* tabmask, int npix2, double* cov, int ncross, int nz, double* nmode, int nz2){
int ncross2 = nbin*(nbin+1)/2;
if (ncross != ncross2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",ncross, ncross2);
if (npix != npix2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix, npix2);
MapExt<double> Maps[nbin];
int j = 0;
int nside =(int)sqrt(npix/12);
for (int m=0; m< nbin; m++){
Maps[m].SetNside(nside, RING);
for (int p=0; p<npix; p++){
Maps[m][p] = tabmap[j];j++;}
}
arr<double> arrmask(tabmask, npix2);
MapExt<double> Mask(arrmask, RING);
CovMat<double> CovMat;
CovMat.FromMap (Maps, nbin, Mask);
arr2<double> Mat = CovMat.getCovMat();
arr<double> nbmode = CovMat.getNbMode();
int nz3 = Mat.size2();
if (nz2 != nz)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",nz, nz2);
if (nz2 != nz3)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",nz3, nz2);
j = 0;
for (int m=0; m< ncross; m++){
for( int i = 0; i < nz; i++ )
{
cov[j] = Mat[m][i];
j++;}
}
for( int i = 0; i < nz; i++ )
nmode[i] = nbmode[i];
}
void _combine_zone(double* tabmap, int nbin, int npix, double* tabmask, int npix2, double* cmap, int npix3, double* W, int nin2, int nz){
if (npix != npix2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix, npix2);
if (npix3 != npix2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix3, npix2);
if (nin2 != nbin)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",nbin, nin2);
int nside =(int)sqrt(npix/12);
arr<double> arrmask(tabmask, npix2);
MapExt<double> Mask(arrmask, RING);
MapExt<double> mapS ;
mapS.SetNside(Mask.Nside(), Mask.Scheme());
mapS.fill(0);
int j = 0;
if (nbin>0) {
for (int iin=0; iin < nbin; iin++) {
MapExt<double> map;
map.SetNside(nside, RING);
for (int p=0; p<npix; p++){
map[p] = tabmap[j];j++;}
mapS.LinCombZone(map, &W[iin*nz], Mask);}
}
arr<double> arrmap (cmap, npix);
arrmap = mapS.Map();
}
......@@ -32,5 +32,10 @@ void fastfill(double *map, int npix, double eps=1e-6, int niter=-1, int startsca
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);
void _make_kp_zone(double *map, int npix,double* covers, int nz, double* mapref, int npixref );
void _make_grid_zone(double *map, int npix , int nside_in);
void _map2cov(double* map, int nin, int npix, double* mask, int npix2, double* cov, int ncross, int nz, double* nmode, int nz2);
void _combine_zone(double* map, int nin, int npix, double* mask, int npix2, double* cmap, int npix3, double* W, int nin2, int nz);
#endif /* !_SPHERELIB_MAP_H_ */
......@@ -53,6 +53,31 @@
%feature("docstring", " Compute the mode coupling matrix.") _mask2mll ;
%feature("autodoc",0) _mask2mll;
// make_kp and grid _zone
%apply (double* INPLACE_ARRAY1, int DIM1) {(double* map, int npix)};
%apply (double* IN_ARRAY1, int DIM1) {(double* mapref, int npixref)};
%apply (double* IN_ARRAY1, int DIM1) {(double* covers, int nz)};
%feature("docstring", "Build kp like zones") _make_kp_zone ;
%feature("autodoc",0) _make_kp_zone;
%feature("docstring", "Build grid like zones") _make_grid_zone ;
%feature("autodoc",0) _make_grid_zone;
// map2cov
%apply (double* IN_ARRAY2, int DIM1, int DIM2) {(double* map, int nin, int npix)};
%apply (double* IN_ARRAY1, int DIM1) {(double* mask, int npix2)};
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {(double *cov, int ncross, int nz)};
%apply (double* INPLACE_ARRAY1, int DIM1) {(double *nmode, int nz2)};
%feature("docstring", "Compute second order statistics of healpix map.") _map2cov ;
%feature("autodoc",0) _map2cov;
// combinezone
%apply (double* IN_ARRAY2, int DIM1, int DIM2) {(double* map, int nin, int npix)};
%apply (double* IN_ARRAY2, int DIM1, int DIM2) {(double* W, int nin2, int nz)};
%apply (double* IN_ARRAY1, int DIM1) {(double* mask, int mpix2)};
%apply (double* INPLACE_ARRAY1, int DIM1) {(double * cmap, int npix3)};
%include "spherelib_map.h"
......
......@@ -3,12 +3,20 @@ 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')
map = healpy.read_map('/home/lejeune/Bureau/map1.fits')
map2 = healpy.read_map('/home/lejeune/Bureau/map2.fits')
mask = healpy.read_map('/home/lejeune/Bureau/mask.fits')
mll = spherelib.mask2mll(mask, 500)
mll[0:3, 0:3]
npix = len(map)
maps = zeros((npix, 2))
maps[:,0] = map
maps[:,1] = map2
(cov, nmode) = spherelib.map2cov (maps, mask)
#mll = spherelib.mask2mll(mask, 500)
#mll[0:3, 0:3]
#healpy.mollview(squeeze(map))
......
......@@ -16,12 +16,13 @@ def set_options(ctx):
def configure(ctx):
print('→ configuring spherelib')
import Options
ctx.env['CXXFLAGS'] = ['-O2']
print('→ omp is ' + str(Options.options.omp))
if Options.options.omp:
Options.options.healpix_target='gcc_omp'
print('→ healpix_target is ' + str(Options.options.healpix_target))
ctx.env['HEALPIX_TARGET'] = str(Options.options.healpix_target)
ctx.env['CXXFLAGS'] = ['-O2']
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_gsl)+"/include")
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_numpy_core)+"/include")
if Options.options.omp:
......@@ -58,7 +59,7 @@ def build(bld):
print('→ building spherelib')
build_healpix(bld)
hp_data = 'HEALPIXDATA="'+bld.path.abspath()+'/../healpix/data"'
libs = ['gsl', 'gslcblas']
import Options
import os
......@@ -72,7 +73,7 @@ def build(bld):
target = 'spherelib',
vnum = get_version(),
includes = ['../lib/src', '../healpix/'+target+'/include'],
defines = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''],
defines = [hp_data],
ccflags = bld.env['CXXFLAGS'],
lib = libs,
uselib = libs,
......@@ -87,9 +88,10 @@ def build(bld):
includes = ['../lib/src', '../healpix/'+target+'/include','spherelib',
'../include', '../lib/src'],
source = ['spherelib/spherelib_map.i','spherelib/spherelib_map.cpp'],
defines = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''],
target = '_spherelib_map',
swig_flags = '-c++ -python -Wall',
ccflags = bld.env['CXXFLAGS'],
defines = [hp_data],
lib = libs,
libpath = [os.path.abspath('../healpix/'+target+'/lib')],
staticlib = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
......@@ -101,8 +103,8 @@ def build(bld):
'../include', '../lib/src'],
source = ['spherelib/spherelib_alm.i','spherelib/spherelib_alm.cpp'],
target = '_spherelib_alm',
defines = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''],
swig_flags = '-c++ -python -Wall',
defines = [hp_data],
lib = libs,
libpath = [os.path.abspath('../healpix/'+target+'/lib')],
staticlib = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
......@@ -114,8 +116,8 @@ def build(bld):
'../include', '../lib/src'],
source = ['spherelib/spherelib_bin.i','spherelib/spherelib_bin.cpp'],
target = '_spherelib_bin',
defines = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''],
swig_flags = '-c++ -python -Wall',
defines = [hp_data],
lib = libs,
libpath = [os.path.abspath('../healpix/'+target+'/lib')],
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