Caps

lenspyx/__init__.py

@@ -1,3 +1,3 @@
-from lenspyx.lensing import alm2lenmap, alm2lenmap_spin, synfast, get_geom
+from lenspyx.lensing import alm2lenmap, alm2lenmap_spin, synfast, get_geom, dlm2angles
 from._version import __version__
 

lenspyx/cachers.py

@@ -24,7 +24,7 @@ def remove(self, fn):
 class cacher_npy(cacher):
     def __init__(self, lib_dir, verbose=False):
         if not os.path.exists(lib_dir):
-            os.makedirs(lib_dir)
+            os.makedirs(lib_dir, exist_ok=True)
         self.lib_dir = lib_dir
         self.verbose = verbose
 

lenspyx/experimental.py

@@ -0,0 +1,90 @@
+import numpy as np
+from ducc0.sht import synthesis_general as syngducc, adjoint_synthesis_general as adjsyngducc
+try:
+    import capsht
+except ImportError:
+    print("capsht not found, you will not be able to use the functions in this module")
+try:
+    from capsht.experimental import synthesis_general_cap, synthesis_general_band, adjoint_synthesis_general_cap, adjoint_synthesis_general_band
+except ImportError:
+    print("synthesis_general_cap or synthesis_general_band not found in capsht.experimental, are you up to date?")
+
+def _epsapo(thtcap, epsilon, lmax, version=1, dl_7=None):
+    #dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2
+    assert version == 1, 'C++ code now only implemented for version 1'
+    if version == 0:
+        if dl_7 is None:
+            dl_7 = 15
+        dl = dl_7 * ((-np.log10(epsilon) / (7 )) ** 1) ** 0.5
+    elif version == 1: 
+        if dl_7 is None:
+            dl_7 = 2*7*np.log(10.)/np.pi
+        dl = dl_7 * (-np.log10(epsilon) / (7. ))
+    else:
+        raise ValueError('version %s not implemented'%version)
+    return np.sqrt(dl / lmax * np.pi / thtcap)
+
+def synthesis_general(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, 
+                      thtcap:float=None, eps_apo:float=None, tht_min:float=None, tht_max:float=None, verbose:bool=False, **kwargs):
+    """Wrapper to capsht synthesis_general function, hiding the choice of eps_apo and SHT algorithm
+
+
+        See ducc0.sht.synthesis_general for arguments, optional arguments and outputs
+
+        relevant keyword, *mode*, *map*, *mmax*
+
+    """
+    if tht_min is not None and tht_max is not None: # attempt at synthesis_general_band
+        eps_apo = eps_apo or 1.2 * _epsapo(tht_max-tht_min, epsilon, lmax)
+        thta_p = tht_min - 0.5 * eps_apo * (tht_max - tht_min)
+        thtb_p = tht_max + 0.5 * eps_apo * (tht_max - tht_min)
+        if (thta_p >= 0.) and (thtb_p <= np.pi):
+            if verbose:
+                print('syng type: band %.1f deg %.1f deg' % (thta_p/np.pi*180, thtb_p/np.pi*180))
+            assert 0, 'fix band to new scheme'
+            return synthesis_general_band(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, 
+                                      thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs)
+        if thta_p < 0.: # Can try synthesis_general_cap later on
+            thtcap = tht_max
+            eps_apo = None
+    if thtcap is not None: # attempt at synthesis_general_cap
+        eps_apo = eps_apo or _epsapo(thtcap, epsilon, lmax)
+        epsilon_nufft = kwargs.pop('epsilon_nufft', epsilon)
+        if verbose:
+            print('syng type: sent to cap %.1f epsapo %.2f' % (thtcap/np.pi*180, eps_apo))
+        return synthesis_general_cap(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon_nufft, thtcap=thtcap, eps_apo=eps_apo, **kwargs)
+    if verbose:
+        print('syng type : general')
+    return syngducc(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon,  **kwargs)
+
+def adjoint_synthesis_general(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, 
+                      thtcap:float=None, eps_apo:float=None, tht_min:float=None, tht_max:float=None, verbose:bool=False, **kwargs):
+    """Wrapper to capsht synthesis_general function, hiding the choice of eps_apo
+
+
+        See ducc0.sht.synthesis_general for arguments, optional arguments and outputs
+
+        relevant keyword, *mode*, *alm*, *mmax*
+
+
+    """
+    if tht_min is not None and tht_max is not None: # attempt at synthesis_general_band
+        eps_apo = eps_apo or 1.2 * _epsapo(tht_max-tht_min, epsilon, lmax)
+        thta_p = tht_min - 0.5 * eps_apo * (tht_max - tht_min)
+        thtb_p = tht_max + 0.5 * eps_apo * (tht_max - tht_min)
+        if (thta_p >= 0.) and (thtb_p <= np.pi):
+            if verbose:
+                print('adjsyng type: band %.1f deg %.1f deg' % (thta_p/np.pi*180, thtb_p/np.pi*180))
+            return adjoint_synthesis_general_band(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, 
+                                      thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs)
+        if thta_p < 0.: # Can try synthesis_general_cap later on
+            thtcap = tht_max
+            eps_apo = None
+    if thtcap is not None: # attempt at synthesis_general_cap
+        eps_apo = eps_apo or _epsapo(thtcap, epsilon, lmax)
+        if verbose:
+            print('adjsyng type: sent to cap %.1f deg, eps_apo %.2f' % (thtcap/np.pi*180, eps_apo))
+        return adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs)
+    if verbose:
+        print('adjsyng type : general')
+    return syngducc(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon,  **kwargs)

