preproc ok

python/spherelib/.#wscript

-lejeune@localhost.localdomain.4210:1275548606
\ No newline at end of file

python/spherelib/__init__.py

+#
+# Project: Spherelib
+# Date: 08/10  
+# Authors: M. Le Jeune
+#
+#
+# Copyright (C) 2008 APC CNRS Universite Paris Diderot <lejeune@apc.univ-paris7.fr>  
 # 
-#  This file is part of Spherelib.
+# 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.
 # 
-#  Spehrelib 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 2 of the License, or
-#  (at your option) any later version.
-# 
-#  Healpy 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 Healpy; if not, write to the Free Software
-#  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-# 
-#  For more information about Healpy, see http://code.google.com/p/healpy
-# 
-"""spherelib blabla
-"""
-
-import warnings
+# 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
+#
 
-try:
-    ImportWarning
-except NameError:
-    class ImportWarning(Warning):
-        pass
+"""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
-import spherelib_map
-
+from spherelib_map import regress
+from fitsfunc2 import *
+from bin import *

python/spherelib/bin.py

+"""
+spectral windows computations. 
+"""
+
+# Copyright (C) 2009 APC CNRS Universite 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
+
+
+from numpy import *
+
+def bin_stat (covmat, bin, nmode=None):
+    """ Perform the averaging of spectral covariance matrices over flat bins. 
+    
+    Matrices are weighted by their number of modes which is supposed to equal to 2.ell+1.
+
+    Parameters
+    ----------
+    covmat : array-like, shape (m,m,lmax+1). Concatenation of lmax+1 symmetric matrix of size m.
+    bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+    nmode: array-like, shape (lmax+1, 1). Number of modes used to build covariance matrices. 
+    Returns
+    -------
+    array-like, shape (m,m,nbin). Binned matrices. 
+    
+    """
+    nbin   = bin.shape[0]
+    ndet   = covmat.shape[0]
+    bincov = zeros((ndet, ndet, nbin))
+    lsize  = zeros((nbin,1))
+
+    for q in range(nbin):
+        lmin   = bin[q,0]
+        lmax   = bin[q,1]
+	rg     = arange(lmin,lmax+1).astype(int32)
+        if nmode is None:
+            nmodes = 2*rg+1
+        else:
+            nmodes = nmode[rg]
+        
+
+	tmpcov = covmat[:,:,rg].reshape(ndet*ndet, len(rg))
+        tmpcov = tmpcov * (ones((ndet*ndet,1))*nmodes)
+        tmpcov = sum(tmpcov, axis=1).reshape(ndet, ndet)
+        
+        ll      = sum(nmodes,axis=0)
+	bincov[:,:,q] = tmpcov * (1.0/ll)
+        lsize[q] = ll
+
+    return [bincov,lsize]
+
+
+def uniformbin (lstart, lstop, deltal):
+    """ Return description of flat bins of constant length. 
+    
+    Parameters
+    ----------
+    lstart : scalar, first multipole value. 
+    lstop  : scalar, last multipole value. 
+    deltal : scalar, length of bins. 
+
+    Returns
+    -------
+    array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+    """
+    L    = lstop - lstart + 1
+    nbin = ceil(L / deltal).astype(int32) 
+    bin  = zeros((nbin, 2))
+    
+    bin[0,0] = lstart
+    for q in range(1,nbin):
+        bin[q  , 0] = bin[q-1, 0] + deltal
+        bin[q-1, 1] = bin[q  , 0] - 1
+        
+    bin[-1,1] = lstop
+    return bin
+
+def wg2bin (lstart, lstop, mind=1):
+    """ Return description of flat bins as defined by WG2.  
+    
+    Parameters
+    ----------
+    lstart : scalar, first multipole value. 
+    lstop  : scalar, last multipole value. 
+    mind   : scalar, minimum bin length. 
+
+    Returns
+    -------
+    array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+    """
+    wg2_lmin = 0
+    wg2_lmax = 4000
+    wg2_ltrans = array([wg2_lmin,10,30,150,420,1200,2500,4000])
+    wg2_deltal = [1,2,5,10,20,50,100]
+    wg2_deltal = array(max (wg2_deltal, mind))
+    
+    ## cut ltrans between lmin and lmax
+    l1     = where(wg2_ltrans>=lstart)[0]
+    l2     = where(wg2_ltrans<=lstop)[0]
+    idx    = intersect1d(l1, l2)
+    ltrans = unique1d(hstack((hstack(([max(lstart-1,0)], wg2_ltrans[idx])), [lstop])))
+
+    ## cut deltal between lmin and lmax
+    if (len(idx)!=0):
+        deltal = hstack((array([wg2_deltal[idx[0]-1]]), wg2_deltal[idx]))
+        deltal = deltal[0:len(ltrans)-1]
+    else:
+        deltal = wg2_deltal[where(wg2_ltrans<lstop)[0]][-1]
+    
+  ## build uniform bins 
+    nbin = size(deltal)
+    if nbin==1:
+        deltal = [deltal]
+    bins = zeros((0,2))
+    for t in range(nbin):
+        b = uniformbin(ltrans[t]+1,ltrans[t+1],deltal[t])
+        bins = vstack((bins,b))
+    return bins
+        
+def get_lmean (bin):
+    """ Return mean multipole values of flat bins. 
+    
+    Parameters
+    ----------
+    bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+
+    Returns
+    -------
+    array-like, shape (nbin, 1), mean multipole values of each bin.
+    """
+    nbin  = len(bin)
+    lmean = zeros((nbin, 1)) 
+    for q in range(nbin):
+        lmean[q] = (bin[q,0] + bin[q,1] ) / 2
+    return lmean
+
+
+def get_nmode (bin):
+    """ Return the number of modes of flat bins. 
+    
+    Number of modes is supposed to equal to \f$\2ell+1\f$ for each mulitpole. 
+
+    Parameters
+    ----------
+    bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+
+    Returns
+    -------
+    array-like, shape (nbin, 1), number of modes of each bin.
+    """ 
+    nbin  = len(bin)
+    nmode = zeros((nbin, 1)) 
+    for q in range(nbin):
+        rg       = arange(bin[q,0] , bin[q,1]+1)
+        nmode[q] = sum (2 * rg +1, axis=0 )
+    return nmode
+
+
+def gaussbeam(fwhm_arcmin, lmax):
+    """Return the transfer window function of a Gaussian beam up to lmax.
+
+    Parameters
+    ----------
+    fwhm_arcmin : scalar. Full width at half maximum in arcminute. 
+    lmax : maximum multipole value. 
+
+    Returns
+    -------
+    array-like, shape (lmax+1, 1), beam transfert function. 
+    """
+    fwhm_rad = (array(fwhm_arcmin)*2*pi/60.) /360.
+    fwhm_rad = reshape(fwhm_rad,(1,size(fwhm_rad)))
+    
+    l = arange(lmax+1)
+    l = reshape(l,(size(l),1))
+    beam = sqrt(exp(-l*(l+1)/8./log(2) * fwhm_rad**2))
+    return beam
+
+def binavg(x, binsize):
+    """Bin a vector into bins of size binsize.
+    """
+    bin  = arange(0,size(x), binsize)
+    nbin = size(bin)
+    res = zeros(nbin)
+    for b in arange(nbin -1):
+        res[b] = mean(x[bin[b]:bin[b+1]-1])
+    res[nbin-1] = mean(x[bin[b]:])
+  
+    return res
+
+def smooth_cl (cl, dloverl=0.1):
+    """ Return smoothed power spectrum.
+
+    Power spectrum is averaged over window of size dloverl (in %)
+    
+    Parameters
+    ----------
+    cl : array-like, shape(lmax+1, 1). Power spectrum. 
+    dloverl : averaging window size. 
+
+    Returns
+    -------
+    array-like, shape(lmax+1, 1). Smoothed power spectrum. 
+    """
+    lmax = size(cl) - 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)
+        clsmooth[l] = mean(cl[llo:lhi])
+    return clsmooth
+
+
+def bin_cl (cl, bin):
+    """ Return binned power spectrum.
+
+    Power spectrum is averaged over flat bins.
+    
+    Parameters
+    ----------
+    cl : array-like, shape(lmax+1, 1). Power spectrum. 
+    bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+
+    Returns
+    -------
+    array-like, shape(nbin, 1). Binned power spectrum. 
+    """
+    nbin = shape(bin)[0]
+    cq   = zeros((nbin, 1))
+    for q in range(nbin):
+        lmin   = bin[q,0]
+        lmax   = bin[q,1]
+	rg     = arange(lmin,lmax+1).astype(int32)
+        lsize = float(size(rg))
+	cq[q] = sum(cl[rg])/lsize
+
+    return squeeze(cq)
+    
+    
+def unbin_cl (cq, bin):
+    """ Return unbinned power spectrum.
+
+    Supposed that power spectrum has been averaged over flat bins.
+
+    Parameters
+    ----------
+    cq : array-like, shape(nbin, 1). Binned power spectrum. 
+    bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
+
+    Returns
+    -------
+    array-like, shape(lmax+1, 1). Unbinned power spectrum. 
+    """
+    nbin = shape(bin)[0]
+    lmax = bin[nbin-1,1]
+    cl   = zeros((lmax+1, 1))
+    for q in range(nbin):
+        lmin   = bin[q,0]
+        lmax   = bin[q,1]
+	rg     = arange(lmin,lmax+1).astype(int32)
+	cl[rg]  = cq[q]
+
+    return squeeze(cl)
+    
+    

python/spherelib/fitsfunc2.py

+""" second order statistics FITS I/O's
+"""
+
+# Copyright (C) 2009 APC CNRS Universite 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
+
+
+
+import healpy as hp
+from numpy import *
+import pyfits
+
+
+def write_alm (filename, alm):
+    """ Write alm to FITS file. 
+
+    FITS format is : 
+    - 1st column : index value 
+    - 2nd column : real part of alm
+    - 3rd column : imaginary part of alm 
+    
+    Parameters
+    ----------
+    filename : string. FITS file name. 
+    alm : complex array, shape(1, (lmax+1)(lmax+2)/2)
+    """
+    
+    lmax   = int(round(((sqrt(1+8*size(alm))-3)/2)))
+
+
+    ##Building index
+    Ind = zeros(size(alm))
+    start = 0 
+    for l in range(lmax+1):
+        stop = start + l 
+        Ind[start:stop+1] = arange(l*l+l+1,(l+1)*(l+1)+1)
+        start = stop +1
+
+    cols   = [pyfits.Column(name='index', format='J', array=array(Ind))]
+    cols.append(pyfits.Column(name='real', format='E', array=array(real(alm))))
+    cols.append(pyfits.Column(name='imaginary', format='E', array=array(imag(alm))))
+            
+    cols  = pyfits.ColDefs(cols)
+    tbhdu = pyfits.new_table(cols)
+    tbhdu.writeto(filename)
+
+
+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 
+    - ...
+
+    Parameters
+    ----------
+    filename : string. FITS file name. 
+    
+    Returns
+    -------
+    covmat : array-like, shape (m, m, N), with N number of covariance matrices. 
+    nmodes : arraly-like, shape(1, N)
+    """
+    
+    hdulist = pyfits.open(filename)
+    ncross  = len(hdulist[1].columns.names)-1
+    lmax    = shape(hdulist[1].data.field(0))[0]-1
+    ndet    = floor(sqrt(ncross*2)) 
+    ndet    = ndet.astype(int32)
+
+    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[j,i,:] = covmat[i,j,:]
+            c = c+1
+
+
+    hdulist.close()    
+    return [covmat, nmodes]
+
+def write_covmat (filename, covmat, nbmodes=None):
+    """ 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 
+    - ...
+    
+    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
+
+    c = 0
+    for i in range(ndet):
+        for j in range(i+1):
+            cols.append(pyfits.Column(name='cross %d' % c, format='E', array=array(covmat[i,j,:])))
+            c = c+1
+            
+    cols  = pyfits.ColDefs(cols)
+    tbhdu = pyfits.new_table(cols)
+    tbhdu.writeto(filename)
+    
+    
+def read_table (filename, startcol=1):
+    """ Read table from FITS file.
+
+    Parameters
+    ----------
+    filename : string. FITS file name. 
+    startcol : number of 1st column read (default 1).
+
+    Returns
+    -------
+    array-like.
+    """
+    
+    hdulist = pyfits.open(filename)
+    ncols   = len(hdulist[1].columns.names)-startcol+1
+    nrows   = shape(hdulist[1].data.field(0))[0]
+
+    mat  = zeros((nrows, ncols))
+    for c in range(ncols):
+        mat[:,c] = hdulist[1].data.field(c+startcol-1)
+
+    hdulist.close()    
+    return mat
+    
+
+def write_table (filename, mat, colname=None):
+    """ Write table to FITS file.
+
+    Parameters
+    ----------
+    filename : string. FITS file name. 
+    mat :  array-like.
+    colname : string. Common prefix name for all columns
+    """
+    
+    ncols   = shape(mat)[1]
+    nrows   = shape(mat)[0]
+
+    if colname is None:
+        colname = []
+        for c in range(ncols):
+            colname.append("column "+str(c))
+
+    cols = []
+    for c in range(ncols):
+        cols.append(pyfits.Column(name=colname[c], format='E', array=array(mat[:,c])))
+
+    cols  = pyfits.ColDefs(cols)
+    tbhdu = pyfits.new_table(cols)
+    tbhdu.writeto(filename)
+    
+
+def read_key (filename, key):
+    """ Read keyword from FITS file.
+
+    Parameters
+    ----------
+    filename : string. FITS file name. 
+    key : string. Keyword name.
+    
+    Returns
+    -------
+    string. Key value as a string. 
+    """
+    
+    hdulist = pyfits.open(filename)
+    val     = hdulist[1].header[key]
+    hdulist.close()    
+    return val
+
+def write_cl(file, cl,startcol=1):
+    """Write power-spectra in a fits file.
+
+    Depending on the shape of cl, write TT or TT EE BB or TT EE BB TE
+
+    Parameters
+    ----------
+    file : string. FITS file name. 
+    cl : array-like, shape (lmax, 1) or (lmax, 3) or (lmax, 4). Power spectrum values. 
+    """
+    col1 = pyfits.Column(name='l',format='1J', array=arange(shape(cl)[0]))
+    if startcol==1:
+        colist= [col1]
+    else:
+        colist = []
+    if size(shape(cl)) == 1:
+        col2 = pyfits.Column(name='Temperature C_l',format='1D', array=cl[:])
+        colist.append(col2)
+    else:
+        col2 = pyfits.Column(name='Temperature C_l',format='1D', array=cl[:,0])
+        colist.append(col2)
+        if shape(cl)[1] >= 3:
+            col3 = pyfits.Column(name='E-mode C_l',format='1D', array=cl[:,1])
+            col4 = pyfits.Column(name='B-mode C_l',format='1D', array=cl[:,2])
+            colist.append(col3)
+            colist.append(col4)
+        if shape(cl)[1] == 4:
+            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.writeto(file,clobber=True)
+
+def read_cl(file,field=1,hdu=1):
+    """Read power-spectra from a fits file
+
+    Parameters
+    ----------
+    file : string. FITS file name. 
+    field : integer, column number. 
+    hdu : integer, hdu number. 
+
+    Returns
+    -------
+    cl : array-like, shape (lmax, 1).
+    """
+    hdulist = pyfits.open(file)
+    cl = hdulist[hdu].data.field(field)
+    return cl
+

python/spherelib/spherelib_map.cpp

@@ -8,7 +8,10 @@ void regress(double *tabmap, double *tabmask, int npix, int order=1)
   arr<double> arrmask(tabmask, npix);
   MapExt<double> map(arrmap, RING);
   MapExt<double> mask(arrmask, RING);
-  //Regress(map, mask);
+  Regress(map, mask, order);
+  map *= mask;
+  arrmap = map.Map();
+  
 }
 
 void fastfill(double *tabmap, int npix,double eps=1e-6, int niter=-1, int startscale=0 )

python/spherelib/test.py

-from spherelib import *
+
 from pylab import *
+import healpy
+import spherelib
+
+#map = healpy.read_map('/home/lejeune/Bureau/unecarte.fits')
+#mask = healpy.read_map('/home/lejeune/Bureau/unmask.fits')