alm ok

python/spherelib/__init__.py

@@ -24,5 +24,8 @@
 
 import map
 from spherelib_map import regress, fastfill, cat2mask, apodize_mask
+from spherelib_alm import _alm2cov
 from fitsfunc2 import *
 from bin import *
+from alm import *
+from utils import *

python/spherelib/alm.py

+"""
+spherical harmonic coefficients 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 pylab import *
+import utils
+import spherelib_alm
+import bin
+
+def alm2cov (list_alm):
+    """ Compute second order statistics of spherical harmonic coefficients.
+    
+    Parameters
+    ----------
+    list_alm : list of alm array of same size. 
+    
+    Returns
+    -------
+    3D array representing symetric covariance matrices
+    vector of number of modes used to compute each matrix. 
+    """
+    
+    nin   = len(list_alm)
+    nalm  = len(list_alm[0])
+    almtab = zeros((nin, nalm))
+    
+    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))
+
+    spherelib_alm._alm2cov (almtab, covtab)
+
+    cov = zeros((nin, nin, lmax+1))
+    for ell in range(lmax+1):
+        cov[:,:,ell] = utils.vec2mat(covtab[:,ell])
+
+        
+    nmode = bin.get_nmode(bin.uniformbin (0, lmax, 1))
+    return cov, nmode
+

python/spherelib/cov.py

+""" second order statistics 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
+
+def unbeamcov (cov, beams, max_inv=10000):
+    """
+    """
+    pass

python/spherelib/fitsfunc2.py

@@ -37,21 +37,29 @@ def write_alm (filename, alm):
     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))
+    lmax   = int(round(((sqrt(1+8*size(alm))-3)/2)))
+    nalm = size(alm)
+    idx = zeros(size(alm), dtype="int32")
     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)
+        idx[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))))
+    l = floor(sqrt(idx-1)).astype(long)
+    m = idx - l**2 - l - 1
+    lmax = l.max()
+    mmax = m.max()
+    i = hp.Alm.getidx(lmax,l,m)
+
+    almr = alm.real[i]
+    almi = alm.imag[i]
+    
+    cols   = [pyfits.Column(name='index', format='J', array=idx)]
+    cols.append(pyfits.Column(name='real', format='D', array=almr))
+    cols.append(pyfits.Column(name='imaginary', format='D', array=almi))
             
     cols  = pyfits.ColDefs(cols)
     tbhdu = pyfits.new_table(cols)

python/spherelib/map.py

