gh-501: change imports to be directly from glass

docs/user/how-glass-works.rst

@@ -30,13 +30,13 @@ which are flat and non-overlapping.
 
 .. plot::
 
-    from glass.shells import redshift_grid, tophat_windows
+    import glass
 
     # create a redshift grid for shell edges
-    zs = redshift_grid(0., 0.5, dz=0.1)
+    zs = glass.redshift_grid(0., 0.5, dz=0.1)
 
     # create the top hat windows
-    ws = tophat_windows(zs)
+    ws = glass.tophat_windows(zs)
 
     # plot each window
     for i, (za, wa, zeff) in enumerate(ws):
@@ -74,17 +74,16 @@ included:
 
 .. plot::
 
-    from glass.shells import (redshift_grid, tophat_windows, linear_windows,
-                              cubic_windows)
+    import glass
 
-    plot_windows = [tophat_windows, linear_windows,
-                    cubic_windows]
+    plot_windows = [glass.tophat_windows, glass.linear_windows,
+                    glass.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)
+    zs = glass.redshift_grid(0., 0.5, dz=0.1)
     zt = np.linspace(0., 0.5, 200)
 
     for ax in axes.flat:

glass/fields.py

@@ -13,17 +13,16 @@
 from gaussiancl import gaussiancl
 from transformcl import cltovar
 
-from glass import grf
+import glass
+import glass.grf
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Generator, Iterable, Iterator, Sequence
     from typing import Any
 
     from numpy.typing import NDArray
 
-    from glass.shells import RadialWindow
-
-    Fields = Sequence[grf.Transformation]
+    Fields = Sequence[glass.grf.Transformation]
     Cls = Sequence[NDArray[Any]]
 
 