lenspyx/qest/qest.py

@@ -14,7 +14,7 @@
 from lenspyx.qest.ivfs import OpFilt
 
 
-def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, verbose=False, get_alm2=None, geometry: Geom or None=None):
+def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, mmax=None, mmax_qlm=None, verbose=False, get_alm2=None, geometry: Geom or None=None):
     """Evaluates a quadratic estimator gradient and curl terms.
 
 
@@ -32,13 +32,14 @@ def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, verbose=False, get_
             glm and clm healpy arrays (gradient and curl terms of the QE estimate)
 
     """
+    assert mmax in [None, lmax_ivf], 'there is a bug with non-trivial mmmax, fix this first'
     qe_list = _get_qes(qe_key, lmax_ivf, cls_weight)
     if geometry is None:
         qe_spin = np.max([qe[0].spin_ou + qe[1].spin_ou for qe in uqe.qe_compress(qe_list)])
         geometry = Geom.get_thingauss_geometry((2 * lmax_ivf + lmax_qlm) // 2 + 1, qe_spin)
-    return _eval_qe(qe_list, get_alm, lmax_qlm, verbose=verbose, get_alm2=get_alm2, geo=geometry)
+    return _eval_qe(qe_list, get_alm, lmax_qlm, mmax_qlm=mmax_qlm, mmax=mmax, verbose=verbose, get_alm2=get_alm2, geo=geometry)
 
-def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, get_alm2=None, mmax_qlm:int or None=None, nthreads=0):
+def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, get_alm2=None, mmax:int or None=None, mmax_qlm:int or None=None, nthreads=0):
     """Evaluation of a QE from its list of leg definitions.
 
         Args:
@@ -105,25 +106,25 @@ def _eval_qe(qe_list:list[uqe.qe], get_alm, lmax_qlm, geo:Geom, verbose=True, ge
             if len(conjugate) > 0:
                 for j in conjugate:
                     print("in-spins conjugate leg and out-spin", qes[j][1].spins_in, qes[j][1].spin_ou)
-        a = q[0](get_alm, geo)
+        a = q[0](get_alm, geo, mmax=mmax)
         if qe_spin:
-            dc += fac1 * a * q[1](get_alm2, geo)
+            dc += fac1 * a * q[1](get_alm2, geo, mmax=mmax)
             for j in conjugate:
-                dc += a.conj() * qes[j][1](get_alm2, geo)
+                dc += a.conj() * qes[j][1](get_alm2, geo, mmax=mmax)
         else: # We must consider the real part only
-            dc += fac1 * (a * q[1](get_alm2, geo)).real
+            dc += fac1 * (a * q[1](get_alm2, geo, mmax=mmax)).real
             for j in conjugate:
-                dc += (a.conj() * qes[j][1](get_alm2, geo)).real
+                dc += (a.conj() * qes[j][1](get_alm2, geo, mmax=mmax)).real
         if symmetrize: # same, swapping alm2 and alm1
-            a = q[0](get_alm2, geo)
+            a = q[0](get_alm2, geo, mmax=mmax)
             if qe_spin:
-                dc += fac1 * a * q[1](get_alm, geo)
+                dc += fac1 * a * q[1](get_alm, geo, mmax=mmax)
                 for j in conjugate:
-                    dc += a.conj() * qes[j][1](get_alm, geo)
+                    dc += a.conj() * qes[j][1](get_alm, geo, mmax=mmax)
             else:
-                dc += fac1 * (a * q[1](get_alm, geo)).real
+                dc += fac1 * (a * q[1](get_alm, geo, mmax=mmax)).real
                 for j in conjugate:
-                    dc += (a.conj() * qes[j][1](get_alm, geo)).real
+                    dc += (a.conj() * qes[j][1](get_alm, geo, mmax=mmax)).real
 
     gclm = geo.adjoint_synthesis(m=dr, spin=qe_spin, lmax=lmax_qlm, mmax=mmax_qlm, nthreads=nthreads)
     if symmetrize:

lenspyx/qest/utils_qe.py

@@ -22,7 +22,7 @@ def __mul__(self, other):
     def __add__(self, other):
         assert self.spin_in == other.spin_in and self.spin_ou == other.spin_ou
         lmax = max(self.get_lmax(), other.get_lmax())
-        cl = np.zeros(lmax + 1, dtype=float)
+        cl = np.zeros(lmax + 1, dtype=self.cl.dtype)
         cl[:len(self.cl)] += self.cl
         cl[:len(other.cl)] += other.cl
         return qeleg(self.spin_in, self.spin_ou, cl)
@@ -50,7 +50,7 @@ def __iadd__(self, other_qe: qeleg):
         self.cls.append(np.copy(other_qe.cl))
         return self
 
-    def __call__(self, get_alm: callable, geometry: Geom, nthreads: int = 0):
+    def __call__(self, get_alm: callable, geometry: Geom, mmax=None, nthreads: int = 0):
         """Returns the spin-weighted real-space map of the estimator.
 
             We first build X_lm in the wanted _{si}X_lm _{so}Y_lm and then convert this alm2map_spin conventions.
