Clean codes

setup.py

 #!/usr/bin/env python
 import sys
-#import hooks
 from distutils.core import setup, Extension
 import numpy
 from Cython.Distutils import build_ext
@@ -33,7 +32,7 @@ if USE_ICC:
         libs += ['gomp', 'iomp5']
         extra += ['-openmp']
 else:
-    extra += ['-O4']
+    extra += ['-O2']
     if USE_OPENMP:
         libs += ['gomp']
         extra += ['-fopenmp']
@@ -47,13 +46,12 @@ libcov = Extension( name="_libcov",
                     extra_compile_args=extra
                     )
 
-
-
 setup(name="xqml",
       description="Code for xQML",
       author="S. Vanneste & M. Tristram",
       version="0.1",
       packages=['xqml'],
       ext_modules=[libcov],
+      install_requires=["numpy","scipy>=1.1.0","healpy>=0.6.1"],
       cmdclass={'build_ext': build_ext}
       )

src/_libcov.pyx

@@ -11,6 +11,7 @@ np.import_array()
 
 
 cdef extern from "libcov.h":
+     void build_Wll( int nl, int npix, double* El, double* Pl, double* Wll)
      void build_dSdC( int nside, int nstokes, int npix, int inl, long *ispec, long *ellbins, long *ipix, double *bl, double* dSdC)
      void dlss( double X, int s1, int s2, int lmax, double *d)
      int polrotangle( double *ri, double *rj, double *cos2a, double *sin2a)
@@ -29,6 +30,13 @@ def dSdC( nside, nstokes,
 		 <double*> np.PyArray_DATA(bl),
 		 <double*> np.PyArray_DATA(dSdC))
 
+def CrossWindow( np.ndarray[double, ndim=3, mode="c"] El not None,
+		 np.ndarray[double, ndim=3, mode="c"] Pl not None,
+		 np.ndarray[double, ndim=1, mode="c"] Wll not None):
+    build_Wll( len(El), len(Pl[0]), 
+	       <double*> np.PyArray_DATA(El),
+	       <double*> np.PyArray_DATA(Pl),
+	       <double*> np.PyArray_DATA(Wll))
 
 def _dlss( X, s1, s2, lmax, np.ndarray[double, ndim=1, mode="c"] d not None):
      dlss( X, s1, s2, lmax, <double*> np.PyArray_DATA(d))

src/libcov.c

 #include "libcov.h"
 
-// void CrossWindowFunction( int nl, int npix, double* El, double* Pl, double* Wll)
-// {
 
-//   memset( Wll, 0., nl*nl * sizeof(double));
+//problem of precision wrt python...
+void build_Wll( int nl, int npix, double* El, double* Pl, double* Wll)
+{
+  int64_t npixtot = npix*npix;
+  
+  memset( Wll, 0., (nl*nl) * sizeof(double));
 
-// #pragma omp parallel default(none) shared(nl, El, Pl, Wll)
-//   {
-// #pragma omp for schedule(dynamic)
-//     for( int l1=0; l1<nl; l1++) {
-//       for( int l2=0; l2<nl; l2++) {
+#pragma omp parallel default(none) shared(nl, npixtot, El, Pl, Wll)
+  {
+#pragma omp for schedule(dynamic)
+    for( int l1=0; l1<nl; l1++) {
+      for( int l2=0; l2<nl; l2++) {
 	
-// 	for( int p=0; p<npix*npix; p++)
-// 	  Wll[l1*nl+l2] += El[l1][
-//       }
-//     }
+	for( int p=0; p<npixtot; p++)
+	  Wll[l1*nl+l2] += El[l1*npixtot+p]*Pl[l2*npixtot+ p];
+	
+      } //loop l2
+    } //loop l1
+    
+  } //end omp parallel
 
-//   }
+//   fprintf( stdout, "E[13] = ");
+//   for( int p=0; p<10; p++) fprintf( stdout, "%f\t", El[13*npixtot+p]);
+//   fprintf( stdout, "\n");
+//   fprintf( stdout, "P[4] = ");
+//   for( int p=0; p<10; p++) fprintf( stdout, "%f\t", Pl[4*npixtot+p]);
+//   fprintf( stdout, "\n");
+//   fprintf( stdout, "E[13]*P[4] = ");
+//   for( int p=0; p<10; p++) fprintf( stdout, "%f\t", El[13*npixtot+p]*Pl[4*npixtot+p]);
+//   fprintf( stdout, "\n");
+//   fprintf( stdout, "Wll[%d,%d]=%f\n", 13, 4, Wll[13*nl+4]);
   
-// }
+}
 
 
 void build_dSdC( int nside, int nstokes, int npix, int nbin, long *ispec, long *ellbins, long *ipix, double *bl, double* dSdC)

