diff --git a/.github/workflows/pull_request.yml b/.github/workflows/pull_request.yml
index 5fba75cc..222e10c0 100644
--- a/.github/workflows/pull_request.yml
+++ b/.github/workflows/pull_request.yml
@@ -10,6 +10,7 @@ on:
 
 jobs:
   check-pr-title:
+    if: ${{ github.actor != 'pre-commit-ci[bot]' }}
     runs-on: ubuntu-latest
     steps:
       - uses: thehanimo/pr-title-checker@v1.4.2
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index f4821173..fe6bd6a3 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,3 +1,7 @@
+ci:
+  autoupdate_commit_msg: "MNT: update pre-commit hooks"
+  autofix_commit_msg: "STY: pre-commit fixes"
+
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.6.0
diff --git a/docs/user/how-glass-works.rst b/docs/user/how-glass-works.rst
index bd94ba95..d873236e 100644
--- a/docs/user/how-glass-works.rst
+++ b/docs/user/how-glass-works.rst
@@ -20,11 +20,11 @@ Radial discretisation
 The discretisation in the radial (line of sight) direction is done in *GLASS*
 using the concept of a :term:`radial window`, which consists of a window
 function :math:`W` that assigns a weight :math:`W(z)` to each redshift
-:math:`z`.  In the *GLASS* code, the :class:`~glass.shells.RadialWindow` named
+:math:`z`.  In the *GLASS* code, the :class:`~glass.RadialWindow` named
 tuple is used to define radial windows.
 
 A sequence :math:`W_1, W_2, \ldots` of such window functions defines the shells
-of the simulation.  For example, the :func:`~glass.shells.tophat_windows`
+of the simulation.  For example, the :func:`~glass.tophat_windows`
 function takes redshift boundaries and returns a sequence of top hat windows,
 which are flat and non-overlapping.
 
@@ -169,7 +169,7 @@ shells:
 
 In short, the requirements say that each shell has an effective redshift which
 partitions the window functions of all other shells. In *GLASS*, it is stored
-as the ``zeff`` attribute of :class:`~glass.shells.RadialWindow`.  Functions
+as the ``zeff`` attribute of :class:`~glass.RadialWindow`.  Functions
 that construct a list of windows for shells should ensure these requirements
 are met.
 
diff --git a/examples/1-basic/density.ipynb b/examples/1-basic/density.ipynb
index 8dfa74b0..f1e778cd 100644
--- a/examples/1-basic/density.ipynb
+++ b/examples/1-basic/density.ipynb
@@ -42,11 +42,8 @@
     "import camb\n",
     "from cosmology import Cosmology\n",
     "\n",
-    "# GLASS imports: matter, random fields, random points, galaxies\n",
-    "import glass.fields\n",
-    "import glass.galaxies\n",
-    "import glass.points\n",
-    "import glass.shells\n",
+    "# import GLASS for matter, random fields, random points, galaxies\n",
+    "import glass\n",
     "\n",
     "# cosmology for the simulation\n",
     "h = 0.7\n",
@@ -68,10 +65,10 @@
     "cosmo = Cosmology.from_camb(pars)\n",
     "\n",
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
+    "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
     "\n",
     "# linear radial window functions\n",
-    "ws = glass.shells.linear_windows(zb)\n",
+    "ws = glass.linear_windows(zb)\n",
     "\n",
     "# load the angular matter power spectra previously computed with CAMB\n",
     "cls = np.load(\"cls.npy\")"