+"""
+Healpix map 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
 
 def rien():
     return 0

python/spherelib/spherelib_alm.cpp

+// 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
+
+#include "xcomplex.h"
+#include "alm.h"
+#include "powspec.h"
+#include "alm_powspec_tools.h"
+#include "alm_healpix_tools.h"
+#include "healpix_base.h"
+#include "healpix_map.h"
+#include "MapExt.h"
+#include "AlmExt.h"
+#include "PowSpecExt.h"
+#include "CovMat.h"
+#include "Python.h"
+#include <numpy/arrayobject.h>
+
+#define dcomplex std::double<complex>
+
+void _alm2cov(npy_cdouble *tabalm, int nbin, int nalm, double* cov, int ncross, int nell)
+{
+
+  int ncross2 = nbin*(nbin+1)/2;
+  if (ncross != ncross2)
+    PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",ncross, ncross2);
+  int nalm2 = (nell)*(nell+1)/2;
+  if (nalm != nalm2)
+    PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",nalm, nalm2);
+
+  int lmax = nell-1;
+  int mmax = lmax;
+  int j=0;
+  AlmExt<double> Alm[nbin];
+  for (int m=0; m< nbin; m++){
+    arr< xcomplex<double> > arralm (nalm);
+    for( int i = 0; i < nalm; i++ ){
+      arralm[i]=xcomplex<double>(tabalm[j].real, tabalm[j].imag);
+      j++;}
+    
+    Alm[m].Set(arralm, lmax, mmax);
+  }
+
+  CovMat<double> CovMat;
+  CovMat.FromAlm (Alm, nbin);
+  
+  arr2<double> Mat = CovMat.getCovMat();
+  j = 0;
+  for (int m=0; m< ncross; m++){
+    for( int i = 0; i < nell; i++ )
+      {
+	cov[j] = Mat[m][i];
+	j++;}
+  }
+}

python/spherelib/spherelib_alm.h

+/*
+** spherelib_alm.h
+// 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
+
+** 
+** Made by Maude Le Jeune
+** Login   <lejeune@localhost.localdomain>
+** 
+** Started on  Tue Jun 15 12:38:54 2010 Maude Le Jeune
+** Last update Tue Jun 15 12:38:54 2010 Maude Le Jeune
+*/
+
+#include <numpy/arrayobject.h>
+
+#ifndef   	SPHERELIB_ALM_H_
+# define   	SPHERELIB_ALM_H_
+
+void _alm2cov (npy_cdouble *alm, int nin, int nalm, double *cov, int ncross, int nell);
+
+#endif 	    /* !SPHERELIB_ALM_H_ */

python/spherelib/spherelib_alm.i

+// 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
+
+
+%module spherelib_alm
+%{
+#define SWIG_FILE_WITH_INIT
+#include "spherelib_alm.h"
+#include <numpy/arrayobject.h>
+%}
+
+%include "numpy.i"
+
+%init 
+%{
+  import_array();
+%}
+
+// Needed (see the documentation of numpy.i)
+%numpy_typemaps(npy_cdouble , NPY_CDOUBLE, int)
+
+// alm2stat
+%apply (npy_cdouble* IN_ARRAY2, int DIM1, int DIM2) {(npy_cdouble *alm, int nin, int nalm)};
+%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {(double *cov, int ncross, int nell)};
+%feature("docstring", "Compute second order statistics of spherical harmonic coefficients.") _alm2cov ;
+%feature("autodoc",0) _alm2cov;
+
+%include "spherelib_alm.h"
+
+
+
+
+

python/spherelib/spherelib_map.cpp

+// 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
+
 #include "mask.h"
 #include "healpix_base.h"
 #include "fastfiller.h"

python/spherelib/spherelib_map.h

+// 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
+
 /*
-** _spherelib_map.h
+** spherelib_map.h
 ** 
 ** Made by Maude Le Jeune
 ** Login   <lejeune@localhost.localdomain>

python/spherelib/spherelib_map.i

+// 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
+
+
 %module spherelib_map
 %{
 #define SWIG_FILE_WITH_INIT

python/spherelib/test.py

@@ -3,7 +3,8 @@ from pylab import *
 import healpy
 import spherelib
 
-#map = healpy.read_map('/home/lejeune/Bureau/unecarte.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/unmask.fits')
 
 #healpy.mollview(squeeze(map))
@@ -12,14 +13,27 @@ import spherelib
 
 #healpy.mollview(map)
 
-map = ones((12*512*512, 1))
+#map = ones((12*512*512, 1))
 
 alm = healpy.map2alm(squeeze(map), 200)
+alm2 = healpy.map2alm(squeeze(map2), 200)
 
-spherelib.write_alm ('alm.fits', alm)
+spherelib.write_alm("a.fits", alm)
+b = healpy.read_alm ("a.fits")
+
+print alm
+print b
+
+cov,nm = spherelib.alm2cov([alm, alm2])
+
+figure(1)
+semilogy(squeeze(cov[1,1,:]))
+
+# 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,radius=50)
+mask = spherelib.cat2mask(mask,theta, phi,50)
+

python/spherelib/utils.py

+"""
+utils. 
+"""
+
+# 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 pylab import *
+
+def vec2mat (V):
+    """ Convert vector to symmetric matrix. 
+    """
+    
+    ncross = size(V)
+    dim    = int(floor(sqrt(ncross*2)))
+    M = zeros((dim, dim))
+    c = 0
+    for i in range(dim):
+        for j in range(i+1):
+            M[i,j] = V[c]
+            M[j,i] = M[i,j]
+            c = c+1
+    return M
+
+def mat2vec (M):
+    """ Convert symmetric matrix to vector. 
+    """
+    dim = shape(M)[0]
+    ncross = dim*(dim+1)/2
+    V = zeros((ncross, 1))
+    c = 0
+    for i in range(dim):
+        for j in range(i+1):
+            V[c] = M[i,j] 
+            c = c+1
+    return squeeze(V)

python/wscript

@@ -50,7 +50,7 @@ def build(bld):
 
 	build_healpix(bld)
 	
-	libs = ['cxxsupport', 'cfitsio', 'fftpack', 'gsl', 'gslcblas', 'healpix_cxx']
+	libs = ['gsl', 'gslcblas']
 	import Options
 	if bld.env["HEALPIX_TARGET"]=="gcc_omp":
 		libs.append('gomp')
@@ -64,24 +64,39 @@ def build(bld):
                 includes     = ['../lib/src', '../healpix/'+target+'/include'], 
 		defines      = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''], 
                 ccflags      = bld.env['CXXFLAGS'],
-                lib          = libs, 
-                libpath      = ['../healpix/'+target+'lib'],
+		lib          = libs,
+		staticlib    = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
+		staticlibpath = ['../healpix/'+target+'lib'],
 		name         = 'spherelib'
 		)
 	cxx_spherelib.install_path = '${PREFIX}/lib'
 
-	swig_spherelib = bld(
+	swig_spherelib_map = bld(
 		features = 'cxx cshlib pyext',
 		includes = ['spherelib', '../include', '../lib/src'],
 		source = ['spherelib/spherelib_map.i','spherelib/spherelib_map.cpp'],
 		target = '_spherelib_map',
 		swig_flags = '-c++ -python -Wall',
-		lib          = libs, 
-		libpath      = ['../healpix/'+target+'lib'],
+		lib          = libs,
+		staticlib    = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
+		staticlibpath = ['../healpix/'+target+'lib'],
 		uselib_local = 'spherelib')
 	
-	swig_spherelib.install_path = '${PYTHONDIR}/spherelib'
+	swig_spherelib_alm = bld(
+		features = 'cxx cshlib pyext',
+		includes = ['spherelib', '../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',
+		lib          = libs,
+		staticlib    = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
+		staticlibpath = ['../healpix/'+target+'lib'],
+		uselib_local = 'spherelib')
 	
+	swig_spherelib_map.install_path = '${PYTHONDIR}/spherelib'
+	swig_spherelib_alm.install_path = '${PYTHONDIR}/spherelib'
+
 	obj = bld(features = 'py')
 	obj.find_sources_in_dirs( ['./spherelib', out+'/default/spherelib'], exts=['.py'] )
 	obj.install_path = '${PYTHONDIR}/spherelib'