@@ -59,10 +59,10 @@ def __call__(self, get_alm: callable, geometry: Geom, nthreads: int = 0):
         if nthreads <= 0:
             nthreads = cpu_count()
         lmax = self.get_lmax()
-        mmax = lmax
+        if mmax is None:
+            mmax = lmax
         ncomp, npix = 1 + (self.spin_ou != 0), geometry.npix()
-        alm_size = Alm.getsize(lmax, mmax)
-        gclm = np.zeros((ncomp, alm_size), dtype=complex)
+        gclm = np.zeros((ncomp, Alm.getsize(lmax, mmax)), dtype=complex)
         for i, (si, cl) in enumerate(zip(self.spins_in, self.cls)):
             assert si in [0, -2, 2], str(si) + ' input spin not implemented'
             alms = [get_alm('e'), get_alm('b')] if abs(si) == 2 else [-get_alm('t'), 0]

lenspyx/remapping/utils_geom.py

@@ -4,7 +4,7 @@
 import ducc0
 from ducc0.misc import GL_thetas, GL_weights
 from ducc0.fft import good_size
-from ducc0.sht.experimental import synthesis, adjoint_synthesis, synthesis_deriv1
+from ducc0.sht import synthesis, adjoint_synthesis, synthesis_deriv1, synthesis_general, adjoint_synthesis_general
 
 def st2mmax(spin, tht, lmax):
     r"""Converts spin, tht and lmax to a maximum effective m, according to libsharp paper polar optimization formula Eqs. 7-8
@@ -49,6 +49,14 @@ def npix(self):
         """
         return int(np.sum(self.nph))
 