@@ -640,7 +639,7 @@ def effective_cls(
     return out
 
 
-def gaussian_fields(shells: Sequence[RadialWindow]) -> Sequence[grf.Normal]:
+def gaussian_fields(shells: Sequence[glass.RadialWindow]) -> Sequence[glass.grf.Normal]:
     """
     Create Gaussian random fields for radial windows *shells*.
 
@@ -654,13 +653,13 @@ def gaussian_fields(shells: Sequence[RadialWindow]) -> Sequence[grf.Normal]:
         A sequence describing the Gaussian random fields.
 
     """
-    return [grf.Normal() for _shell in shells]
+    return [glass.grf.Normal() for _shell in shells]
 
 
 def lognormal_fields(
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     shift: Callable[[float], float] | None = None,
-) -> Sequence[grf.Lognormal]:
+) -> Sequence[glass.grf.Lognormal]:
     """
     Create lognormal fields for radial windows *shells*.  If *shifts* is
     given, it must be a callable that returns a lognormal shift (i.e.
@@ -681,7 +680,7 @@ def lognormal_fields(
     if shift is None:
         shift = lambda _z: 1.0  # noqa: E731
 
-    return [grf.Lognormal(shift(shell.zeff)) for shell in shells]
+    return [glass.grf.Lognormal(shift(shell.zeff)) for shell in shells]
 
 
 def compute_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
@@ -711,7 +710,7 @@ def compute_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
 
     gls = []
     for i, j, cl in enumerate_spectra(spectra):
-        gl = grf.compute(cl, fields[i], fields[j])
+        gl = glass.grf.compute(cl, fields[i], fields[j])
         gls.append(gl)
     return gls
 
@@ -757,7 +756,7 @@ def solve_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
             monopole = 0.0 if cl[0] == 0 else None
 
             # call solver
-            gl, _cl_out, info = grf.solve(cl, t1, t2, pad=pad, monopole=monopole)
+            gl, _cl_out, info = glass.grf.solve(cl, t1, t2, pad=pad, monopole=monopole)
 
             # warn if solver didn't converge
             if info == 0:

glass/galaxies.py

@@ -25,19 +25,18 @@
 import healpix
 import numpy as np
 
-from glass.core.array import broadcast_leading_axes, cumulative_trapezoid
+import glass
+import glass.core.array
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
 
     from cosmology import Cosmology
 
-    from glass.shells import RadialWindow
-
 
 def redshifts(
     n: int | NDArray[np.float64],
-    w: RadialWindow,
+    w: glass.RadialWindow,
     *,
     rng: np.random.Generator | None = None,
 ) -> NDArray[np.float64]:
@@ -118,7 +117,7 @@ def redshifts_from_nz(
         rng = np.random.default_rng()
 
     # bring inputs' leading axes into common shape
-    dims, *rest = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
+    dims, *rest = glass.core.array.broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
     count_out, z_out, nz_out = rest
 
     # list of results for all dimensions
@@ -130,7 +129,7 @@ def redshifts_from_nz(
     # go through extra dimensions; also works if dims is empty
     for k in np.ndindex(dims):
         # compute the CDF of each galaxy population
-        cdf = cumulative_trapezoid(nz_out[k], z_out[k], dtype=float)
+        cdf = glass.core.array.cumulative_trapezoid(nz_out[k], z_out[k], dtype=float)
         cdf /= cdf[-1]
 
         # sample redshifts and store result

glass/grf/_solver.py

@@ -5,13 +5,13 @@
 import numpy as np
 from transformcl import cltocorr, corrtocl
 
+import glass.grf
+
 if TYPE_CHECKING:
     from typing import Any
 
     from numpy.typing import NDArray
 
-from glass.grf._core import Transformation, corr, dcorr, icorr
-
 
 def _relerr(dx: NDArray[Any], x: NDArray[Any]) -> float:
     """Compute the relative error max(|dx/x|)."""
@@ -21,8 +21,8 @@ def _relerr(dx: NDArray[Any], x: NDArray[Any]) -> float:
 
 def solve(  # noqa: PLR0912, PLR0913
     cl: NDArray[Any],
-    t1: Transformation,
-    t2: Transformation | None = None,
+    t1: glass.grf.Transformation,
+    t2: glass.grf.Transformation | None = None,
     *,
     pad: int = 0,
     initial: NDArray[Any] | None = None,
@@ -91,7 +91,7 @@ def solve(  # noqa: PLR0912, PLR0913
         raise ValueError(msg)
 
     if initial is None:
-        gl = corrtocl(icorr(t1, t2, cltocorr(cl)))
+        gl = corrtocl(glass.grf.icorr(t1, t2, cltocorr(cl)))
     else:
         gl = np.zeros(n)
         gl[: len(initial)] = initial[:n]
@@ -100,7 +100,7 @@ def solve(  # noqa: PLR0912, PLR0913
         gl[0] = monopole
 
     gt = cltocorr(np.pad(gl, (0, pad)))
-    rl = corrtocl(corr(t1, t2, gt))
+    rl = corrtocl(glass.grf.corr(t1, t2, gt))
     fl = rl[:n] - cl
     if monopole is not None:
         fl[0] = 0
@@ -116,7 +116,7 @@ def solve(  # noqa: PLR0912, PLR0913
             break
 
         ft = cltocorr(np.pad(fl, (0, pad)))
-        dt = dcorr(t1, t2, gt)
+        dt = glass.grf.dcorr(t1, t2, gt)
         xl = -corrtocl(ft / dt)[:n]
         if monopole is not None:
             xl[0] = 0
@@ -126,7 +126,7 @@ def solve(  # noqa: PLR0912, PLR0913
         while True:
             gl_ = gl + xl
             gt_ = cltocorr(np.pad(gl_, (0, pad)))
-            rl_ = corrtocl(corr(t1, t2, gt_))
+            rl_ = corrtocl(glass.grf.corr(t1, t2, gt_))
             fl_ = rl_[:n] - cl
             if monopole is not None:
                 fl_[0] = 0

glass/grf/_transformations.py

@@ -8,7 +8,7 @@
 
 from numpy.typing import NDArray  # noqa: TC002
 
-from glass.grf._core import corr, dcorr, icorr
+from glass.grf import corr, dcorr, icorr
 
 
 @dataclass

glass/lensing.py

@@ -43,7 +43,7 @@
 
     from cosmology import Cosmology
 
-    from glass.shells import RadialWindow
+    import glass
 
 
 @overload
@@ -351,7 +351,7 @@ def shear_from_convergence(
     transform.
 
     .. deprecated:: 2023.6
-       Use the more general :func:`from_convergence` function instead.
+       Use the more general :func:`glass.from_convergence` function instead.
 
     Parameters
     ----------
@@ -420,7 +420,7 @@ def __init__(self, cosmo: Cosmology) -> None:
         self.kappa2: NDArray[np.float64] | None = None
         self.kappa3: NDArray[np.float64] | None = None
 
-    def add_window(self, delta: NDArray[np.float64], w: RadialWindow) -> None:
+    def add_window(self, delta: NDArray[np.float64], w: glass.RadialWindow) -> None:
         """
         Add a mass plane from a window function to the convergence.
 
