From 8dfe40f1e5657ec4ffb955d1ad8b7e76a322c85a Mon Sep 17 00:00:00 2001 From: Nicolas Tessore Date: Sun, 30 Jul 2023 16:18:52 +0100 Subject: [PATCH] ENH(shells): radial windows for linear and cubic interpolation (#119) Adds window functions `linear_windows()` and `cubic_windows()` that correspond to linear and cubic spline interpolation, respectively. Documentation is updated to give an overview of available window functions. Closes: #118 Added: New `linear_windows()` and `cubic_windows()` window functions for shells. --- docs/user/how-glass-works.rst | 45 ++++++++++++- glass/shells.py | 120 ++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+), 1 deletion(-) diff --git a/docs/user/how-glass-works.rst b/docs/user/how-glass-works.rst index 8fbfc003..bd94ba95 100644 --- a/docs/user/how-glass-works.rst +++ b/docs/user/how-glass-works.rst @@ -44,7 +44,6 @@ which are flat and non-overlapping. plt.fill_between(za, np.zeros_like(wa), wa, alpha=0.5) plt.annotate(f'shell {i+1}', (zeff, 0.5), ha='center', va='center') - plt.ylim(0., 2.) plt.xlabel('redshift $z$') plt.ylabel('window function $W(z)$') plt.tight_layout() @@ -63,6 +62,48 @@ then simulates the (radially) continuous field :math:`F(z)` as the (radially) discretised fields :math:`F_i`. +.. _user-window-functions: + +Window functions +^^^^^^^^^^^^^^^^ + +*GLASS* supports arbitrary window functions (although the computation of +:ref:`line-of-sight integrals ` makes some assumptions). +The following :ref:`window functions ` are +included: + +.. plot:: + + from glass.shells import (redshift_grid, tophat_windows, linear_windows, + cubic_windows) + + plot_windows = [tophat_windows, linear_windows, + cubic_windows] + nr = (len(plot_windows)+1)//2 + + fig, axes = plt.subplots(nr, 2, figsize=(8, nr*3), layout="constrained", + squeeze=False, sharex=False, sharey=True) + + zs = redshift_grid(0., 0.5, dz=0.1) + zt = np.linspace(0., 0.5, 200) + + for ax in axes.flat: + ax.axis(False) + for windows, ax in zip(plot_windows, axes.flat): + ws = windows(zs) + wt = np.zeros_like(zt) + ax.axis(True) + ax.set_title(windows.__name__) + for i, (za, wa, zeff) in enumerate(ws): + wt += np.interp(zt, za, wa, left=0., right=0.) + ax.fill_between(za, np.zeros_like(wa), wa, alpha=0.5) + ax.plot(zt, wt, c="k", lw=2) + for ax in axes.flat: + ax.set_xlabel("redshift $z$") + for ax in axes[:, 0]: + ax.set_ylabel("window function $W(z)$") + + Angular discretisation ---------------------- @@ -89,6 +130,8 @@ difference between the two is whether or not the pixel window function is applied to the spherical harmonic expansion of the fields. +.. _user-los-integrals: + Line-of-sight integrals ----------------------- diff --git a/glass/shells.py b/glass/shells.py index f16e9667..02105f59 100644 --- a/glass/shells.py +++ b/glass/shells.py @@ -10,11 +10,15 @@ matter shells, i.e. the radial discretisation of the light cone. +.. _reference-window-functions: + Window functions ---------------- .. autoclass:: RadialWindow .. autofunction:: tophat_windows +.. autofunction:: linear_windows +.. autofunction:: cubic_windows Window function tools @@ -151,6 +155,10 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3, ws : (N,) list of :class:`RadialWindow` List of window functions. + See Also + -------- + :ref:`user-window-functions` + ''' if len(zbins) < 2: raise ValueError('zbins must have at least two entries') @@ -173,6 +181,118 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3, return ws +def linear_windows(zgrid: ArrayLike1D, dz: float = 1e-3, + weight: Optional[WeightFunc] = None, + ) -> List[RadialWindow]: + '''Linear interpolation window functions. + + Uses the *N+2* given redshifts as nodes to construct *N* triangular + window functions between the first and last node. These correspond + to linear interpolation of radial functions. The redshift spacing + of the windows is approximately equal to ``dz``. + + An optional weight function :math:`w(z)` can be given using + ``weight``; it is applied to the triangular windows. + + The resulting windows functions are :class:`RadialWindow` instances. + Their effective redshifts correspond to the given nodes. + + Parameters + ---------- + zgrid : (N+2,) array_like + Redshift grid for the triangular window functions. + dz : float, optional + Approximate spacing of the redshift grid. + weight : callable, optional + If given, a weight function to be applied to the window + functions. + + Returns + ------- + ws : (N,) list of :class:`RadialWindow` + List of window functions. + + See Also + -------- + :ref:`user-window-functions` + + ''' + if len(zgrid) < 3: + raise ValueError('nodes must have at least 3 entries') + if zgrid[0] != 0: + warnings.warn('first triangular window does not start at z=0') + + ws = [] + for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]): + n = max(round((zmid - zmin)/dz), 2) - 1 + m = max(round((zmax - zmid)/dz), 2) + z = np.concatenate([np.linspace(zmin, zmid, n, endpoint=False), + np.linspace(zmid, zmax, m)]) + w = np.concatenate([np.linspace(0., 1., n, endpoint=False), + np.linspace(1., 0., m)]) + if weight is not None: + w *= weight(z) + ws.append(RadialWindow(z, w, zmid)) + return ws + + +def cubic_windows(zgrid: ArrayLike1D, dz: float = 1e-3, + weight: Optional[WeightFunc] = None, + ) -> List[RadialWindow]: + '''Cubic interpolation window functions. + + Uses the *N+2* given redshifts as nodes to construct *N* cubic + Hermite spline window functions between the first and last node. + These correspond to cubic spline interpolation of radial functions. + The redshift spacing of the windows is approximately equal to + ``dz``. + + An optional weight function :math:`w(z)` can be given using + ``weight``; it is applied to the cubic spline windows. + + The resulting windows functions are :class:`RadialWindow` instances. + Their effective redshifts correspond to the given nodes. + + Parameters + ---------- + zgrid : (N+2,) array_like + Redshift grid for the cubic spline window functions. + dz : float, optional + Approximate spacing of the redshift grid. + weight : callable, optional + If given, a weight function to be applied to the window + functions. + + Returns + ------- + ws : (N,) list of :class:`RadialWindow` + List of window functions. + + See Also + -------- + :ref:`user-window-functions` + + ''' + if len(zgrid) < 3: + raise ValueError('nodes must have at least 3 entries') + if zgrid[0] != 0: + warnings.warn('first cubic spline window does not start at z=0') + + ws = [] + for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]): + n = max(round((zmid - zmin)/dz), 2) - 1 + m = max(round((zmax - zmid)/dz), 2) + z = np.concatenate([np.linspace(zmin, zmid, n, endpoint=False), + np.linspace(zmid, zmax, m)]) + u = np.linspace(0., 1., n, endpoint=False) + v = np.linspace(1., 0., m) + w = np.concatenate([u**2*(3-2*u), v**2*(3-2*v)]) + if weight is not None: + w *= weight(z) + ws.append(RadialWindow(z, w, zmid)) + return ws + + def restrict(z: ArrayLike1D, f: ArrayLike1D, w: RadialWindow ) -> Tuple[np.ndarray, np.ndarray]: '''Restrict a function to a redshift window.