From e827c506345cc08298768ba610510980ed0fc2b6 Mon Sep 17 00:00:00 2001 From: jcarron Date: Mon, 17 Feb 2025 14:26:29 +0100 Subject: [PATCH 01/18] minor --- lenspyx/tests/helper.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/lenspyx/tests/helper.py b/lenspyx/tests/helper.py index 5eb80d0..6b38751 100644 --- a/lenspyx/tests/helper.py +++ b/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): From 7b255715f215b7168eb1ca1dcb769cae5adf36de Mon Sep 17 00:00:00 2001 From: jcarron Date: Fri, 21 Feb 2025 16:10:23 +0100 Subject: [PATCH 02/18] misc --- lenspyx/remapping/utils_geom.py | 1 + lenspyx/tests/test_adjsyng_cap.py | 64 +++++++++++++++ lenspyx/tests/test_syng_cap.py | 63 +++++++++++++++ lenspyx/utils_cap.py | 126 ++++++++++++++++++++++++++++++ lenspyx/utils_hp.py | 10 ++- 5 files changed, 260 insertions(+), 4 deletions(-) create mode 100644 lenspyx/tests/test_adjsyng_cap.py create mode 100644 lenspyx/tests/test_syng_cap.py create mode 100644 lenspyx/utils_cap.py diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py index bb6ad1f..1413236 100644 --- a/lenspyx/remapping/utils_geom.py +++ b/lenspyx/remapping/utils_geom.py @@ -176,6 +176,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) diff --git a/lenspyx/tests/test_adjsyng_cap.py b/lenspyx/tests/test_adjsyng_cap.py new file mode 100644 index 0000000..c4eb278 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/lenspyx/tests/test_syng_cap.py b/lenspyx/tests/test_syng_cap.py new file mode 100644 index 0000000..231c1a5 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/lenspyx/utils_cap.py b/lenspyx/utils_cap.py new file mode 100644 index 0000000..25b4a74 --- /dev/null +++ b/lenspyx/utils_cap.py @@ -0,0 +1,126 @@ +import numpy as np +from lenspyx.remapping.utils_geom import Geom +import pylab as pl + +class args_default: + def __init__(self): + self.lmax_len = 4000 + self.dlmax_gl = 500 + self.spin = 2 + self.nt = 4 + self.HL = 0 + self.alloc = 0 + self.tracemalloc = False + self.epsilon = 7 + self.gonly = False + self.dlmax = 500 + self.whiten = True # make the spectra white before interpolation + self.apofct = 0 + self.dl=7 + self.adapt_mmax = True + + +def fskycap(thtcap): + return (1. - np.cos(thtcap)) * 0.5 + + +class transition01: + """Helpers for transition functions equal to 1 at 0 and 0 at 1""" + @staticmethod + def _eval01(x, version): + assert np.all( (x < 1) & (x > 0)) + if version in ['ES', 'es', 3]: + # Exponential of semi-circle, with b giving 1/2 at 1/2 + beta = np.log(0.5) / (np.sqrt(3/4.) - 1.) + return np.exp(beta * (np.sqrt(1. - x ** 2) - 1)) + if version in ['KB', 'kb', 4]: # Kaiser-bessel + b = 5.74 + return np.i0(b * np.sqrt(1. - x ** 2)) / np.i0(b) + if version in [0]: # Smooth typical bump function. this is 1/2 at 1/2 + fx = np.exp(-1. / x) + f1x = np.exp(-1. / (1 - x)) + return f1x / (fx + f1x) + if version in ['Hann', 1]: + return np.cos(x * np.pi * 0.5) ** 2 + if version in [2]: # another smooth bump + return np.exp(-1. /( 1 - x ** 2) + 1.) + + @staticmethod + def eval(x, version): + ret = np.zeros(x.size) + ret[np.where(x <= 0.)] = 1. + ret[np.where(x >= 1.)] = 0. + i = np.where( (x > 0) & (x < 1)) + ret[i] = transition01._eval01(x[i], version) + return ret + + @staticmethod + def eval_fft(npts, version): + x = np.arange(npts) * ( (2 * np.pi) / npts ) - np.pi + return np.fft.fft(transition01.eval(x, version)) + + + @staticmethod + def bump(x, x1, dx, version): + """bump, equal to 1 on [-1, 1], and 0 outside of [-1-dx, 1+dx] + + """ + assert dx >= 0, dx + ret = np.zeros(x.size) + ax = np.abs(x) + ret[np.where(ax <= x1)] = 1. + ret[np.where(ax >= (x1 + dx))] = 0. + i = np.where( (ax > x1) & (ax < (x1 + dx)) ) + ret[i] = transition01.eval( (ax[i] - x1) / dx, version) + return ret + + @staticmethod + def bump_fft(npts, x1, dx, version): + x = np.arange(npts) * ( (2 * np.pi) / npts ) - np.pi + bp = transition01.bump(x, x1, dx, version) + return np.fft.fft(bp) + +def Nc(lmax, dl, thetacap, dtheta): + """ + + Args: + lmax: band-limit of the alm array + dl: Effective band-limit of the apodization window (for white spectra, dl=7 gives single precision results?) + thetacap: co-latitude angle defining the cap + dtheta: angular distance dedicated to apodization + + This should have a mimimum at sqrt(dl / lmax * (thetcap /pi)) + + Returns: + + """ + return (lmax + np.pi * dl / dtheta) * (thetacap + dtheta) / np.pi + +# parametrize interval +def eps_opti(lmax, thetacap, dl=7): + """Guess of optimum choice of apodization length, given the alm band-limit lmax and coordinate of the cap, and the [0, pi) window bandlimit""" + guess = np.sqrt(dl / lmax * np.pi / thetacap * np.pi) + return min(guess, (np.pi / thetacap - 1.)*(1.-1e-13)) # Cant overshoot + +def examine(ref, diff, geom:Geom, epsilon_ref, label=''): + """Look at difference to the ref map ring by ring + + + """ + rad2deg = 180 / np.pi + # extracts rings and plot rms dev + rms = np.zeros(geom.theta.size) + thta = np.min(geom.theta) + thtb = np.max(geom.theta) + dtheta = 0.1 * (thtb - thta) + for i, ir in enumerate(np.argsort(geom.theta)): + pix = geom.rings2pix(geom, [ir]) + rms[i] = np.sqrt(np.mean(diff[pix] ** 2)) /np.sqrt(np.mean(ref[pix] ** 2)) + pl.semilogy(np.sort(geom.theta) * rad2deg, rms, label=label) + pl.axhline(epsilon_ref,c='k') + pl.axvline(thta * rad2deg, c='grey') + pl.axvline(thtb * rad2deg, c='grey') + pl.xlim( max((thta-1.2 * dtheta) * rad2deg, 0.), min((thtb + 1.2 * dtheta) * rad2deg, np.pi * rad2deg)) + pl.xlabel(r'$\theta$ [deg]') + pl.ylabel(r'rel dev. (rms)') + return rms \ No newline at end of file diff --git a/lenspyx/utils_hp.py b/lenspyx/utils_hp.py index d751844..02ccddd 100644 --- a/lenspyx/utils_hp.py +++ b/lenspyx/utils_hp.py @@ -94,10 +94,10 @@ def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64): return alm -def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = np.float64): +def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np.float64): """Creates Gaussian field alms from input cl dictionary - Parametersseed + Parameters ---------- cls : dict The power spectra of the maps (e.g. as coming from CAMB) @@ -105,7 +105,9 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n Maximum multipole simulated mmax: int Maximum m defining the alm layout, defaults to lmax if None or < 0 - rlm_dtype(optional, defaults to np.float64): + seed: (optional, defaults to None) + Random generator seed + rlm_dtype:(optional, defaults to np.float64) Precision of real components of the array (e.g. np.float32 for single precision output array) Returns @@ -174,7 +176,7 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type = n elif rlm_dtype == np.float64: dtype_complex = np.complex128 else: - assert 0, "please either choose np.float32 (single precission), or np.float64 (double precission) as rlm_dtype" + assert 0, "please either choose np.float32 (single precision), or np.float64 (double precision) as rlm_dtype" alms = np.zeros((len(labels_wgrad), phases[0].size), dtype=dtype_complex) # for L in Ls: #L @ L.T is full matrx for i, f in enumerate(labels): From c39836f5d8d4ebfbf40d1e25f4d23b1f9e1b25dc Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 25 Feb 2025 15:35:34 +0100 Subject: [PATCH 03/18] dlm2angles and minors --- lenspyx/lensing.py | 42 ++++++++++++++++++++++++++++++++- lenspyx/remapping/utils_geom.py | 1 - lenspyx/wigners/wigners.py | 2 +- 3 files changed, 42 insertions(+), 3 deletions(-) diff --git a/lenspyx/lensing.py b/lenspyx/lensing.py index d850323..93465d6 100644 --- a/lenspyx/lensing.py +++ b/lenspyx/lensing.py @@ -1,12 +1,15 @@ from __future__ import print_function, annotations from os import cpu_count import numpy as np +from numpy.random import default_rng + +from ducc0.misc import get_deflected_angles + from lenspyx.remapping.utils_geom import Geom from lenspyx.remapping.deflection_029 import deflection from lenspyx import cachers from lenspyx.utils_hp import almxfl, Alm from lenspyx.utils import timer -from numpy.random import default_rng def get_geom(geometry: tuple[str, dict]=('healpix', {'nside':2048})): @@ -22,6 +25,43 @@ def get_geom(geometry: tuple[str, dict]=('healpix', {'nside':2048})): return geo(**geometry[1]) +def dlm2angles(dlms:np.ndarray, geometry:Geom, mmax=None, nthreads: int=0, calc_rotation=False): + r"""Returns pointing information from lensing deflection harmonic coefficients + + Args: + dlms: The spin-1 deflection, in the form of one or two arrays. + + The two arrays are the gradient and curl deflection healpy alms: + + :math:`\sqrt{L(L+1)}\phi_{LM}` with :math:`\phi` the lensing potential + + :math:`\sqrt{L(L+1)}\Omega_{LM}` with :math:`\Omega` the lensing curl potential + + + The curl can be omitted if zero, resulting in principle in slightly faster execution + + mmax(optional): maximum m value of dlms, if not equal to lmax + geometry(optional): desired sphere pixelization (here only iso-latitude rings) + nthreads(optional): number of threads to use (defaults to os.cpu_count()) + calc_rotation(optional): also computes the angle gamma by which to rotate non-zero spin fields after deflection + e.g. :math:`{2}_{P} \rightarrow e^{2 i \gamma}{}_{2}P` + + Returns: + array of shape (npix, 2 (or 3 if calc_rotation is set)) with co-latitude, longitude and rotation angle + :math:`\theta, \phi, \gamma` + + + """ + if nthreads <= 0: + nthreads = cpu_count() + dlms2d = np.atleast_2d(dlms) + lmax = Alm.getlmax(dlms2d[0].size, mmax=mmax) + assert dlms[0].size == Alm.getsize(lmax, mmax), ('Inconsistent input lmax and mmax', (lmax, mmax)) + tht, phi0, nph, ofs = geometry.theta, geometry.phi0, geometry.nph, geometry.ofs + d1 = geometry.synthesis(dlms2d, 1, lmax, mmax, nthreads) + return get_deflected_angles(theta=tht, phi0=phi0, nphi=nph, ringstart=ofs, deflect=d1.T, calc_rotation=calc_rotation, nthreads=nthreads) + + def alm2lenmap(alm, dlms, geometry: tuple[str, dict]=('healpix', {'nside':2048}), epsilon=1e-7, verbose=0, nthreads: int=0, pol=True): r"""Computes lensed CMB maps from their alm's and deflection field alm's. diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py index 1413236..a8578d4 100644 --- a/lenspyx/remapping/utils_geom.py +++ b/lenspyx/remapping/utils_geom.py @@ -6,7 +6,6 @@ from ducc0.fft import good_size from ducc0.sht.experimental import synthesis, adjoint_synthesis, synthesis_deriv1 - 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 diff --git a/lenspyx/wigners/wigners.py b/lenspyx/wigners/wigners.py index 9cc83db..f4bb51c 100644 --- a/lenspyx/wigners/wigners.py +++ b/lenspyx/wigners/wigners.py @@ -7,7 +7,7 @@ """ from __future__ import annotations import numpy as np -from ducc0.sht.experimental import alm2leg, leg2alm +from ducc0.sht import alm2leg, leg2alm from ducc0.misc import GL_thetas, GL_weights GL_cache = {} From 52135ebf4eadcb08e361b88a146c6c815a1f765d Mon Sep 17 00:00:00 2001 From: jcarron Date: Fri, 21 Feb 2025 16:22:36 +0100 Subject: [PATCH 04/18] minor --- lenspyx/remapping/deflection_028.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lenspyx/remapping/deflection_028.py b/lenspyx/remapping/deflection_028.py index 59fcffb..ee41b08 100644 --- a/lenspyx/remapping/deflection_028.py +++ b/lenspyx/remapping/deflection_028.py @@ -65,6 +65,9 @@ def __init__(self, lens_geom:Geom, dglm, mmax_dlm:int or None, numthreads:int=0, """ + assert np.iscomplexobj(dglm), 'harmonic coefficients must be complex' + if dclm is not None: + assert np.iscomplexobj(dclm), 'harmonic coefficients must be complex' lmax = Alm.getlmax(dglm.size, mmax_dlm) if mmax_dlm is None: mmax_dlm = lmax From 52b82150f6ef1d58466f12451265fa60e32794a4 Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 25 Feb 2025 15:44:27 +0100 Subject: [PATCH 05/18] doc --- docs/lensing.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/lensing.rst b/docs/lensing.rst index 3db3730..bd11687 100644 --- a/docs/lensing.rst +++ b/docs/lensing.rst @@ -5,4 +5,4 @@ lenspyx.lensing See examples/demo_lenspyx.ipynb for example calls of this couple of routines .. automodule:: lenspyx.lensing - :members: alm2lenmap, alm2lenmap_spin, synfast + :members: alm2lenmap, alm2lenmap_spin, synfast, dlm2angles From 154c608ccc2058dd90afe455fc3ec9c40850b78e Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 25 Feb 2025 15:47:14 +0100 Subject: [PATCH 06/18] doc --- lenspyx/lensing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenspyx/lensing.py b/lenspyx/lensing.py index 93465d6..422ad39 100644 --- a/lenspyx/lensing.py +++ b/lenspyx/lensing.py @@ -44,7 +44,7 @@ def dlm2angles(dlms:np.ndarray, geometry:Geom, mmax=None, nthreads: int=0, calc_ geometry(optional): desired sphere pixelization (here only iso-latitude rings) nthreads(optional): number of threads to use (defaults to os.cpu_count()) calc_rotation(optional): also computes the angle gamma by which to rotate non-zero spin fields after deflection - e.g. :math:`{2}_{P} \rightarrow e^{2 i \gamma}{}_{2}P` + e.g. :math:`{}_{2}{P}(\hat n) \rightarrow e^{2 i \gamma(\hat n)}{}_{2}P(\hat n')` Returns: array of shape (npix, 2 (or 3 if calc_rotation is set)) with co-latitude, longitude and rotation angle From a014b688b54002f7ffeb7417fa5c0a37c188d08b Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 25 Feb 2025 15:48:40 +0100 Subject: [PATCH 07/18] doc --- lenspyx/lensing.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lenspyx/lensing.py b/lenspyx/lensing.py index 422ad39..86023e8 100644 --- a/lenspyx/lensing.py +++ b/lenspyx/lensing.py @@ -47,8 +47,7 @@ def dlm2angles(dlms:np.ndarray, geometry:Geom, mmax=None, nthreads: int=0, calc_ e.g. :math:`{}_{2}{P}(\hat n) \rightarrow e^{2 i \gamma(\hat n)}{}_{2}P(\hat n')` Returns: - array of shape (npix, 2 (or 3 if calc_rotation is set)) with co-latitude, longitude and rotation angle - :math:`\theta, \phi, \gamma` + angles: array of shape (npix, 2 (or 3 if calc_rotation is set)) with co-latitude, longitude and rotation angles :math:`\theta, \phi, \gamma` """ From dc2f50f07bda6b89c15248eaace4de7e4062f404 Mon Sep 17 00:00:00 2001 From: jcarron Date: Mon, 10 Mar 2025 14:33:04 +0100 Subject: [PATCH 08/18] minor --- lenspyx/remapping/utils_geom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py index a8578d4..056a948 100644 --- a/lenspyx/remapping/utils_geom.py +++ b/lenspyx/remapping/utils_geom.py @@ -197,7 +197,7 @@ def map2alm(self, m:np.ndarray, lmax:int, mmax:int, nthreads:int, zbounds=(-1., @staticmethod def rings2pix(geom:Geom, rings:np.ndarray[int]): - return np.concatenate([geom.ofs[ir] + np.arange(geom.nph[ir], dtype=int) for ir in rings]) + return np.concatenate([geom.ofs[ir] + np.arange(geom.nph[ir], dtype=np.uint64) for ir in rings]) @staticmethod def phis(geom:Geom, ir): From 8aa9b6a4f5a335aa2ae0c69e3abbc5863d0ea093 Mon Sep 17 00:00:00 2001 From: jcarron Date: Wed, 12 Mar 2025 08:52:27 +0100 Subject: [PATCH 09/18] minor --- lenspyx/lensing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenspyx/lensing.py b/lenspyx/lensing.py index 86023e8..c56b3e3 100644 --- a/lenspyx/lensing.py +++ b/lenspyx/lensing.py @@ -57,7 +57,7 @@ def dlm2angles(dlms:np.ndarray, geometry:Geom, mmax=None, nthreads: int=0, calc_ lmax = Alm.getlmax(dlms2d[0].size, mmax=mmax) assert dlms[0].size == Alm.getsize(lmax, mmax), ('Inconsistent input lmax and mmax', (lmax, mmax)) tht, phi0, nph, ofs = geometry.theta, geometry.phi0, geometry.nph, geometry.ofs - d1 = geometry.synthesis(dlms2d, 1, lmax, mmax, nthreads) + d1 = geometry.synthesis(dlms2d, 1, lmax, mmax, nthreads, mode='STANDARD' if dlms2d.shape[0] == 2 else 'GRAD_ONLY') return get_deflected_angles(theta=tht, phi0=phi0, nphi=nph, ringstart=ofs, deflect=d1.T, calc_rotation=calc_rotation, nthreads=nthreads) From 3df10665e7ab33e019ad02da84bde8cd5b564265 Mon Sep 17 00:00:00 2001 From: jcarron Date: Fri, 28 Mar 2025 14:55:21 -0700 Subject: [PATCH 10/18] experimental syng_cap --- lenspyx/__init__.py | 2 +- lenspyx/remapping/utils_geom.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lenspyx/__init__.py b/lenspyx/__init__.py index 5eebc5d..cff7c8d 100644 --- a/lenspyx/__init__.py +++ b/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__ diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py index 056a948..adb9fb9 100644 --- a/lenspyx/remapping/utils_geom.py +++ b/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 From 4961ecb02a838817bdc30cbfb953e69cfc3e277e Mon Sep 17 00:00:00 2001 From: jcarron Date: Mon, 31 Mar 2025 17:52:24 +0200 Subject: [PATCH 11/18] adj_syng_cap --- lenspyx/experimental.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 lenspyx/experimental.py diff --git a/lenspyx/experimental.py b/lenspyx/experimental.py new file mode 100644 index 0000000..394bafe --- /dev/null +++ b/lenspyx/experimental.py @@ -0,0 +1,41 @@ +import numpy as np +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 +except ImportError: + print("synthesis_general_cap not found in capsht.experimental, are you up to date?") + +def _epsapo(thtcap, epsilon, lmax, dl_7=20): + dl = int(np.round(dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2)) + return np.sqrt(dl / lmax * np.pi / thtcap), + +def synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): + """Wrapper to capsht synthesis_general_cap function, hiding the choice of eps_apo + + If thtcap is not specified, it is set to np.max(loc[:, 0]) + + See ducc0.sht.synthesis_general for arguments, optional arguments and outputs + + """ + thtcap = kwargs.pop('thtcap', None) + eps_apo = kwargs.pop('eps_apo', None) + if thtcap is None: thtcap = np.max(loc[:, 0]) + if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) + return capsht.experimental.synthesis_general_cap(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) + +def adjoint_synthesis_general_cap(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): + """Wrapper to capsht adjoint_synthesis_general_cap function, hiding the choice of eps_apo + + If thtcap is not specified, it is set to np.max(loc[:, 0]) + + See ducc0.sht.adjoint_synthesis_general for arguments, optional arguments and outputs + + """ + thtcap = kwargs.pop('thtcap', None) + eps_apo = kwargs.pop('eps_apo', None) + if thtcap is None: thtcap = np.max(loc[:, 0]) + if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) + return capsht.experimental.adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) \ No newline at end of file From ab750de3fad8e8d6aaee48dd1bde561d59a46fe5 Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 1 Apr 2025 11:19:32 +0200 Subject: [PATCH 12/18] typo --- lenspyx/experimental.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lenspyx/experimental.py b/lenspyx/experimental.py index 394bafe..93f702c 100644 --- a/lenspyx/experimental.py +++ b/lenspyx/experimental.py @@ -9,8 +9,8 @@ print("synthesis_general_cap not found in capsht.experimental, are you up to date?") def _epsapo(thtcap, epsilon, lmax, dl_7=20): - dl = int(np.round(dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2)) - return np.sqrt(dl / lmax * np.pi / thtcap), + dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2 + return np.sqrt(dl / lmax * np.pi / thtcap) def synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): """Wrapper to capsht synthesis_general_cap function, hiding the choice of eps_apo From 564fd30e62758cd35c78c194bfb37a78417fd0bb Mon Sep 17 00:00:00 2001 From: jcarron Date: Thu, 3 Apr 2025 18:59:11 +0200 Subject: [PATCH 13/18] minors --- lenspyx/qest/qest.py | 26 +++++++++++++------------- lenspyx/qest/utils_qe.py | 8 ++++---- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/lenspyx/qest/qest.py b/lenspyx/qest/qest.py index 70ac64d..864e652 100644 --- a/lenspyx/qest/qest.py +++ b/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. @@ -36,9 +36,9 @@ def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, verbose=False, get_ 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 +105,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: diff --git a/lenspyx/qest/utils_qe.py b/lenspyx/qest/utils_qe.py index d51d722..73c6c73 100644 --- a/lenspyx/qest/utils_qe.py +++ b/lenspyx/qest/utils_qe.py @@ -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] From 2b4c18ac483d6480eb2b259382bc743335096006 Mon Sep 17 00:00:00 2001 From: jcarron Date: Thu, 24 Apr 2025 21:47:01 +0200 Subject: [PATCH 14/18] minor --- lenspyx/cachers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenspyx/cachers.py b/lenspyx/cachers.py index 3a1a354..d69620d 100644 --- a/lenspyx/cachers.py +++ b/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 From be2b1e6ad8f3a15e9ba6560bb7ff3d8991dabd8d Mon Sep 17 00:00:00 2001 From: jcarron Date: Wed, 7 May 2025 16:03:54 +0200 Subject: [PATCH 15/18] minors --- lenspyx/qest/qest.py | 1 + lenspyx/qest/utils_qe.py | 2 +- lenspyx/utils_hp.py | 7 +++++-- setup.py | 2 +- 4 files changed, 8 insertions(+), 4 deletions(-) diff --git a/lenspyx/qest/qest.py b/lenspyx/qest/qest.py index 864e652..760b754 100644 --- a/lenspyx/qest/qest.py +++ b/lenspyx/qest/qest.py @@ -32,6 +32,7 @@ def eval_qe(qe_key, lmax_ivf, cls_weight, get_alm, lmax_qlm, mmax=None, mmax_qlm 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)]) diff --git a/lenspyx/qest/utils_qe.py b/lenspyx/qest/utils_qe.py index 73c6c73..ed5c79b 100644 --- a/lenspyx/qest/utils_qe.py +++ b/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) diff --git a/lenspyx/utils_hp.py b/lenspyx/utils_hp.py index 02ccddd..06ed811 100644 --- a/lenspyx/utils_hp.py +++ b/lenspyx/utils_hp.py @@ -94,7 +94,7 @@ def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64): return alm -def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np.float64): +def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np.float64, rngen=None): """Creates Gaussian field alms from input cl dictionary Parameters @@ -109,6 +109,7 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np. Random generator seed rlm_dtype:(optional, defaults to np.float64) Precision of real components of the array (e.g. np.float32 for single precision output array) + rngen: (optional, defaults to None) Allows to use a custom random generator, seed is ignored in this case. Returns ------- @@ -116,6 +117,8 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np. harmonic coefficients of Gaussian field with lmax, mmax parameters """ + if rngen is not None: + assert hasattr(rngen, 'standard_normal'), 'rngen must have standard_normal method' lmax_cls = np.max([len(cl) - 1 for cl in cls.values()]) if lmax is None: lmax = lmax_cls @@ -158,7 +161,7 @@ def synalms(cls: dict, lmax:int, mmax:int or None, seed=None, rlm_dtype:type=np. m[:] = np.dot(v, np.dot(np.diag(np.sqrt(t)), v.T)) # Build phases: alm_size = Alm.getsize(lmax, mmax) - rng = default_rng(seed) + rng = default_rng(seed) if rngen is None else rngen phases = 1j * rng.standard_normal((ncomp, alm_size), dtype=rlm_dtype) phases += rng.standard_normal((ncomp, alm_size), dtype=rlm_dtype) phases *= np.sqrt(0.5) diff --git a/setup.py b/setup.py index ef75415..62c17c8 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,7 @@ setup( name='lenspyx', version=__version__, - packages=['lenspyx', 'lenspyx.remapping', 'lenspyx.tests', 'lenspyx.wigners'], + packages=['lenspyx', 'lenspyx.remapping', 'lenspyx.tests', 'lenspyx.wigners', 'lenspyx.qest'], url='https://github.com/carronj/lenspyx', author='Julien Carron', data_files=[('lenspyx/data/cls', ['lenspyx/data/cls/FFP10_wdipole_lensedCls.dat', From 67235fb3b00f60a6074cdd629f6a416d2872e2e4 Mon Sep 17 00:00:00 2001 From: jcarron Date: Tue, 26 Aug 2025 22:21:56 +0200 Subject: [PATCH 16/18] mmax --- lenspyx/remapping/utils_geom.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/lenspyx/remapping/utils_geom.py b/lenspyx/remapping/utils_geom.py index adb9fb9..9544cb0 100644 --- a/lenspyx/remapping/utils_geom.py +++ b/lenspyx/remapping/utils_geom.py @@ -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 @@ -343,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 From 75937ee74b57d2c68e07e0aba53681c149142737 Mon Sep 17 00:00:00 2001 From: jcarron Date: Sun, 8 Feb 2026 11:26:18 +0100 Subject: [PATCH 17/18] minors --- lenspyx/experimental.py | 122 +++++++++++++++++++++++++++++++++++++--- lenspyx/utils_hp.py | 26 ++++++--- 2 files changed, 133 insertions(+), 15 deletions(-) diff --git a/lenspyx/experimental.py b/lenspyx/experimental.py index 93f702c..537341e 100644 --- a/lenspyx/experimental.py +++ b/lenspyx/experimental.py @@ -1,18 +1,83 @@ 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 + from capsht.experimental import synthesis_general_cap, synthesis_general_band, adjoint_synthesis_general_cap, adjoint_synthesis_general_band except ImportError: - print("synthesis_general_cap not found in capsht.experimental, are you up to date?") + print("synthesis_general_cap or synthesis_general_band not found in capsht.experimental, are you up to date?") -def _epsapo(thtcap, epsilon, lmax, dl_7=20): - dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2 +def _epsapo(thtcap, epsilon, lmax, dl_7=15): + #dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2 + dl = dl_7 * ((-np.log10(epsilon) / (7 )) ** 1) ** 0.5 return np.sqrt(dl / lmax * np.pi / thtcap) -def synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): +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)) + 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) + 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, 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) + +def _synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): """Wrapper to capsht synthesis_general_cap function, hiding the choice of eps_apo If thtcap is not specified, it is set to np.max(loc[:, 0]) @@ -26,7 +91,7 @@ def synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) return capsht.experimental.synthesis_general_cap(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) -def adjoint_synthesis_general_cap(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): +def _adjoint_synthesis_general_cap(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): """Wrapper to capsht adjoint_synthesis_general_cap function, hiding the choice of eps_apo If thtcap is not specified, it is set to np.max(loc[:, 0]) @@ -38,4 +103,47 @@ def adjoint_synthesis_general_cap(map: np.ndarray, spin: int, lmax: int, loc: np eps_apo = kwargs.pop('eps_apo', None) if thtcap is None: thtcap = np.max(loc[:, 0]) if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) - return capsht.experimental.adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) \ No newline at end of file + return capsht.experimental.adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) + +def _synthesis_general_band(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): + """Wrapper to capsht synthesis_general_band function, hiding the choice of eps_apo + + If thtcap is not specified, it is set to np.max(loc[:, 0]) + + See ducc0.sht.synthesis_general for arguments, optional arguments and outputs + + """ + tht_min = kwargs.pop('tht_min', None) + tht_max = kwargs.pop('tht_max', None) + eps_apo = kwargs.pop('eps_apo', None) + if tht_max is None: tht_max = np.max(loc[:, 0]) + if tht_min is None: tht_min = np.min(loc[:, 0]) + if eps_apo is None: eps_apo = 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): + return capsht.experimental.synthesis_general_band(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, + thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs) + return syngducc(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, **kwargs) + +def _adjoint_synthesis_general_band(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): + """Wrapper to capsht adjoint_synthesis_general_band function, hiding the choice of eps_apo + + If thtcap is not specified, it is set to np.max(loc[:, 0]) + + See ducc0.sht.adjoint_synthesis_general for arguments, optional arguments and outputs + + """ + tht_min = kwargs.pop('tht_min', None) + tht_max = kwargs.pop('tht_max', None) + eps_apo = kwargs.pop('eps_apo', None) + if tht_max is None: tht_max = np.max(loc[:, 0]) + if tht_min is None: tht_min = np.min(loc[:, 0]) + if eps_apo is None: eps_apo = 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): + return capsht.experimental.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) + else: + return adjsyngducc(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, **kwargs) \ No newline at end of file diff --git a/lenspyx/utils_hp.py b/lenspyx/utils_hp.py index 06ed811..2a77fc0 100644 --- a/lenspyx/utils_hp.py +++ b/lenspyx/utils_hp.py @@ -9,7 +9,7 @@ def almxfl(alm:np.ndarray, fl:np.ndarray, mmax:int or None, inplace:bool): Parameters ---------- alm : array - The alm to multiply + The alm, or alms to multiply fl : array The function (at l=0..fl.size-1) by which alm must be multiplied. mmax : None or int @@ -24,20 +24,20 @@ def almxfl(alm:np.ndarray, fl:np.ndarray, mmax:int or None, inplace:bool): if inplace is True. """ - lmax = Alm.getlmax(alm.size, mmax) + lmax = Alm.getlmax(alm.shape[-1], mmax) if mmax is None or mmax < 0: mmax = lmax assert fl.size > lmax, (fl.size, lmax) if inplace: for m in range(mmax + 1): b = m * (2 * lmax + 1 - m) // 2 + m - alm[b:b + lmax - m + 1] *= fl[m:lmax+1] + alm[...,b:b + lmax - m + 1] *= fl[m:lmax+1] return else: ret = np.empty_like(alm) for m in range(mmax + 1): b = m * (2 * lmax + 1 - m) // 2 + m - ret[b:b + lmax - m + 1] = alm[b:b + lmax - m + 1] * fl[m:lmax+1] + ret[...,b:b + lmax - m + 1] = alm[...,b:b + lmax - m + 1] * fl[m:lmax+1] return ret @@ -63,7 +63,7 @@ def gauss_beam(fwhm:float, lmax:int): return bl -def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64): +def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64, seed=None, rngen=None, alms_r:np.ndarray=None): """Creates a Gaussian field alm from input cl array Parameters @@ -83,14 +83,24 @@ def synalm(cl:np.ndarray, lmax:int, mmax:int or None, rlm_dtype=np.float64): harmonic coefficients of Gaussian field with lmax, mmax parameters """ + if rngen is not None: + assert hasattr(rngen, 'standard_normal'), 'rngen must have standard_normal method' assert lmax + 1 <= cl.size if mmax is None or mmax < 0: mmax = lmax alm_size = Alm.getsize(lmax, mmax) - alm = rng.standard_normal(alm_size, dtype=rlm_dtype) + 1j * rng.standard_normal(alm_size, dtype=rlm_dtype) - almxfl(alm, np.sqrt(cl[:lmax+1] * 0.5), mmax, True) + rng = default_rng(seed) if rngen is None else rngen real_idcs = Alm.getidx(lmax, np.arange(lmax + 1, dtype=int), 0) - alm[real_idcs] = alm[real_idcs].real * np.sqrt(2.) + if alms_r is None: + alm = rng.standard_normal(alm_size, dtype=rlm_dtype) + 1j * rng.standard_normal(alm_size, dtype=rlm_dtype) + almxfl(alm, np.sqrt(cl[:lmax+1] * 0.5), mmax, True) + else: + assert alms_r.shape == (2, alm_size), alms_r.shape + rng.standard_normal((2, alm_size), dtype=rlm_dtype, out=alms_r) + almxfl(alms_r[0], np.sqrt(cl[:lmax+1] * 0.5), mmax, True) + almxfl(alms_r[1], np.sqrt(cl[:lmax+1] * 0.5), mmax, True) + alm = alms_r + alm[...,real_idcs] = alm[...,real_idcs].real * np.sqrt(2.) return alm From 0dae1ce65c68181a222bc6882e25542fb636e0f2 Mon Sep 17 00:00:00 2001 From: jcarron Date: Mon, 23 Feb 2026 19:46:53 +0100 Subject: [PATCH 18/18] new epsapo scheme --- lenspyx/experimental.py | 91 ++++++++--------------------------------- 1 file changed, 16 insertions(+), 75 deletions(-) diff --git a/lenspyx/experimental.py b/lenspyx/experimental.py index 537341e..f931849 100644 --- a/lenspyx/experimental.py +++ b/lenspyx/experimental.py @@ -9,9 +9,19 @@ 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, dl_7=15): +def _epsapo(thtcap, epsilon, lmax, version=1, dl_7=None): #dl = dl_7 * ((- np.log10(epsilon) + 1) / (7 + 1)) ** 2 - dl = dl_7 * ((-np.log10(epsilon) / (7 )) ** 1) ** 0.5 + 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, @@ -31,6 +41,7 @@ def synthesis_general(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, ep 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 @@ -38,9 +49,10 @@ def synthesis_general(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, ep 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, thtcap=thtcap, eps_apo=eps_apo, **kwargs) + 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) @@ -75,75 +87,4 @@ def adjoint_synthesis_general(map: np.ndarray, spin: int, lmax: int, loc: np.nda 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) - -def _synthesis_general_cap(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): - """Wrapper to capsht synthesis_general_cap function, hiding the choice of eps_apo - - If thtcap is not specified, it is set to np.max(loc[:, 0]) - - See ducc0.sht.synthesis_general for arguments, optional arguments and outputs - - """ - thtcap = kwargs.pop('thtcap', None) - eps_apo = kwargs.pop('eps_apo', None) - if thtcap is None: thtcap = np.max(loc[:, 0]) - if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) - return capsht.experimental.synthesis_general_cap(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) - -def _adjoint_synthesis_general_cap(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): - """Wrapper to capsht adjoint_synthesis_general_cap function, hiding the choice of eps_apo - - If thtcap is not specified, it is set to np.max(loc[:, 0]) - - See ducc0.sht.adjoint_synthesis_general for arguments, optional arguments and outputs - - """ - thtcap = kwargs.pop('thtcap', None) - eps_apo = kwargs.pop('eps_apo', None) - if thtcap is None: thtcap = np.max(loc[:, 0]) - if eps_apo is None: eps_apo = _epsapo(thtcap, epsilon, lmax) - return capsht.experimental.adjoint_synthesis_general_cap(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, thtcap=thtcap, eps_apo=eps_apo, **kwargs) - -def _synthesis_general_band(alm: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): - """Wrapper to capsht synthesis_general_band function, hiding the choice of eps_apo - - If thtcap is not specified, it is set to np.max(loc[:, 0]) - - See ducc0.sht.synthesis_general for arguments, optional arguments and outputs - - """ - tht_min = kwargs.pop('tht_min', None) - tht_max = kwargs.pop('tht_max', None) - eps_apo = kwargs.pop('eps_apo', None) - if tht_max is None: tht_max = np.max(loc[:, 0]) - if tht_min is None: tht_min = np.min(loc[:, 0]) - if eps_apo is None: eps_apo = 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): - return capsht.experimental.synthesis_general_band(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, - thta=tht_min, thtb=tht_max, eps_apo=eps_apo, **kwargs) - return syngducc(alm=alm, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, **kwargs) - -def _adjoint_synthesis_general_band(map: np.ndarray, spin: int, lmax: int, loc: np.ndarray, epsilon: float, **kwargs): - """Wrapper to capsht adjoint_synthesis_general_band function, hiding the choice of eps_apo - - If thtcap is not specified, it is set to np.max(loc[:, 0]) - - See ducc0.sht.adjoint_synthesis_general for arguments, optional arguments and outputs - - """ - tht_min = kwargs.pop('tht_min', None) - tht_max = kwargs.pop('tht_max', None) - eps_apo = kwargs.pop('eps_apo', None) - if tht_max is None: tht_max = np.max(loc[:, 0]) - if tht_min is None: tht_min = np.min(loc[:, 0]) - if eps_apo is None: eps_apo = 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): - return capsht.experimental.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) - else: - return adjsyngducc(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, **kwargs) \ No newline at end of file + return syngducc(map=map, spin=spin, lmax=lmax, loc=loc, epsilon=epsilon, **kwargs) \ No newline at end of file