+    def mmax(self, spin:int, lmax:int):
+        """Safe value of mmax for a given lmax and spin, given the latitude range of the geometry
+
+
+        """
+        return min(int(st2mmax(spin, np.pi * 0.5 - np.min(np.abs(np.pi * 0.5 - self.theta)), lmax)) + 1, lmax)
+
+
     def fsky(self):
         """Fractional area of the sky covered by the pixelization
 
@@ -175,6 +183,7 @@ def adjoint_synthesis(self, m: np.ndarray, spin:int, lmax:int, mmax:int, nthread
                                  nthreads=nthreads, ringstart=self.ofs, alm=alm,  **kwargs)
 
     def alm2map_spin(self, gclm:np.ndarray, spin:int, lmax:int, mmax:int, nthreads:int, zbounds=(-1., 1.), **kwargs):
+        # FIXME: method only here for backwards compatiblity
         # FIXME: method only here for backwards compatiblity
         assert zbounds[0] == -1 and zbounds[1] == 1., zbounds
         return self.synthesis(gclm, spin, lmax, mmax, nthreads, **kwargs)
@@ -342,7 +351,7 @@ def get_tgl_geometry(lmax:int, smax:int, good_size_real=True):
 
         """
         return Geom.get_thingauss_geometry(lmax, smax, good_size_real=good_size_real)