src/libcov.h

@@ -2,7 +2,6 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
-/* #include <chealpix.h> */
 
 #define OK 0
 #define NOK 1
@@ -12,6 +11,9 @@
 #define dSdC(L,I,J) dSdC[ (L)*(npixtot*npixtot) + (I)*npixtot + (J)]
 #define halfpi M_PI/2.
 
+void build_Wll( int nl, int npix, double* El, double* Pl, double* Wll);
+
+
 void build_dSdC( int nside, int nstokes, int npix, int inl, long *ispec, long *ellbins, long *ipix, double *bl, double* dSdC);
 /**
     Compute the Pl = dS/dCl matrices.

xqml/estimators.py

 """
-Set of routines to ...
-
-Author: Vanneste
+Set of routines generating the estimators for xQML code
 """
 from __future__ import division
 
 import numpy as np
+import _libcov as clibcov
+import timeit
+import threading
 
 
+#
 def Pl(ds_dcb):
     """
     Reshape ds_dcb (nspec, nbins) into Pl (nspec * nbins)
@@ -41,7 +43,7 @@ def Pl(ds_dcb):
     nnpix = np.shape(ds_dcb)[-1]
     return np.copy(ds_dcb).reshape(2 * (np.shape(ds_dcb)[1]), nnpix, nnpix)
 
-
+#
 def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, corr=False):
     """
     Compute correlation matrix S = sum_l Pl*Cl
@@ -101,9 +103,10 @@ def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, corr=False):
     return S
 
 
-def El(invCAA, invCBB, Pl):
+def El(invCAA, invCBB, Pl, thread=False, verbose=False):
     """
     Compute El = CAA^-1.Pl.CBB^-1
+    (Note: El is not symmetric for cross)
 
     Parameters
     ----------
@@ -136,16 +139,40 @@ def El(invCAA, invCBB, Pl):
 
     """
 
-    El = [np.dot(np.dot(invCAA, P), invCBB) for P in Pl]
+    tstart = timeit.default_timer()
+    nl = len(Pl)
+    npix = len(Pl[0])
+
+    if thread:
+        #Note: longer than with list...
+        El = np.ndarray( np.shape(Pl))
+        def CPC( l):
+            El[l] = np.dot(np.dot(invCAA, Pl[l]), invCBB)
+        procs = []
+        for l in range(nl):
+            proc = threading.Thread(target=CPC,args=(l,))
+            procs.append(proc)
+
+        for proc in procs:
+            proc.start()
+        for proc in procs:
+            proc.join()
+
+    else:
+        El = [np.dot(np.dot(invCAA, P), invCBB) for P in Pl]
+
+    if verbose:
+        print( "Construct El (nl=%d): %.1f sec" % (nl,timeit.default_timer()-tstart))
 
     return El
 
 
-def CrossWindowFunction(El, Pl):
+
+def CrossWindowFunction(El, Pl, openMP=False, thread=False, verbose=False):
     """
     Compute mode-mixing matrix (Tegmark's window matrix)
     Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
-    Use the trick with matrices: Trace[A.B] = sum(A x B) (where . is the matrix product and x the elementwise mult)
+    Use the trick with matrices: Trace[A.B] = sum(A x B^T) (where . is the matrix product and x the elementwise mult)
     
     Parameters
     ----------
@@ -169,13 +196,41 @@ def CrossWindowFunction(El, Pl):
      [134 478 822]]
     """
     nl = len(El)
-    
-    # No transpose because E symm
-    Wll = np.asarray( [np.sum(E * P) for E in El for P in Pl] ).reshape(nl,nl)
-    
+
+    tstart = timeit.default_timer()
+
+    if openMP:
+        #pb of precision (sum of +/- big numbers)
+        Wll = np.ndarray( nl*nl)
+        clibcov.CrossWindow( np.asarray(El), Pl, Wll)
+        Wll = Wll.reshape(nl,nl)
+    elif thread:
+        #gain a factor 2.5 on npix=600
+        Wll = np.ndarray( (nl, nl))
+        def EP( l1, l2):
+            Wll[l1,l2] = np.sum( El[l1]*Pl[l2])
+        procs = []
+        for l1 in range(nl):
+            for l2 in range(nl):                
+                proc = threading.Thread(target=EP, args=(l1,l2,))
+                procs.append(proc)
+
+        for proc in procs:
+                proc.start()
+        for proc in procs:
+                proc.join()
+        
+    else:
+        # No transpose because P symm
+        Wll = np.asarray( [np.sum(E * P) for E in El for P in Pl] ).reshape(nl,nl)
+
+    if verbose:
+        print( "Construct Wll (nl=%d): %.1f sec" % (nl,timeit.default_timer()-tstart))
+
     return Wll
 
 
+
 def CrossWindowFunctionLong(invCAA, invCBB, Pl):
     """
     Compute mode-mixing matrix (Tegmark's window matrix)
@@ -214,6 +269,7 @@ def CrossWindowFunctionLong(invCAA, invCBB, Pl):
     return Wll
 
 
+
 def CrossGisherMatrix(El, CAB):
     """
     Compute matrix GAB = Trace[El.CAB.El.CAB]
@@ -241,9 +297,10 @@ def CrossGisherMatrix(El, CAB):
     nl = len(El)
 
     El_CAB = [np.dot(CAB, E) for E in El]
-    GAB = np.asarray([np.sum(Ei * Ej.T) for Ei in El_CAB for Ej in El_CAB]).reshape(nl,nl)
+    GAB = [np.sum(Ei * Ej.T) for Ei in El_CAB for Ej in El_CAB]
     
-    return GAB
+    return np.asarray(GAB).reshape(nl,nl)
+
 
 
 def CrossGisherMatrixLong(El, CAB):
@@ -273,10 +330,8 @@ def CrossGisherMatrixLong(El, CAB):
     """
     lmax = len(El)
     lrange = np.arange(lmax)
-    GAB = np.asarray(
-        [np.sum(np.dot(CAB, El[il]) * np.dot(CAB, El[jl]).T)
-         for il in lrange for jl in lrange]).reshape(lmax, lmax)
-    return GAB
+    GAB = [np.sum(np.dot(CAB, El[il]) * np.dot(CAB, El[jl]).T) for il in lrange for jl in lrange]
+    return np.asarray(GAB).reshape(lmax, lmax)
 
 
 def yQuadEstimator(dA, dB, El):
@@ -291,7 +346,7 @@ def yQuadEstimator(dA, dB, El):
         Pixels dataset B
     El : ndarray of floats
         Quadratic parameter matrices such that yl = dA.El.dB.T
-
+    
     Returns
     ----------
     >>> dA = np.arange(12)
@@ -372,6 +427,7 @@ def CovAB(invWll, GAB):
     return covAB
 
 
+
 if __name__ == "__main__":
     """
     Run the doctest using

xqml/libcov.py

 """
 Set of routines to ...
-
-Author: Vanneste
 """
 from __future__ import division
 

xqml/xqml.py

 """
 (Cross-) power spectra estimation using the QML method.
 [Vanneste et al. 2018, arXiv:1807.02484]
-
-Author: Vanneste
 """
 from __future__ import division
 
@@ -10,7 +8,7 @@ import sys
 import timeit
 import string
 
-from scipy import linalg
+#from scipy import linalg
 
 import numpy as np
 import healpy as hp
@@ -41,7 +39,7 @@ class xQML(object):
         mask : 1D array of booleans
             Mask defining the region of interest (of value True)
         bins : 1D array of floats (nbin+1)
-            lower multipole bin
+            lower multipole for each bin
         clth : ndarray of floats
             Array containing fiducial CMB spectra (unbinned)
         lmax : int
@@ -128,7 +126,7 @@ class xQML(object):
                                           self.spec, pixwin=self.pixwin, verbose=verbose, MC=MC, openMP=openMP)
         return( self.Pl, self.S)
 
-    def construct_esti(self, NA, NB=None, verbose=False):
+    def construct_esti(self, NA, NB=None, verbose=False, thread=False):
         """
         Compute the inverse of the datasets pixel covariance matrices C,
         the quadratic matrix parameter E, and inverse of the window
@@ -154,11 +152,11 @@ class xQML(object):
         # Invert (signalB + noise) matrix
         invCb = pd_inv(self.S + self.NB)
 
-        # Compute El = Ca^-11.Pl.Cb^-1 (long)
-        self.El = El(invCa, invCb, self.Pl)
+        # Compute El = Ca^-1.Pl.Cb^-1 (long)
+        self.El = El(invCa, invCb, self.Pl, thread=thread, verbose=verbose)
         
-        # Finally compute invW by inverting (long)
-        self.invW = linalg.inv(CrossWindowFunction(self.El, self.Pl))
+        # Finally compute invW by inverting (longer)
+        self.invW = np.linalg.inv(CrossWindowFunction(self.El, self.Pl, thread=thread, verbose=verbose))
         
         # Compute bias for auto
 #        if not self.cross:

xqml/xqml_utils.py

 """
-Set of routines to ...
-
-Author: Vanneste
+Utils for xQML code
 """
 from __future__ import division
 
@@ -13,6 +11,10 @@ from scipy import linalg, sparse
 
 
 def pd_inv(a):
+    '''
+    Quicker to inverse for pixel-pixel matrices
+    (almost 5x quicker)
+    '''
     n = a.shape[0]
     I = np.identity(n)
     return linalg.solve(a, I, sym_pos = True, overwrite_b = True)
@@ -62,57 +64,6 @@ def getstokes(spec):
     return stokes, spec, istokes, ispecs
 
 
-def ComputeSizeDs_dcb(nside, fsky, deltal=1, unit="Gb"):
-    """
-    ???
-
-    Parameters
-    ----------
-    nside : ???
-        ???
-
-    Returns
-    ----------
-    ???
-    """
-    if units[0].lower() == "g":
-        toGB = 1024. * 1024. * 1024.
-    else:
-        toGB = 1.
-
-    nspec = 2
-    nell = nspec * (3.*nside-1) /deltal
-    npixtot = 2*12*nside**2*fsky
-    sizeds_dcb = (npixtot)**2 * (nell + 5)
-
-    return( 8.* sizeds_dcb / toGB)
-#    print("size (Gb) = " + str(sizeds_dcb))
-
-
-def get_colors(num_colors):
-    """
-    ???
-
-    Parameters
-    ----------
-    num_colors : ???
-        ???
-
-    Returns
-    ----------
-    ???
-    """
-    import colorsys
-    colors = []
-    ndeg = 250.
-    for i in np.arange(0., ndeg, ndeg / num_colors):
-        hue = i/360.
-        lightness = 0.5
-        saturation = 0.7
-        colors.append(colorsys.hls_to_rgb(hue, lightness, saturation))
-    return np.array(colors)
-
-
 def progress_bar(i, n):
     """
     ???
@@ -142,94 +93,8 @@ def progress_bar(i, n):
         sys.stdout.flush()
 
 
-def check_symmetric(a, tol=1e-8):
-    """
-    ???
-
-    Parameters
-    ----------
-    a : ???
-        ???
-
-    Returns
-    ----------
-    ???
-    """
-    return np.allclose(a, a.T, atol=tol)
-
-
-def randomword(length):
-    """
-    ???
-
-    Parameters
-    ----------
-    length : ???
-        ???
-
-    Returns
-    ----------
-    ???
-    """
-    return ''.join(rd.choice(string.lowercase) for i in range(length))
 
 
-def cov_from_maps(maps0, maps1):
-    """
-    ???
-
-    Parameters
-    ----------
-    maps0 : ???
-        ???
-
-    Returns
-    ----------
-    ???
-    """
-    sh = np.shape(maps0)
-    npix = sh[1]
-    nbmc = sh[0]
-    covmc = np.zeros((npix, npix))
-    mm0 = np.mean(maps0, axis=0)
-    # print(mm0)
-    mm1 = np.mean(maps1, axis=0)
-    # print(mm1)
-    themaps0 = np.zeros((nbmc, npix))
-    themaps1 = np.zeros((nbmc, npix))
-    start = timeit.default_timer()
-    for i in np.arange(npix):
-        progress_bar(i, npix, timeit.default_timer()-start)
-        themaps0[:, i] = maps0[:, i] - mm0[i]
-        themaps1[:, i] = maps1[:, i] - mm1[i]
-    print('hy')
-    start = timeit.default_timer()
-    for i in np.arange(npix):
-        progress_bar(i, npix, timeit.default_timer()-start)
-        for j in np.arange(npix):
-            covmc[i, j] = np.mean(themaps0[:, i] * themaps1[:, j])
-    return(covmc)
-
-
-def IsInvertible(F):
-    """
-    ???
-
-    Parameters
-    ----------
-    F : ???
-        ???
-    ...
-
-    Returns
-    ----------
-    ??? : ???
-        ???
-    """
-    eps = np.finfo(F.dtype).eps
-    print("Cond Numb = ", np.linalg.cond(F), "Matrix eps=", eps)
-    return np.linalg.cond(F) > np.finfo(F.dtype).eps
-
 
 
 def GetBinningMatrix( ellbins, lmax, norm=False):
@@ -342,12 +207,12 @@ def GetBinningMatrix( ellbins, lmax, norm=False):
 
 
 class symarray(np.ndarray):
-
-    # wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage.
-    # Usage:
-    # If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array 
-    # that is a packed version of A.  To convert it back, just wrap the flat list in Sym().  Note that Sym(Sym(A).packed)
-
+    '''
+    wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage.
+    Usage:
+    If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array 
+    that is a packed version of A.  To convert it back, just wrap the flat list in Sym().  Note that Sym(Sym(A).packed)
+    '''
     def __new__(cls, input_array):
         obj = np.asarray(input_array).view(cls)