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

Merge branch 'master' of gitlab.in2p3.fr:spherelib/spherelib

parents 6bc812fa 1e0f89b5
......@@ -53,7 +53,7 @@ def get_binMatrix(bins, btype='flat', lmax=None):
return Mbins
def bin_covariance(covmat, bins, nmode=None):
def bin_covariance(covmat, bins):
""" Perform the averaging of spectral covariance matrices over flat bins.
Matrices are weighted by their number of modes ie 2.ell+1.
......@@ -69,27 +69,24 @@ def bin_covariance(covmat, bins, nmode=None):
nbin = bins.shape[0]
ndet = covmat.shape[0]
bincov = np.zeros((ndet, ndet, nbin))
lsize = np.zeros((nbin))
if bins.shape[1]==2:
for q in range(nbin):
lmin = bins[q,0]
lmax = bins[q,1]
rg = np.arange(lmin,lmax+1).astype(int32)
nmodes = 2*rg+1 if nmode is None else nmode[rg]
nmodes = 2*rg+1
tmpcov = covmat[:,:,rg].reshape(ndet*ndet, len(rg))
tmpcov = tmpcov*(np.ones((ndet*ndet,1))*nmodes)
tmpcov = np.sum(tmpcov, axis=1).reshape(ndet, ndet)
ll = np.sum(nmodes,axis=0)
bincov[:,:,q] = tmpcov*(1.0/(1.0*ll))
lsize[q] = ll
return [bincov,lsize]
else: # use binning matrix
for q in range(nbin):
tmpcov = covmat.reshape(ndet*ndet, -1)
tmpcov = tmpcov*(np.ones((ndet*ndet,1))*bins[q,:])
tmpcov = np.sum(tmpcov, axis=1).reshape(ndet, ndet)
bincov[:,:,q] = tmpcov
return bincov
return bincov
def uniformbin(lstart, lstop, deltal):
""" Return description of flat bins of constant length.
......
......@@ -22,6 +22,8 @@ import bin
import map
import healpy
def beam_correction (cov, beams, max_inv=10000):
""" Correct covariance matrices from beam effect
......@@ -49,7 +51,7 @@ def beam_correction (cov, beams, max_inv=10000):
return ib
def mask_correction (cov, mask, corr="fsky", nmode=None, mll=None):
def mask_correction (cov, mask, corr="fsky", nmode=None, mll=None,fast=False):
"""
Parameters
......@@ -62,7 +64,7 @@ def mask_correction (cov, mask, corr="fsky", nmode=None, mll=None):
-------
scalar, fsky value.
square matrix (lmax+1,lmax+1) mll'.
"""
"""
if corr=="fsky":
fsky = mean(double(mask)**2)
cov /= fsky
......@@ -75,10 +77,19 @@ def mask_correction (cov, mask, corr="fsky", nmode=None, mll=None):
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,:]))
cov[j,i,:] = cov[i,j,:]
if fast:
cls = np.array([squeeze(cov[i,j,:] ) for i in range(nmap) for j in range(i,nmap)])
bin.spherelib_bin._mllsolve(cls, mll)
k = 0
for i in range(nmap):
for j in range(i,nmap):
cov[i,j,:] = cls[k]; k+=1 #linalg.solve(mll , squeeze(cov[i,j,:]))
cov[j,i,:] = cov[i,j,:]
else:
for i in range(nmap):
for j in range(i,nmap):
cov[i,j,:] = linalg.solve(mll , squeeze(cov[i,j,:]))
cov[j,i,:] = cov[i,j,:]
if nmode is not None:
nmode *= fsky
return fsky
......
......@@ -251,7 +251,8 @@ def write_cl(file, cl,startcol=1):
col5 = pyfits.Column(name='T-E cross-corr.',format='1D', array=cl[:,3])
colist.append(col5)
cols = pyfits.ColDefs(colist)
tbhdu = pyfits.new_table(cols)
#tbhdu = pyfits.new_table(cols)
tbhdu = pyfits.TableHDU.from_columns(cols)
tbhdu.writeto(file,clobber=True)
def read_cl(file,field=1,hdu=1):
......
......@@ -16,6 +16,11 @@
#include "arr.h"
#include "window.h"
#include "Python.h"
#include <iostream>
#include <Eigen/Dense>
#include <omp.h>
using namespace std;
using namespace Eigen;
void _make_prolate(int *tabbands, int nband, double* arrbins, int nwin, int nell, float* tabthetas, int nwin2,int rule){
......@@ -56,6 +61,34 @@ void _make_spline (int* tabbands, int nband, double* arrbins, int nwin, int nell
for (int j=0; j<nell; j++){
arrbins[k] = (double)bins[i][j];
k++;}
}
void _mllsolve (double* cls, int ncl, int nell, double* mll, int nell1, int nell2){
if (nell != nell1)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)", nell, nell1);
if (nell != nell2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)", nell, nell2);
MatrixXf A(nell,nell);
int k = 0;
for (int j=0; j<nell; j++){
for (int i=0; i<nell; i++){
A(j,i) = mll[k];
k++; } }
#pragma omp parallel for
for (int n=0; n<ncl; n++){
VectorXf b(nell);
VectorXf x(nell);
for (int j=0; j<nell; j++){
b(j) = cls[n*nell+j];
}
//x = A.colPivHouseholderQr().solve(b);
x = A.fullPivLu().solve(b);
for (int j=0; j<nell; j++)
cls[n*nell+j] = x(j);
}
}
......@@ -29,6 +29,5 @@
void _make_prolate(int* bands, int nband, double* bins, int nwin, int nell, float* thetas, int nwin2,int rule);
void _make_spline (int* bands, int nband, double* bins, int nwin, int nell, int order);
void _mllsolve (double* cls, int ncl, int nell, double* mll, int nell1, int nell2);
#endif /* !SPHERELIB_BIN_H_ */
......@@ -41,5 +41,11 @@
%feature("docstring", "Build spline window functions") _make_spline ;
%feature("autodoc",0) _make_spline;
// mll
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {(double* cls, int ncl, int nell)};
%apply (double* IN_ARRAY2, int DIM1, int DIM2) {(double* mll, int nell1, int nell2 )};
%feature("docstring", "Solve for mll") _mllsolve ;
%feature("autodoc",0) _mllsolve;
%include "spherelib_bin.h"
......@@ -12,6 +12,7 @@ def options(ctx):
ctx.add_option('--with_gsl', action='store', default='/usr', help='location of gsl')
ctx.add_option('--with_fftw3', action='store', default='/usr', help='location of fftw3')
ctx.add_option('--with_numpy_core', action='store', default=None, help='location of numpy core headers')
ctx.add_option('--with_eigen', action='store', default="/usr/include/eigen3", help='location of eigen headers')
ctx.tool_options('compiler_cxx')
ctx.tool_options('python')
......@@ -35,13 +36,13 @@ def configure(ctx):
ctx.env["HEALPIX_EXT_LIB"] = ctx.env["HEALPIX_EXT_PREFIX"]+'/lib'
ctx.env["HEALPIX_EXT_INC"] = ctx.env["HEALPIX_EXT_PREFIX"]+'/include'
ctx.env['CXXFLAGS'].append("-I"+ctx.env["HEALPIX_EXT_INC"])
ctx.check(header_name=['alm.h'],compile_mode='cxx', mandatory=True)
ctx.check(header_name=['alm.h'],compile_mode='cxx', mandatory=True)
ctx.check_cxx(lib=['cxxsupport', 'fftpack', 'healpix_cxx','sharp','c_utils'],libpath=ctx.env["HEALPIX_EXT_LIB"], uselib_store='healpix',mandatory=True)
print('→ find cfitsio from cfitsio_prefix')
ctx.env["CFITSIO_PREFIX"] = str(Options.options.cfitsio_prefix)
ctx.env["CFITSIO_EXT_LIB"] = ctx.env["CFITSIO_PREFIX"]
ctx.env["CFITSIO_EXT_INC"] = ctx.env["CFITSIO_PREFIX"]
ctx.env["CFITSIO_EXT_LIB"] = ctx.env["CFITSIO_PREFIX"]+"/lib"
ctx.env["CFITSIO_EXT_INC"] = ctx.env["CFITSIO_PREFIX"]+"/include"
ctx.env['CXXFLAGS'].append("-I"+ctx.env["CFITSIO_EXT_INC"])
ctx.check(header_name=['fitsio.h'],compile_mode='cxx', mandatory=True)
ctx.check_cxx(lib=['cfitsio'], libpath=ctx.env["CFITSIO_EXT_LIB"], uselib_store='cfitsio', mandatory=True)
......@@ -57,9 +58,12 @@ def configure(ctx):
except:
ctx.fatal('→ can\'t find numpy core headers')
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_eigen))
ctx.check(header_name=['Eigen/Dense'],compile_mode='cxx', mandatory=True)
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_gsl)+"/include")
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_fftw3)+"/include")
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_numpy_core)+"/include")
ctx.env['CXXFLAGS'].append("-I"+str(Options.options.with_numpy_core)+"/include")
ctx.env['CXXFLAGS'].append('-fPIC')
ctx.env['CXXFLAGS'].append('-fopenmp')
ctx.env['CXXFLAGS'].append('-fpermissive')
......
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