-
+    
 
 class pbounds:
     """Class to regroup simple functions handling sky maps longitude truncation

lenspyx/tests/helper.py

@@ -21,14 +21,19 @@ def _extend_cl(cl:np.ndarray, lmax):
         ret[:lmax_cl+1] = cl
         return ret
 
-def syn_alms(spin, lmax_unl=5120, ctyp=np.complex128):
+def syn_alms(spin, lmax_unl=5120, ctyp=np.complex128, white=False):
     ncomp = 1 + (abs(spin) > 0)
     mmax_unl = lmax_unl
     rtyp = lenspyx.remapping.deflection_028.rtype[ctyp]
     eblm = np.empty( (ncomp, Alm.getsize(lmax_unl, mmax_unl)), dtype=ctyp)
-    eblm[0] = synalm(_extend_cl(cls_unl['ee' if abs(spin) > 0 else 'tt'][:lmax_unl + 1], lmax_unl), lmax_unl, mmax_unl, rlm_dtype=rtyp)
+    clgg = _extend_cl(cls_unl['ee' if abs(spin) > 0 else 'tt'][:lmax_unl + 1], lmax_unl)
+    clcc = _extend_cl(cls_unl['bb'][:lmax_unl + 1], lmax_unl)
+    if white:
+        clgg = np.ones_like(clgg)
+        clcc = np.ones_like(clcc)
+    eblm[0] = synalm(clgg, lmax_unl, mmax_unl, rlm_dtype=rtyp)
     if ncomp > 1:
-       eblm[1] = synalm(_extend_cl(cls_unl['bb'][:lmax_unl + 1], lmax_unl), lmax_unl, mmax_unl, rlm_dtype=rtyp)
+       eblm[1] = synalm(clcc, lmax_unl, mmax_unl, rlm_dtype=rtyp)
     return eblm
 
 def syn_dlm(lmax_unl=5120, ctyp=np.complex128):

lenspyx/tests/test_adjsyng_cap.py

@@ -0,0 +1,64 @@
+import numpy as np
+from lenspyx.tests.helper import syn_ffi_ducc_29
+import pylab as pl
+from duccjc.sht import adjoint_synthesis_general_ringweight as adj_syng_w
+from ducc0.sht import adjoint_synthesis_general as adj_syng
+from lenspyx.utils_hp import alm_copy, alm2cl
+from lenspyx.remapping.utils_geom import st2mmax
+from lenspyx.utils_cap import fskycap, eps_opti, args_default
+
+if __name__ == '__main__':
+    """This tests adjoint_synthesis_general_cap
+
+        This generates delfected positions according to LCDM,and compare output alm of full-sky vs cap routines
+
+    """
+    args = args_default()
+    args.whiten = True
+
+    thta = 0 / 180 * np.pi
+    thtb = 50 / 180 * np.pi # cap size
+    #dtheta = 7. / 180 * np.pi
+
+
+    ffi, geom = syn_ffi_ducc_29(lmax_len=args.lmax_len, dlmax=args.dlmax, dlmax_gl=args.dlmax_gl, nthreads=args.nt,
+                             verbosity=1, epsilon=10 ** (-args.epsilon))
+    lmax_unl, mmax_unl = args.lmax_len + args.dlmax, args.lmax_len + args.dlmax
+
+    sht_mode = 'STANDARD' if (args.spin == 0 or not args.gonly) else 'GRAD ONLY'
+    ncomp = 1 + (args.spin > 0) * (not args.gonly)
+    ptg = ffi._get_ptg() # locations
+    assert ptg.shape[-1] == 2, ptg.shape
+
+
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=False, northsouth_sym=False)
+    slic = slice(np.min(geom_trc.ofs), np.max(geom_trc.ofs + geom_trc.nph))
+    ptg_sliced = ptg[slic, :]
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=True, northsouth_sym=False)
+    m = np.random.standard_normal((ncomp, ptg_sliced.shape[0]))
+    # Full-sky adjoint synthesis_general, at higher accuracy
+    eblm = adj_syng(map=m, lmax=lmax_unl, mmax=mmax_unl, loc=ptg_sliced, spin=args.spin, epsilon=ffi.epsilon * 0.1,
+                  nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity * 1)
+
+    if thta == 0:    # capped synthesis_general
+        thtcap = np.max(ptg_sliced[:, 0]) * 1.0001
+        eps_apo = eps_opti(lmax_unl, thtcap, dl=args.dl)
+        new_mmax = min(int(st2mmax(args.spin, thtcap * (1. + eps_apo), lmax_unl)) + 1, mmax_unl) if args.adapt_mmax else mmax_unl
+        print("Testing capped syng %2.f deg, fsky %.2f"%( (thtcap / np.pi * 180), fskycap(thtcap)))
+        print('tht cap in deg %.0f, eps apo %.2f'%(thtcap/np.pi * 180, eps_apo))
+        print('induced dlmax %s'%(int(lmax_unl *eps_apo)))
+        print('reduction in mmax by a factor %.1f'%(mmax_unl/new_mmax))
+
+        eblm_2= adj_syng_w(map=m, lmax=lmax_unl, mmax=new_mmax, loc=ptg_sliced,
+                spin=args.spin, epsilon=ffi.epsilon,
+                nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity, thtcap=thtcap, eps_apo=eps_apo, apofct=args.apofct)
+        # Resize alms to full-sky adjoint synthesis general
+        eblm_2 = np.array([alm_copy(alm, new_mmax, lmax_unl, mmax_unl) for alm in eblm_2])
+        for i, (ref, diff) in enumerate(zip(eblm, eblm_2 - eblm)):
+            cldiff = alm2cl(diff, diff, lmax_unl, mmax_unl, lmax_unl)
+            clref = alm2cl(ref, ref, lmax_unl, mmax_unl, lmax_unl)
+            pl.semilogy(np.sqrt(cldiff[2:]/clref[2:]), label='comp %s capped, weighted (dl=%s)'%(i, args.dl))
+            pl.xlabel(r'$\ell$')
+            pl.ylabel(r'$C_\ell$')
+    pl.legend()
+    pl.show()

lenspyx/tests/test_syng_cap.py

@@ -0,0 +1,63 @@
+import numpy as np
+from lenspyx.tests.helper import syn_ffi_ducc_29, syn_alms
+
+import pylab as pl
+from duccjc.sht import synthesis_general_ringweight as syng_w, adjoint_synthesis_general_ringweight as adj_syng_w
+from ducc0.sht import synthesis_general as syng, adjoint_synthesis_general as adj_syng
+from lenspyx.utils_hp import alm_copy
+from lenspyx.remapping import deflection_029, deflection_028 as deflection_28
+from lenspyx.remapping.utils_geom import st2mmax
+from lenspyx.utils_cap import fskycap, eps_opti, examine, args_default
+
+if __name__ == '__main__':
+    """This tests synthesis_general_cap
+
+        This generates delfected positions according to LCDM,and compare output maps of full-sky vs cap routines
+
+    """
+    args = args_default()
+    thta = 0 / 180 * np.pi
+    thtb = 50 / 180 * np.pi # cap size
+#    args.whiten = False
+
+
+    ffi, geom = syn_ffi_ducc_29(lmax_len=args.lmax_len, dlmax=args.dlmax, dlmax_gl=args.dlmax_gl, nthreads=args.nt,
+                             verbosity=1, epsilon=10 ** (-args.epsilon))
+    lmax_unl, mmax_unl = args.lmax_len + args.dlmax, args.lmax_len + args.dlmax
+
+    eblm = syn_alms(args.spin, lmax_unl=lmax_unl, ctyp=np.complex64 if ffi.single_prec else np.complex128, white=args.whiten)
+    eblm = np.atleast_2d(eblm)
+    if args.gonly:
+        eblm = eblm[:1]
+    sht_mode = deflection_28.ducc_sht_mode(eblm, args.spin)
+    ptg = ffi._get_ptg() # locations
+    assert ptg.shape[-1] == 2, ptg.shape
+
+
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=False, northsouth_sym=False)
+    slic = slice(np.min(geom_trc.ofs), np.max(geom_trc.ofs + geom_trc.nph))
+    ptg_sliced = ptg[slic, :]
+    geom_trc = ffi.geom.restrict(thta, thtb, update_ringstart=True, northsouth_sym=False)
+
+    # Full-sky synthesis_general, at higher accuracy
+    values = syng(lmax=lmax_unl, mmax=mmax_unl, alm=eblm, loc=ptg, spin=args.spin, epsilon=ffi.epsilon * 0.1,
+                  nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity * 1)
+
+    if thta == 0:    # capped synthesis_general
+        thtcap = np.max(ptg_sliced[:, 0]) * 1.0001
+        eps_apo = eps_opti(lmax_unl, thtcap, dl=args.dl)
+        new_mmax = min(int(st2mmax(args.spin, thtcap * (1. + eps_apo), lmax_unl)) + 1, mmax_unl) if args.adapt_mmax else mmax_unl
+        print("Testing capped syng %2.f deg, fsky %.2f"%( (thtcap / np.pi * 180), fskycap(thtcap)))
+        print('tht cap in deg %.0f, eps apo %.2f'%(thtcap/np.pi * 180, eps_apo))
+        print('induced dlmax %s'%(int(lmax_unl *eps_apo)))
+        print('reduction in mmax by a factor %.1f'%(mmax_unl/new_mmax))
+
+        gclm = np.array([alm_copy(alm, mmax_unl, lmax_unl, new_mmax) for alm in eblm])
+        values_w = syng_w(lmax=lmax_unl, mmax=new_mmax, alm=gclm, loc=ptg_sliced,
+                ringweights=np.array([1.]), spin=args.spin, epsilon=ffi.epsilon,
+                nthreads=ffi.sht_tr, mode=sht_mode, verbose=ffi.verbosity, thtcap=thtcap, eps_apo=eps_apo, apofct=args.apofct)
+        print(np.max(np.abs(values_w - values[:, slic])))
+        for i, (val, diff) in enumerate(zip(values[:, slic], values_w-values[:, slic])):
+            examine(val, diff, geom_trc, ffi.epsilon, label='comp %s capped, weighted (dl=%s)'%(i, args.dl))
+    pl.legend()
+    pl.show()