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

minor bugs

parent 1c2cf39d
......@@ -52,12 +52,12 @@ template <typename T> Status KPMask (arr<T> covers, MapExt<T> &Map, const MapExt
int lmax = 2*KMap.Nside();
AlmExt<T> uKAlm (lmax,lmax);
arr<T> Bl (lmax+1);
GaussBeam (Bl, 120);
GaussBeam (Bl, 480);
PowSpecExt<T> Beam (1, lmax);
Beam.Set (Bl);
MapExt<T> Mask ; Mask.SetNside(KMap.Nside(), RING); Mask.fill(-1);
MapExt<T> dKMap ; dKMap.SetNside(64, RING);
MapExt<T> dKMap ; dKMap.SetNside(16, RING);
MapExt<T> uKMap ; uKMap.SetNside(KMap.Nside(), RING);
for (int c=0; c<ncovers; c++)
......
......@@ -73,7 +73,7 @@ void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tab
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
}
// fit is ok, set to new value
if (chisqrdof < maxxi2 && guess[0]>0){
if (chisqrdof < maxxi2 && guess[0]>0 && err[0]!=INFINITY){
addGaussianPatternOnSky(hmap, pointing(guess[1],guess[2]), fwhm, guess[0], err[0]/10, true);
tabcat[icol+0] = guess[1];
tabcat[icol+1] = guess[2];
......
......@@ -316,6 +316,7 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
noisespec.resize(nbins,((double)sigmanoise*xsize*xsize));
}
for(int i = start; i < stop; i++){
printf("patch number %d\n", i);
sky->getMap(i, sqmap, true);
//if(show){
// ostringstream outfilename;
......
......@@ -58,7 +58,7 @@ def bin_stat (covmat, bin, nmode=None):
tmpcov = sum(tmpcov, axis=1).reshape(ndet, ndet)
ll = sum(nmodes,axis=0)
bincov[:,:,q] = tmpcov * (1.0/ll)
bincov[:,:,q] = tmpcov * (1.0/(1.0*ll))
lsize[q] = ll
return [bincov,lsize]
......
......@@ -70,9 +70,8 @@ def read_covmat (filename):
""" Read covariance matrices from FITS file.
FITS format is :
- 1st column : number of modes
- 2nd column : 1x1 auto spectra of 1st channel
- 3rd column : 2x1 cross spectra
- 1st column : 1x1 auto spectra of 1st channel
- 2nd column : 2x1 cross spectra
- ...
Parameters
......@@ -82,7 +81,6 @@ def read_covmat (filename):
Returns
-------
covmat : array-like, shape (m, m, N), with N number of covariance matrices.
nmodes : arraly-like, shape(1, N)
"""
hdulist = pyfits.open(filename)
......@@ -91,45 +89,38 @@ def read_covmat (filename):
ndet = floor(sqrt(ncross*2))
ndet = ndet.astype(int32)
nmodes = hdulist[1].data.field(0)
#nmodes = hdulist[1].data.field(0)
covmat = zeros((ndet, ndet, lmax+1))
c = 0
for i in range(ndet):
for j in range(i+1):
covmat[i,j,:] = reshape(hdulist[1].data.field(c+1),(1,1,lmax+1))
covmat[i,j,:] = reshape(hdulist[1].data.field(c),(1,1,lmax+1))
covmat[j,i,:] = covmat[i,j,:]
c = c+1
hdulist.close()
return [covmat, nmodes]
return covmat
def write_covmat (filename, covmat, nbmodes=None):
def write_covmat (filename, covmat):
""" Write covariance matrices to FITS file.
FITS format is :
- 1st column : number of modes
- 2nd column : 1x1 auto spectra of 1st channel
- 3rd column : 2x1 cross spectra
- 1st column : 1x1 auto spectra of 1st channel
- 2nd column : 2x1 cross spectra
- ...
Parameters
----------
filename : string. FITS file name.
covmat : array-like, shape (m, m, N), with N number of covariance matrices.
nbmodes : arraly-like, shape(1, N)
"""
ndet = covmat.shape[0]
lmax = covmat.shape[2]-1
if nbmodes is None:
cols = [pyfits.Column(name='multipole', format='D', array=array(range(lmax+1)))]
else:
cols = [pyfits.Column(name='multipole', format='D', array=nbmodes)]
ncross = ndet*(ndet+1)/2
cols = []
c = 0
for i in range(ndet):
for j in range(i+1):
......
......@@ -241,15 +241,15 @@ def read_table (file):
"""
if __piolib(file):
try:
tab = pio.ReadTAB2DObject(file, 'PIODOUBLE', '')
tab = pio.ReadVECTObject(file, 'PIODOUBLE', '')
grpname = pio.GetGrpName(file)
grp = pio.OpenTAB2DGrp(grpname, "r")
ncols = pio.GetTAB2DColumnGrp(grp)
nrows = len(tab)/ncols
grp = pio.OpenVECTGrp(grpname, "r")
(ncols,com)= pio.ReadKeywordObject("m", 'PIOINT', file, grp)
(nrows,com)= pio.ReadKeywordObject("Q", 'PIOINT', file, grp)
cov = np.zeros((nrows, ncols))
for c in range(ncols):
cov[:,c] = np.squeeze(tab[c::ncols])
val = pio.CloseTAB2DGrp(grp)
val = pio.CloseVECTGrp(grp)
except pio.pioError,val:
return str(val)
else:
......@@ -269,28 +269,25 @@ def read_covmat (file):
"""
if __piolib(file):
try:
tab = pio.ReadTAB2DObject(file, 'PIODOUBLE', '')
tab = pio.ReadVECTObject(file, 'PIODOUBLE', '')
grpname = pio.GetGrpName(file)
grp = pio.OpenTAB2DGrp(grpname, "r")
ncols = pio.GetTAB2DColumnGrp(grp)
nrows = len(tab)/ncols
nmodes = np.squeeze(tab[0::ncols])
cov = np.zeros((nrows, ncols-1))
for c in range(ncols-1):
cov[:,c] = np.squeeze(tab[c+1::ncols])
val = pio.CloseTAB2DGrp(grp)
ndet = np.int(np.sqrt(ncols*2))
stats = np.zeros ((ndet, ndet, nrows))
for q in range(nrows):
stats[:,:,q] = sp.vec2mat(cov[q,:])
cov = stats
grp = pio.OpenVECTGrp(grpname, "r")
(ncols,com)= pio.ReadKeywordObject("m", 'PIOINT', file, grp)
(nrows,com)= pio.ReadKeywordObject("Q", 'PIOINT', file, grp)
cov =((ncols, ncols, nrows))
c = 0
for i in range(ncols):
for j in range(ncols):
stats[i,j,q] = np.squeeze(tab[c::ncols])
c = c+1
val = pio.CloseVECTGrp(grp)
except pio.pioError,val:
return str(val)
else:
[cov,nmodes] = sp.read_covmat(file)
return [cov,nmodes]
cov = sp.read_covmat(file)
return cov
def write_covmat(file, cov, nbmodes=None):
def write_covmat(file, cov):
""" Write cov mat to file.
Parameters
......@@ -302,30 +299,24 @@ def write_covmat(file, cov, nbmodes=None):
if __piolib(file):
try:
nFreqs = np.shape(cov)[0]
ncross = nFreqs*(nFreqs+1)/2
Q = np.shape(cov)[2]
R = np.zeros((Q,ncross))
for q in range(Q):
R[q,:] = sp.mat2vec(cov[:,:,q])
grpname = pio.GetGrpName(file)
checkgroup (grpname, "TAB2D", str(ncross+1))
grp = pio.OpenTAB2DGrp(grpname, "w")
val = pio.CreateTAB2DObject(file,'PIODOUBLE')
nr = np.shape(R)[0]
nc = np.shape(R)[1]+1
if nbmodes is None:
w = np.zeros((Q))
else:
w = np.squeeze(nbmodes)
val = pio.WriteTAB2DObject(w,file, 'PIODOUBLE',"tab=0:0,0:%d"%(nr-1), grp)
val = pio.WriteTAB2DObject(R,file, 'PIODOUBLE',"tab=1:%d,0:%d"%((nc-1),(nr-1)) , grp)
val = pio.CloseTAB2DGrp(grp)
pio.CloseTAB2DGrp(grp)
checkgroup (grpname, "VECT", str(nFreqs*nFreqs*Q))
grp = pio.OpenVECTGrp(grpname, "w")
val = pio.CreateVECTObject(file,'PIODOUBLE')
val = pio.WriteKeywordObject(nFreqs, "number of channels", "m", 'PIOINT', file, grp)
val = pio.WriteKeywordObject(Q, "number of bins", "Q", 'PIOINT', file, grp)
tab = np.zeros((0))
for i in range(ncols):
for j in range(ncols):
tab = np.hstack ((tab, cov[i,j,:]))
pio.WriteVECTObject (tab, file, 'PIODOUBLE','begin=0 &\n end='+str(nFreqs*nFreqs*Q),grp)
val = pio.CloseVECTGrp(grp)
pio.CloseVECTGrp(grp)
except pio.pioError,val:
return str(val)
else:
sp.write_covmat(file, cov, nbmodes=nbmodes)
sp.write_covmat(file, cov)
......@@ -342,11 +333,13 @@ def write_table(file, tab):
nr = np.shape(tab)[0]
nc = np.shape(tab)[1]
grpname = pio.GetGrpName(file)
checkgroup (grpname, "TAB2D", str(nc))
grp = pio.OpenTAB2DGrp(grpname, "w")
val = pio.CreateTAB2DObject(file,'PIODOUBLE')
val = pio.WriteTAB2DObject(tab,file, 'PIODOUBLE',"tab=0:%d,0:%d"%((nc-1),(nr-1)) , grp)
pio.CloseTAB2DGrp(grp)
checkgroup (grpname, "VECT", str(nc*nr))
grp = pio.OpenVECTGrp(grpname, "w")
val = pio.CreateVECTObject(file,'PIODOUBLE')
val = pio.WriteKeywordObject(nc, "number of channels", "m", 'PIOINT', file, grp)
val = pio.WriteKeywordObject(nr, "number of bins", "Q", 'PIOINT', file, grp)
pio.WriteVECTObject (np.reshape(tab, nr*nc), file, 'PIODOUBLE','begin=0 &\n end='+str(nr*nc),grp)
pio.CloseVECTGrp(grp)
except pio.pioError,val:
return str(val)
else:
......
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