@@ -527,7 +527,7 @@ def wlens(self) -> float:
 
 
 def multi_plane_matrix(
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     cosmo: Cosmology,
 ) -> NDArray[np.float64]:
     """
@@ -555,7 +555,7 @@ def multi_plane_matrix(
 
 def multi_plane_weights(
     weights: NDArray[np.float64],
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     cosmo: Cosmology,
 ) -> NDArray[np.float64]:
     """

glass/observations.py

@@ -34,7 +34,7 @@
 import healpy as hp
 import numpy as np
 
-from glass.core.array import cumulative_trapezoid
+import glass.core.array
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
@@ -266,7 +266,7 @@ def equal_dens_zbins(
     # first compute the cumulative integral (by trapezoidal rule)
     # then normalise: the first z is at CDF = 0, the last z at CDF = 1
     # interpolate to find the z values at CDF = i/nbins for i = 0, ..., nbins
-    cuml_nz = cumulative_trapezoid(nz, z)
+    cuml_nz = glass.core.array.cumulative_trapezoid(nz, z)
     cuml_nz /= cuml_nz[[-1]]
     zbinedges = np.interp(np.linspace(0, 1, nbins + 1), cuml_nz, z)
 

glass/points.py

@@ -36,23 +36,22 @@
 import healpix
 import numpy as np
 
-from glass.core.array import broadcast_first, broadcast_leading_axes, trapezoid_product
+import glass
+import glass.core.array
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Generator
 
     from numpy.typing import NDArray
 
-    from glass.shells import RadialWindow
-
 
 ARCMIN2_SPHERE = 60**6 // 100 / np.pi
 
 
 def effective_bias(
     z: NDArray[np.float64],
     bz: NDArray[np.float64],
-    w: RadialWindow,
+    w: glass.RadialWindow,
 ) -> float | NDArray[np.double]:
     r"""
     Effective bias parameter from a redshift-dependent bias function.
@@ -86,7 +85,7 @@ def effective_bias(
 
     """
     norm = np.trapezoid(w.wa, w.za)
-    return trapezoid_product((z, bz), (w.za, w.wa)) / norm
+    return glass.core.array.trapezoid_product((z, bz), (w.za, w.wa)) / norm
 
 
 def linear_bias(
@@ -189,7 +188,8 @@ def positions_from_delta(  # noqa: PLR0912, PLR0913, PLR0915
     bias_model
         The bias model to apply. If a string, refers to a function in
         the :mod:`~glass.points` module, e.g. ``'linear'`` for
-        :func:`linear_bias()` or ``'loglinear'`` for :func:`loglinear_bias`.
+        :func:`glass.linear_bias()` or ``'glass.loglinear'`` for
+        :func:`glass.loglinear_bias`.
     remove_monopole
         If true, the monopole of the density contrast
         after biasing is fixed to zero.
@@ -232,7 +232,7 @@ def positions_from_delta(  # noqa: PLR0912, PLR0913, PLR0915
         inputs.append((bias, 0))
     if vis is not None:
         inputs.append((vis, 1))
-    dims, *rest = broadcast_leading_axes(*inputs)
+    dims, *rest = glass.core.array.broadcast_leading_axes(*inputs)
     ngal, delta, *rest = rest
     if bias is not None:
         bias, *rest = rest
@@ -409,7 +409,7 @@ def position_weights(
     """
     # bring densities and bias into the same shape
     if bias is not None:
-        densities, bias = broadcast_first(densities, bias)
+        densities, bias = glass.core.array.broadcast_first(densities, bias)
     # normalise densities after shape has been fixed
     densities = densities / np.sum(densities, axis=0)
     # apply bias after normalisation