@@ -99,10 +96,10 @@
    "outputs": [],
    "source": [
     "# compute Gaussian cls for lognormal fields for 3 correlated shells\n",
-    "gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
+    "gls = glass.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
     "\n",
     "# generator for lognormal matter fields\n",
-    "matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)"
+    "matter = glass.generate_lognormal(gls, nside, ncorr=3)"
    ]
   },
   {
@@ -131,7 +128,7 @@
     "dndz = np.full_like(z, 0.01)\n",
     "\n",
     "# distribute the dN/dz over the linear window functions\n",
-    "ngal = glass.shells.partition(z, dndz, ws)"
+    "ngal = glass.partition(z, dndz, ws)"
    ]
   },
   {
@@ -166,12 +163,9 @@
     "# simulate and add galaxies in each matter shell to cube\n",
     "for i, delta_i in enumerate(matter):\n",
     "    # simulate positions from matter density\n",
-    "    for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(\n",
-    "        ngal[i],\n",
-    "        delta_i,\n",
-    "    ):\n",
+    "    for gal_lon, gal_lat, gal_count in glass.positions_from_delta(ngal[i], delta_i):\n",
     "        # sample redshifts uniformly in shell\n",
-    "        gal_z = glass.galaxies.redshifts(gal_count, ws[i])\n",
+    "        gal_z = glass.redshifts(gal_count, ws[i])\n",
     "\n",
     "        # add counts to cube\n",
     "        z1 = gal_z * np.cos(np.deg2rad(gal_lon)) * np.cos(np.deg2rad(gal_lat))\n",
diff --git a/examples/1-basic/lensing.ipynb b/examples/1-basic/lensing.ipynb
index 5c5ba2ef..e32ff2aa 100644
--- a/examples/1-basic/lensing.ipynb
+++ b/examples/1-basic/lensing.ipynb
@@ -46,11 +46,8 @@
     "import camb\n",
     "from cosmology import Cosmology\n",
     "\n",
-    "# GLASS imports\n",
-    "import glass.fields\n",
-    "import glass.galaxies\n",
-    "import glass.lensing\n",
-    "import glass.shells\n",
+    "# import GLASS\n",
+    "import glass\n",
     "\n",
     "# cosmology for the simulation\n",
     "h = 0.7\n",
@@ -72,10 +69,10 @@
     "cosmo = Cosmology.from_camb(pars)\n",
     "\n",
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
+    "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
     "\n",
     "# linear radial window functions\n",
-    "ws = glass.shells.linear_windows(zb)\n",
+    "ws = glass.linear_windows(zb)\n",
     "\n",
     "# load the angular matter power spectra previously computed with CAMB\n",
     "cls = np.load(\"cls.npy\")"
@@ -104,10 +101,10 @@
    "source": [
     "# compute Gaussian cls for lognormal fields for 3 correlated shells\n",
     "# putting nside here means that the HEALPix pixel window function is applied\n",
-    "gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
+    "gls = glass.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
     "\n",
     "# generator for lognormal matter fields\n",
-    "matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)"
+    "matter = glass.generate_lognormal(gls, nside, ncorr=3)"
    ]
   },
   {
@@ -132,7 +129,7 @@
    "outputs": [],
    "source": [
     "# this will compute the convergence field iteratively\n",
-    "convergence = glass.lensing.MultiPlaneConvergence(cosmo)"
+    "convergence = glass.MultiPlaneConvergence(cosmo)"
    ]
   },
   {
@@ -162,7 +159,7 @@
     "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n",
     "\n",
     "# distribute dN/dz over the radial window functions\n",
-    "ngal = glass.shells.partition(z, dndz, ws)"
+    "ngal = glass.partition(z, dndz, ws)"
    ]
   },
   {
@@ -202,7 +199,7 @@
     "    kappa_i = convergence.kappa\n",
     "\n",
     "    # compute shear field\n",
-    "    gamm1_i, gamm2_i = glass.lensing.shear_from_convergence(kappa_i)\n",
+    "    gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n",
     "\n",
     "    # add to mean fields using the galaxy number density as weight\n",
     "    kappa_bar += ngal[i] * kappa_i\n",
diff --git a/examples/1-basic/matter.ipynb b/examples/1-basic/matter.ipynb
index 24d4c724..45585bc2 100644
--- a/examples/1-basic/matter.ipynb
+++ b/examples/1-basic/matter.ipynb
@@ -42,9 +42,8 @@
     "import camb\n",
     "from cosmology import Cosmology\n",
     "\n",
-    "# these are the GLASS imports: matter and random fields\n",
-    "import glass.fields\n",
-    "import glass.shells\n",
+    "# import GLASS for matter and random fields\n",
+    "import glass\n",
     "\n",
     "# cosmology for the simulation\n",
     "h = 0.7\n",
@@ -67,16 +66,16 @@
     "cosmo = Cosmology.from_camb(pars)\n",
     "\n",
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
+    "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
     "\n",
     "# load precomputed angular matter power spectra\n",
     "cls = np.load(\"cls.npy\")\n",
     "\n",
     "# compute Gaussian cls for lognormal fields with 3 correlated shells\n",
-    "gls = glass.fields.lognormal_gls(cls, ncorr=3, nside=nside)\n",
+    "gls = glass.lognormal_gls(cls, ncorr=3, nside=nside)\n",
     "\n",
     "# this generator will yield the matter fields in each shell\n",
-    "matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)"
+    "matter = glass.generate_lognormal(gls, nside, ncorr=3)"
    ]
   },
   {
diff --git a/examples/1-basic/photoz.ipynb b/examples/1-basic/photoz.ipynb
index ff3b6dbb..9e500cbe 100644
--- a/examples/1-basic/photoz.ipynb
+++ b/examples/1-basic/photoz.ipynb
@@ -38,9 +38,8 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
-    "# GLASS imports: matter shells, galaxies, random points, and observational\n",
-    "import glass.galaxies\n",
-    "import glass.observations\n",
+    "# import GLASS for matter shells, galaxies, random points, and observational\n",
+    "import glass\n",
     "\n",
     "# how many arcmin2 over the entire sphere\n",
     "from glass.core.constants import ARCMIN2_SPHERE\n",
@@ -53,7 +52,7 @@
     "\n",
     "# parametric galaxy redshift distribution\n",
     "z = np.linspace(0, 3, 301)\n",
-    "dndz = n_arcmin2 * glass.observations.smail_nz(z, 1.0, 2.2, 1.5)\n",
+    "dndz = n_arcmin2 * glass.smail_nz(z, 1.0, 2.2, 1.5)\n",
     "\n",
     "# compute the over galaxy number density on the sphere\n",
     "ngal = np.trapz(dndz, z)"
@@ -86,10 +85,10 @@
     "n = np.random.poisson(ngal * ARCMIN2_SPHERE)\n",
     "\n",
     "# sample n true redshifts\n",
-    "ztrue = glass.galaxies.redshifts_from_nz(n, z, dndz)\n",
+    "ztrue = glass.redshifts_from_nz(n, z, dndz)\n",
     "\n",
     "# sample n photometric redshifts\n",
-    "zphot = glass.galaxies.gaussian_phz(ztrue, phz_sigma_0)"
+    "zphot = glass.gaussian_phz(ztrue, phz_sigma_0)"
    ]
   },
   {
@@ -146,7 +145,7 @@
    "metadata": {},
    "source": [
     "Now define a number of photometric redshift bins.  They are chosen by the\n",
-    ":func:`~glass.observations.equal_dens_zbins` function to produce the same\n",
+    ":func:`~glass.equal_dens_zbins` function to produce the same\n",
     "number of galaxies in each bin.\n",
     "\n"
    ]
@@ -166,7 +165,7 @@
    "outputs": [],
    "source": [
     "nbins = 5\n",
-    "zbins = glass.observations.equal_dens_zbins(z, dndz, nbins)"
+    "zbins = glass.equal_dens_zbins(z, dndz, nbins)"
    ]
   },
   {
@@ -175,7 +174,7 @@
    "source": [
     "After the photometric bins are defined, make histograms of the *true* redshift\n",
     "distribution $n(z)$ using the *photometric* redshifts for binning.  Use\n",
-    "the :func:`~glass.observations.tomo_nz_gausserr()` function to also plot the\n",
+    "the :func:`~glass.tomo_nz_gausserr()` function to also plot the\n",
     "expected tomographic redshift distributions with the same model.\n",
     "\n"
    ]
@@ -210,7 +209,7 @@
     }
    ],
    "source": [
-    "tomo_nz = glass.observations.tomo_nz_gausserr(z, dndz, phz_sigma_0, zbins)\n",
+    "tomo_nz = glass.tomo_nz_gausserr(z, dndz, phz_sigma_0, zbins)\n",
     "tomo_nz *= ARCMIN2_SPHERE * (z[-1] - z[0]) / 40\n",
     "\n",
     "for (z1, z2), nz in zip(zbins, tomo_nz):\n",
diff --git a/examples/1-basic/shells.ipynb b/examples/1-basic/shells.ipynb
index 4b901ac5..f284077b 100644
--- a/examples/1-basic/shells.ipynb
+++ b/examples/1-basic/shells.ipynb
@@ -44,8 +44,8 @@
     "import camb\n",
     "from cosmology import Cosmology\n",
     "\n",
+    "import glass\n",
     "import glass.ext.camb\n",
-    "import glass.shells\n",
     "\n",
     "# cosmology for the simulation\n",
     "h = 0.7\n",
@@ -67,10 +67,10 @@
     "cosmo = Cosmology.from_camb(pars)\n",
     "\n",
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zgrid = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
+    "zgrid = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
     "\n",
     "# triangular radial windows, equivalent to linear interpolation of n(z)\n",
-    "shells = glass.shells.linear_windows(zgrid)\n",
+    "shells = glass.linear_windows(zgrid)\n",
     "\n",
     "# compute angular matter power spectra with CAMB\n",
     "cls = glass.ext.camb.matter_cls(pars, lmax, shells)"
diff --git a/examples/2-advanced/cosmic_shear.ipynb b/examples/2-advanced/cosmic_shear.ipynb
index a4b28845..146ddeeb 100644
--- a/examples/2-advanced/cosmic_shear.ipynb
+++ b/examples/2-advanced/cosmic_shear.ipynb
@@ -42,12 +42,7 @@
     "from cosmology import Cosmology\n",
     "\n",
     "# GLASS modules: cosmology and everything in the glass namespace\n",
-    "import glass.fields\n",
-    "import glass.galaxies\n",
-    "import glass.lensing\n",
-    "import glass.points\n",
-    "import glass.shapes\n",
-    "import glass.shells\n",
+    "import glass\n",
     "from glass.core.constants import ARCMIN2_SPHERE\n",
     "\n",
     "# cosmology for the simulation\n",
@@ -70,10 +65,10 @@
     "cosmo = Cosmology.from_camb(pars)\n",
     "\n",
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
+    "zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n",
     "\n",
     "# linear window function for shells\n",
-    "ws = glass.shells.linear_windows(zb)\n",
+    "ws = glass.linear_windows(zb)\n",
     "\n",
     "# load the angular matter power spectra previously computed with CAMB\n",
     "cls = np.load(\"../1-basic/cls.npy\")"
@@ -102,10 +97,10 @@
    "source": [
     "# compute Gaussian cls for lognormal fields for 3 correlated shells\n",
     "# putting nside here means that the HEALPix pixel window function is applied\n",
-    "gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
+    "gls = glass.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
     "\n",
     "# generator for lognormal matter fields\n",
-    "matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)"
+    "matter = glass.generate_lognormal(gls, nside, ncorr=3)"
    ]
   },
   {
@@ -130,7 +125,7 @@
    "outputs": [],
    "source": [
     "# this will compute the convergence field iteratively\n",
-    "convergence = glass.lensing.MultiPlaneConvergence(cosmo)"
+    "convergence = glass.MultiPlaneConvergence(cosmo)"
    ]
   },
   {
@@ -167,7 +162,7 @@
     "dndz *= n_arcmin2 / np.trapz(dndz, z)\n",
     "\n",
     "# distribute dN/dz over the radial window functions\n",
-    "ngal = glass.shells.partition(z, dndz, ws)"
+    "ngal = glass.partition(z, dndz, ws)"
    ]
   },
   {
@@ -213,21 +208,16 @@
     "    # compute the lensing maps for this shell\n",
     "    convergence.add_window(delta_i, ws[i])\n",
     "    kappa_i = convergence.kappa\n",
-    "    gamm1_i, gamm2_i = glass.lensing.shear_from_convergence(kappa_i)\n",
+    "    gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n",
     "\n",
     "    # generate galaxy positions uniformly over the sphere\n",
-    "    for gal_lon, gal_lat, gal_count in glass.points.uniform_positions(ngal[i]):\n",
+    "    for gal_lon, gal_lat, gal_count in glass.uniform_positions(ngal[i]):\n",
     "        # generate galaxy ellipticities from the chosen distribution\n",
-    "        gal_eps = glass.shapes.ellipticity_intnorm(gal_count, sigma_e)\n",
+    "        gal_eps = glass.ellipticity_intnorm(gal_count, sigma_e)\n",
     "\n",
     "        # apply the shear fields to the ellipticities\n",
-    "        gal_she = glass.galaxies.galaxy_shear(\n",
-    "            gal_lon,\n",
-    "            gal_lat,\n",
-    "            gal_eps,\n",
-    "            kappa_i,\n",
-    "            gamm1_i,\n",
-    "            gamm2_i,\n",
+    "        gal_she = glass.galaxy_shear(\n",
+    "            gal_lon, gal_lat, gal_eps, kappa_i, gamm1_i, gamm2_i\n",
     "        )\n",
     "\n",
     "        # map the galaxy shears to a HEALPix map; this is opaque but works\n",
diff --git a/examples/2-advanced/stage_4_galaxies.ipynb b/examples/2-advanced/stage_4_galaxies.ipynb
index 598d4bed..7f2269bd 100644
--- a/examples/2-advanced/stage_4_galaxies.ipynb
+++ b/examples/2-advanced/stage_4_galaxies.ipynb
@@ -45,14 +45,8 @@
     "from cosmology import Cosmology\n",
     "\n",
     "# GLASS modules: cosmology and everything in the glass namespace\n",
+    "import glass\n",
     "import glass.ext.camb\n",
-    "import glass.fields\n",
-    "import glass.galaxies\n",
-    "import glass.lensing\n",
-    "import glass.observations\n",
-    "import glass.points\n",
-    "import glass.shapes\n",
-    "import glass.shells\n",
     "\n",
     "# cosmology for the simulation\n",
     "h = 0.7\n",
@@ -96,20 +90,20 @@
    "outputs": [],
    "source": [
     "# shells of 200 Mpc in comoving distance spacing\n",
-    "zb = glass.shells.distance_grid(cosmo, 0.0, 3.0, dx=200.0)\n",
+    "zb = glass.distance_grid(cosmo, 0.0, 3.0, dx=200.0)\n",
     "\n",
     "# linear window functions for shells\n",
-    "ws = glass.shells.linear_windows(zb)\n",
+    "ws = glass.linear_windows(zb)\n",
     "\n",
     "# compute the angular matter power spectra of the shells with CAMB\n",
     "cls = glass.ext.camb.matter_cls(pars, lmax, ws)\n",
     "\n",
     "# compute Gaussian cls for lognormal fields for 3 correlated shells\n",
     "# putting nside here means that the HEALPix pixel window function is applied\n",
-    "gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
+    "gls = glass.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)\n",
     "\n",
     "# generator for lognormal matter fields\n",
-    "matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)"
+    "matter = glass.generate_lognormal(gls, nside, ncorr=3)"
    ]
   },
   {
@@ -134,7 +128,7 @@
    "outputs": [],
    "source": [
     "# this will compute the convergence field iteratively\n",
-    "convergence = glass.lensing.MultiPlaneConvergence(cosmo)"
+    "convergence = glass.MultiPlaneConvergence(cosmo)"
    ]
   },
   {
@@ -163,15 +157,15 @@
     "\n",
     "# true redshift distribution following a Smail distribution\n",
     "z = np.arange(0.0, 3.0, 0.01)\n",
-    "dndz = glass.observations.smail_nz(z, z_mode=0.9, alpha=2.0, beta=1.5)\n",
+    "dndz = glass.smail_nz(z, z_mode=0.9, alpha=2.0, beta=1.5)\n",
     "dndz *= n_arcmin2\n",
     "\n",
     "# distribute dN/dz over the radial window functions\n",
-    "ngal = glass.shells.partition(z, dndz, ws)\n",
+    "ngal = glass.partition(z, dndz, ws)\n",
     "\n",
     "# compute tomographic redshift bin edges with equal density\n",
     "nbins = 10\n",
-    "zbins = glass.observations.equal_dens_zbins(z, dndz, nbins=nbins)\n",
+    "zbins = glass.equal_dens_zbins(z, dndz, nbins=nbins)\n",
     "\n",
     "# photometric redshift error\n",
     "sigma_z0 = 0.03\n",
@@ -220,7 +214,7 @@
     }
    ],
    "source": [
-    "vis = glass.observations.vmap_galactic_ecliptic(nside)\n",
+    "vis = glass.vmap_galactic_ecliptic(nside)\n",
     "\n",
     "# checking the mask:\n",
     "hp.mollview(vis, title=\"Stage IV Space Survey-like Mask\", unit=\"Visibility\")\n",
@@ -278,35 +272,27 @@
     "    # compute the lensing maps for this shell\n",
     "    convergence.add_window(delta_i, ws[i])\n",
     "    kappa_i = convergence.kappa\n",
-    "    gamm1_i, gamm2_i = glass.lensing.shear_from_convergence(kappa_i)\n",
+    "    gamm1_i, gamm2_i = glass.shear_from_convergence(kappa_i)\n",
     "\n",
     "    # generate galaxy positions from the matter density contrast\n",
-    "    for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(\n",
-    "        ngal[i],\n",
-    "        delta_i,\n",
-    "        bias,\n",
-    "        vis,\n",
+    "    for gal_lon, gal_lat, gal_count in glass.positions_from_delta(\n",
+    "        ngal[i], delta_i, bias, vis\n",
     "    ):\n",
     "        # generate random redshifts over the given shell\n",
-    "        gal_z = glass.galaxies.redshifts(gal_count, ws[i])\n",
+    "        gal_z = glass.redshifts(gal_count, ws[i])\n",
     "\n",
     "        # generator photometric redshifts using a Gaussian model\n",
-    "        gal_phz = glass.galaxies.gaussian_phz(gal_z, sigma_z0)\n",
+    "        gal_phz = glass.gaussian_phz(gal_z, sigma_z0)\n",
     "\n",
     "        # attach tomographic bin IDs to galaxies, based on photometric redshifts\n",
     "        gal_zbin = np.digitize(gal_phz, np.unique(zbins)) - 1\n",
     "\n",
     "        # generate galaxy ellipticities from the chosen distribution\n",
-    "        gal_eps = glass.shapes.ellipticity_intnorm(gal_count, sigma_e)\n",
+    "        gal_eps = glass.ellipticity_intnorm(gal_count, sigma_e)\n",
     "\n",
     "        # apply the shear fields to the ellipticities\n",
-    "        gal_she = glass.galaxies.galaxy_shear(\n",
-    "            gal_lon,\n",
-    "            gal_lat,\n",
-    "            gal_eps,\n",
-    "            kappa_i,\n",
-    "            gamm1_i,\n",
-    "            gamm2_i,\n",
+    "        gal_she = glass.galaxy_shear(\n",
+    "            gal_lon, gal_lat, gal_eps, kappa_i, gamm1_i, gamm2_i\n",
     "        )\n",
     "\n",
     "        # make a mini-catalogue for the new rows\n",
@@ -365,7 +351,7 @@
    ],
    "source": [
     "# split dndz using the same Gaussian error model assumed in the sampling\n",
-    "tomo_nz = glass.observations.tomo_nz_gausserr(z, dndz, sigma_z0, zbins)\n",
+    "tomo_nz = glass.tomo_nz_gausserr(z, dndz, sigma_z0, zbins)\n",
     "\n",
     "# redshift distribution of tomographic bins & input distributions\n",
     "plt.figure()\n",
diff --git a/glass/__init__.py b/glass/__init__.py
index 71bcc518..5fdd889c 100644
--- a/glass/__init__.py
+++ b/glass/__init__.py
@@ -1,6 +1,71 @@
 """GLASS package."""
 
 import contextlib
+from importlib.metadata import PackageNotFoundError
 
-with contextlib.suppress(ModuleNotFoundError):
-    from ._version import __version__, __version_tuple__  # noqa: F401
+with contextlib.suppress(PackageNotFoundError):
+    from ._version import __version__, __version_tuple__
+
+from glass.fields import (
+    cls2cov,
+    effective_cls,
+    gaussian_gls,
+    generate_gaussian,
+    generate_lognormal,
+    getcl,
+    iternorm,
+    lognormal_gls,
+    multalm,
+    transform_cls,
+)
+from glass.galaxies import (
+    galaxy_shear,
+    gaussian_phz,
+    redshifts,
+    redshifts_from_nz,
+)
+from glass.lensing import (
+    MultiPlaneConvergence,
+    deflect,
+    from_convergence,
+    multi_plane_matrix,
+    multi_plane_weights,
+    shear_from_convergence,
+)
+from glass.observations import (
+    equal_dens_zbins,
+    fixed_zbins,
+    gaussian_nz,
+    smail_nz,
+    tomo_nz_gausserr,
+    vmap_galactic_ecliptic,
+)
+from glass.points import (
+    effective_bias,
+    linear_bias,
+    loglinear_bias,
+    position_weights,
+    positions_from_delta,
+    uniform_positions,
+)
+from glass.shapes import (
+    ellipticity_gaussian,
+    ellipticity_intnorm,
+    ellipticity_ryden04,
+    triaxial_axis_ratio,
+)
+from glass.shells import (
+    RadialWindow,
+    combine,
+    cubic_windows,
+    density_weight,
+    distance_grid,
+    distance_weight,
+    linear_windows,
+    partition,
+    redshift_grid,
+    restrict,
+    tophat_windows,
+    volume_weight,
+)
+from glass.user import load_cls, save_cls, write_catalog
diff --git a/glass/fields.py b/glass/fields.py
index 7929df00..94089c36 100644
--- a/glass/fields.py
+++ b/glass/fields.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Random fields (:mod:`glass.fields`).
-====================================
+Random fields.
+=============
 
-.. currentmodule:: glass.fields
+.. currentmodule:: glass
 
-The :mod:`glass.fields` module provides functionality for simulating random
+The following functions provide functionality for simulating random
 fields on the sphere.  This is done in the form of HEALPix maps.
 
 Functions
diff --git a/glass/galaxies.py b/glass/galaxies.py
index 6e096048..3d834a02 100644
--- a/glass/galaxies.py
+++ b/glass/galaxies.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Galaxies (:mod:`glass.galaxies`).
-================================
+Galaxies.
+========
 
-.. currentmodule:: glass.galaxies
+.. currentmodule:: glass
 
-The :mod:`glass.galaxies` module provides functionality for simulating galaxies
+The following functions provide functionality for simulating galaxies
 as typically observed in a cosmological galaxy survey.
 
 Functions
@@ -23,15 +23,15 @@
 
 from typing import TYPE_CHECKING
 
-import healpix
-import numpy as np
-
-from .core.array import broadcast_leading_axes, cumtrapz
-
 if TYPE_CHECKING:
     from numpy.typing import ArrayLike
 
-    from .shells import RadialWindow
+    from glass.shells import RadialWindow
+
+import healpix
+import numpy as np
+
+from glass.core.array import broadcast_leading_axes, cumtrapz
 
 
 def redshifts(
@@ -51,7 +51,7 @@ def redshifts(
     n : int or array_like
         Number of redshifts to sample.  If an array is given, the
         results are concatenated.
-    w : :class:`~glass.shells.RadialWindow`
+    w : :class:`~glass.RadialWindow`
         Radial window function.
     rng : :class:`~numpy.random.Generator`, optional
         Random number generator.  If not given, a default RNG is used.
@@ -236,7 +236,7 @@ def gaussian_phz(
 
     See Also
     --------
-    glass.observations.tomo_nz_gausserr :
+    glass.tomo_nz_gausserr :
         Create tomographic redshift distributions assuming the same
         model.
 
diff --git a/glass/lensing.py b/glass/lensing.py
index f03075f8..a9b87bff 100644
--- a/glass/lensing.py
+++ b/glass/lensing.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Lensing (:mod:`glass.lensing`).
-==============================
+Lensing.
+========
 
-.. currentmodule:: glass.lensing
+.. currentmodule:: glass
 
-The :mod:`glass.lensing` module provides functionality for simulating
+The following functions/classes provide functionality for simulating
 gravitational lensing by the matter distribution in the universe.
 
 Iterative lensing
@@ -45,7 +45,7 @@
 
     from cosmology import Cosmology
 
-    from .shells import RadialWindow
+    from glass.shells import RadialWindow
 
 
 def from_convergence(  # noqa: PLR0913
@@ -398,7 +398,7 @@ def multi_plane_weights(
     weights : array_like
         Relative weight of each shell.  The first axis must broadcast
         against the number of shells, and is normalised internally.
-    shells : list of :class:`~glass.shells.RadialWindow`
+    shells : list of :class:`~glass.RadialWindow`
         Window functions of the shells.
     cosmo : Cosmology
         Cosmology instance.
diff --git a/glass/observations.py b/glass/observations.py
index 46980626..530ead88 100644
--- a/glass/observations.py
+++ b/glass/observations.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Observations (:mod:`glass.observations`).
-========================================
+Observations.
+============
 
-.. currentmodule:: glass.observations
+.. currentmodule:: glass
 
-The :mod:`glass.observations` module provides functionality for simulating
+The following functions provide functionality for simulating
 observational effects of surveys.
 
 
@@ -36,7 +36,7 @@
 import healpy as hp
 import numpy as np
 
-from .core.array import cumtrapz
+from glass.core.array import cumtrapz
 
 if TYPE_CHECKING:
     from numpy.typing import ArrayLike
diff --git a/glass/points.py b/glass/points.py
index 00db9c6e..1bef536b 100644
--- a/glass/points.py
+++ b/glass/points.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Random points (:mod:`glass.points`).
-===================================
+Random points.
+=============
 
-.. currentmodule:: glass.points
+.. currentmodule:: glass
 
-The :mod:`glass.points` module provides functionality for simulating point
+The following functions provide functionality for simulating point
 processes on the sphere and sampling random positions.
 
 Sampling
@@ -34,8 +34,8 @@
 import healpix
 import numpy as np
 
-from .core.array import broadcast_first, broadcast_leading_axes, trapz_product
-from .core.constants import ARCMIN2_SPHERE
+from glass.core.array import broadcast_first, broadcast_leading_axes, trapz_product
+from glass.core.constants import ARCMIN2_SPHERE
 
 
 def effective_bias(z, bz, w):
@@ -50,7 +50,7 @@ def effective_bias(z, bz, w):
     ----------
     z, bz : array_like
         Redshifts and values of the bias function :math:`b(z)`.
-    w : :class:`~glass.shells.RadialWindow`
+    w : :class:`~glass.RadialWindow`
         The radial window function :math:`w(z)`.
 
     Returns
diff --git a/glass/shapes.py b/glass/shapes.py
index 05e5be53..d64d15f2 100644
--- a/glass/shapes.py
+++ b/glass/shapes.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Observed shapes (:mod:`glass.shapes`).
-=====================================
+Observed shapes.
+===============
 
-.. currentmodule:: glass.shapes
+.. currentmodule:: glass
 
-The :mod:`glass.shapes` module provides functionality for simulating the
+The following functions provide functionality for simulating the
 observed shapes of objects, such as e.g. galaxies.
 
 Ellipticity
diff --git a/glass/shells.py b/glass/shells.py
index b36c85a9..cd31f819 100644
--- a/glass/shells.py
+++ b/glass/shells.py
@@ -1,12 +1,12 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-Shells (:mod:`glass.shells`).
-============================
+Shells.
+======
 
-.. currentmodule:: glass.shells
+.. currentmodule:: glass
 
-The :mod:`glass.shells` module provides functions for the definition of
+The following functions provide functionality for the definition of
 matter shells, i.e. the radial discretisation of the light cone.
 
 
@@ -53,7 +53,7 @@
 
 import numpy as np
 
-from .core.array import ndinterp
+from glass.core.array import ndinterp
 
 if TYPE_CHECKING:
     from numpy.typing import ArrayLike
@@ -96,7 +96,7 @@ class RadialWindow(NamedTuple):
     immutable (however, the array entries may **not** be immutable; do
     not change them in place)::
 
-        >>> from glass.shells import RadialWindow
+        >>> from glass import RadialWindow
         >>> w1 = RadialWindow(..., ..., zeff=0.1)
         >>> w1.zeff = 0.15
         Traceback (most recent call last):
@@ -528,7 +528,7 @@ def partition_nnls(
     requires SciPy.
 
     """
-    from .core.algorithm import nnls
+    from glass.core.algorithm import nnls
 
     # make sure nothing breaks
     sumtol = max(sumtol, 1e-4)
diff --git a/glass/user.py b/glass/user.py
index acfae20d..e099cae8 100644
--- a/glass/user.py
+++ b/glass/user.py
@@ -1,13 +1,13 @@
 # author: Nicolas Tessore <n.tessore@ucl.ac.uk>
 # license: MIT
 """
-User utilities (:mod:`glass.user`).
-==================================
+User utilities.
+==============
 
-.. currentmodule:: glass.user
+.. currentmodule:: glass
 
-The :mod:`glass.user` module contains convenience functions for users of the
-library.
+The following functions/classes provide convenience functionality for users
+of the library.
 
 
 Input and Output
diff --git a/pyproject.toml b/pyproject.toml
index f0d743d2..8acc4db9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -95,10 +95,10 @@ src = [
     "glass",
 ]
 lint.ignore = [
-    "ANN001", # TEMPORARY
-    "ANN002", # TEMPORARY
-    "ANN003", # TEMPORARY
-    "ANN201", # TEMPORARY
+    "ANN001", # TODO: missing-type-function-argument
+    "ANN002", # TODO: missing-type-args
+    "ANN003", # TODO: missing-type-kwargs
+    "ANN201", # TODO: missing-return-type-undocumented-public-function
     "COM812", # missing-trailing-comma (ruff-format recommended)
     "D203", # one-blank-line-before-class
     "D212", # blank-line-before-class
@@ -121,7 +121,9 @@ lint.isort = {known-first-party = [
     "camb",
     "cosmology",
 ]}}
-lint.per-file-ignores = {"docs*" = [
+lint.per-file-ignores = {"__init__.py" = [
+    "F401", # unused-import
+], "docs*" = [
     "D100", # undocumented-public-module
     "INP001", # implicit-namespace-package,
 ], "examples*" = [