diff --git a/.flake8 b/.flake8 deleted file mode 100644 index e7866dee..00000000 --- a/.flake8 +++ /dev/null @@ -1,2 +0,0 @@ -[flake8] -ignore = E226,E501,E741 diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 00000000..82406740 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +# applied ruff-format - https://github.com/glass-dev/glass/pull/199 +c894d9c7d349e44ac7fd50ba35859695758bbb95 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index fa9e274f..c2692843 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -16,12 +16,13 @@ concurrency: jobs: - style: - name: Style - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - run: pipx run flake8 glass +# TODO: replace with a more general pre-commit linter +# style: +# name: Style +# runs-on: ubuntu-latest +# steps: +# - uses: actions/checkout@v4 +# - run: pipx run flake8 glass tests: name: Tests diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0b2e536a..24cd4ab1 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,16 +1,14 @@ repos: -- repo: https://github.com/pycqa/flake8 - rev: 7.1.1 - hooks: - - id: flake8 - additional_dependencies: - - flake8-print -- repo: https://github.com/pappasam/toml-sort - rev: v0.23.1 - hooks: - - id: toml-sort-fix - args: - - --all - - --in-place - - --spaces-indent-inline-array=4 - - --trailing-comma-inline-array + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.6.3 + hooks: + - id: ruff-format + - repo: https://github.com/pappasam/toml-sort + rev: v0.23.1 + hooks: + - id: toml-sort-fix + args: + - --all + - --in-place + - --spaces-indent-inline-array=4 + - --trailing-comma-inline-array diff --git a/docs/conf.py b/docs/conf.py index 038320a6..498c7440 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -10,9 +10,9 @@ import datetime from importlib import metadata -project = 'GLASS' -author = 'GLASS developers' -copyright = f'2022-{datetime.date.today().year} {author}' +project = "GLASS" +author = "GLASS developers" +copyright = f"2022-{datetime.date.today().year} {author}" version = metadata.version("glass") release = version @@ -23,21 +23,21 @@ # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ - 'numpydoc', - 'sphinx.ext.intersphinx', - 'sphinxcontrib.katex', - 'matplotlib.sphinxext.plot_directive', - 'nbsphinx', + "numpydoc", + "sphinx.ext.intersphinx", + "sphinxcontrib.katex", + "matplotlib.sphinxext.plot_directive", + "nbsphinx", "sphinx.ext.viewcode", ] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] # -- Options for HTML output ------------------------------------------------- @@ -45,15 +45,15 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'furo' +html_theme = "furo" # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] -html_logo = '_static/logo.png' -html_favicon = '_static/favicon.ico' +html_logo = "_static/logo.png" +html_favicon = "_static/favicon.ico" html_css_files = [] @@ -62,14 +62,14 @@ # This config value contains the locations and names of other projects that # should be linked to in this documentation. intersphinx_mapping = { - 'python': ('https://docs.python.org/3', None), - 'numpy': ('https://numpy.org/doc/stable/', None), + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), } # -- autodoc ----------------------------------------------------------------- -autodoc_typehints = 'none' +autodoc_typehints = "none" # -- numpydoc ---------------------------------------------------------------- @@ -89,12 +89,12 @@ plot_html_show_formats = False # File formats to generate (default: ['png', 'hires.png', 'pdf']). -plot_formats = [('svg', 150), ('png', 150)] +plot_formats = [("svg", 150), ("png", 150)] # A dictionary containing any non-standard rcParams that should be applied # before each plot (default: {}). plot_rcparams = { - 'axes.facecolor': (1.0, 1.0, 1.0, 1.0), - 'savefig.facecolor': (1.0, 1.0, 1.0, 0.5), - 'savefig.transparent': False, + "axes.facecolor": (1.0, 1.0, 1.0, 1.0), + "savefig.facecolor": (1.0, 1.0, 1.0, 0.5), + "savefig.transparent": False, } diff --git a/examples/1-basic/density.ipynb b/examples/1-basic/density.ipynb index 9cf61845..a21cf90e 100644 --- a/examples/1-basic/density.ipynb +++ b/examples/1-basic/density.ipynb @@ -58,20 +58,21 @@ "nside = lmax = 128\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "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.)\n", + "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n", "\n", "# linear radial window functions\n", "ws = glass.shells.linear_windows(zb)\n", "\n", "# load the angular matter power spectra previously computed with CAMB\n", - "cls = np.load('cls.npy')" + "cls = np.load(\"cls.npy\")" ] }, { @@ -124,7 +125,7 @@ "outputs": [], "source": [ "# constant galaxy density distribution\n", - "z = np.linspace(0., 1., 100)\n", + "z = np.linspace(0.0, 1.0, 100)\n", "dndz = np.full_like(z, 0.01)\n", "\n", "# distribute the dN/dz over the linear window functions\n", @@ -158,22 +159,24 @@ "source": [ "# make a cube for galaxy number in redshift\n", "zcub = np.linspace(-zb[-1], zb[-1], 21)\n", - "cube = np.zeros((zcub.size-1,)*3)\n", + "cube = np.zeros((zcub.size - 1,) * 3)\n", "\n", "# simulate and add galaxies in each matter shell to cube\n", "for i, delta_i in enumerate(matter):\n", - "\n", " # simulate positions from matter density\n", - " for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(ngal[i], delta_i):\n", - "\n", + " for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(\n", + " ngal[i], delta_i\n", + " ):\n", " # sample redshifts uniformly in shell\n", " gal_z = glass.galaxies.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", - " z2 = gal_z*np.sin(np.deg2rad(gal_lon))*np.cos(np.deg2rad(gal_lat))\n", - " z3 = gal_z*np.sin(np.deg2rad(gal_lat))\n", - " (i, j, k), c = np.unique(np.searchsorted(zcub[1:], [z1, z2, z3]), axis=1, return_counts=True)\n", + " z1 = gal_z * np.cos(np.deg2rad(gal_lon)) * np.cos(np.deg2rad(gal_lat))\n", + " z2 = gal_z * np.sin(np.deg2rad(gal_lon)) * np.cos(np.deg2rad(gal_lat))\n", + " z3 = gal_z * np.sin(np.deg2rad(gal_lat))\n", + " (i, j, k), c = np.unique(\n", + " np.searchsorted(zcub[1:], [z1, z2, z3]), axis=1, return_counts=True\n", + " )\n", " cube[i, j, k] += c" ] }, @@ -217,19 +220,28 @@ ], "source": [ "# positions of grid cells of the cube\n", - "z = (zcub[:-1] + zcub[1:])/2\n", + "z = (zcub[:-1] + zcub[1:]) / 2\n", "z1, z2, z3 = np.meshgrid(z, z, z)\n", "\n", "# plot the galaxy distribution in pseudo-3D\n", "fig = plt.figure()\n", - "ax = fig.add_subplot(111, projection='3d', proj_type='ortho')\n", + "ax = fig.add_subplot(111, projection=\"3d\", proj_type=\"ortho\")\n", "norm = LogNorm(vmin=np.min(cube[cube > 0]), vmax=np.max(cube), clip=True)\n", - "for i in range(len(zcub)-1):\n", + "for i in range(len(zcub) - 1):\n", " v = norm(cube[..., i])\n", " c = plt.cm.inferno(v)\n", - " c[..., -1] = 0.2*v\n", - " ax.plot_surface(z1[..., i], z2[..., i], z3[..., i], rstride=1, cstride=1,\n", - " facecolors=c, linewidth=0, shade=False, antialiased=False)\n", + " c[..., -1] = 0.2 * v\n", + " ax.plot_surface(\n", + " z1[..., i],\n", + " z2[..., i],\n", + " z3[..., i],\n", + " rstride=1,\n", + " cstride=1,\n", + " facecolors=c,\n", + " linewidth=0,\n", + " shade=False,\n", + " antialiased=False,\n", + " )\n", "fig.tight_layout()\n", "plt.show()" ] diff --git a/examples/1-basic/lensing.ipynb b/examples/1-basic/lensing.ipynb index fad5bc3c..f7a6c3b5 100644 --- a/examples/1-basic/lensing.ipynb +++ b/examples/1-basic/lensing.ipynb @@ -62,20 +62,21 @@ "nside = lmax = 256\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "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.)\n", + "zb = glass.shells.distance_grid(cosmo, 0.0, 1.0, dx=200.0)\n", "\n", "# linear radial window functions\n", "ws = glass.shells.linear_windows(zb)\n", "\n", "# load the angular matter power spectra previously computed with CAMB\n", - "cls = np.load('cls.npy')" + "cls = np.load(\"cls.npy\")" ] }, { @@ -156,7 +157,7 @@ "# localised redshift distribution\n", "# the actual density per arcmin2 does not matter here, it is never used\n", "z = np.linspace(0.0, 1.0, 101)\n", - "dndz = np.exp(-(z - 0.5)**2/(0.1)**2)\n", + "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)" @@ -186,13 +187,12 @@ "outputs": [], "source": [ "# the integrated convergence and shear field over the redshift distribution\n", - "kappa_bar = np.zeros(12*nside**2)\n", - "gamm1_bar = np.zeros(12*nside**2)\n", - "gamm2_bar = np.zeros(12*nside**2)\n", + "kappa_bar = np.zeros(12 * nside**2)\n", + "gamm1_bar = np.zeros(12 * nside**2)\n", + "gamm2_bar = np.zeros(12 * nside**2)\n", "\n", "# main loop to simulate the matter fields iterative\n", "for i, delta_i in enumerate(matter):\n", - "\n", " # add lensing plane from the window function of this shell\n", " convergence.add_window(delta_i, ws[i])\n", "\n", @@ -255,26 +255,29 @@ ], "source": [ "# get the angular power spectra of the lensing maps\n", - "sim_cls = hp.anafast([kappa_bar, gamm1_bar, gamm2_bar],\n", - " pol=True, lmax=lmax, use_pixel_weights=True)\n", + "sim_cls = hp.anafast(\n", + " [kappa_bar, gamm1_bar, gamm2_bar], pol=True, lmax=lmax, use_pixel_weights=True\n", + ")\n", "\n", "# get the expected cls from CAMB\n", "pars.min_l = 1\n", "pars.set_for_lmax(lmax)\n", - "pars.SourceWindows = [camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type='lensing')]\n", + "pars.SourceWindows = [\n", + " camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type=\"lensing\")\n", + "]\n", "theory_cls = camb.get_results(pars).get_source_cls_dict(lmax=lmax, raw_cl=True)\n", "\n", "# get the HEALPix pixel window function, since the lensing fields have it\n", "pw = hp.pixwin(nside, lmax=lmax)\n", "\n", "# plot the realised and expected cls\n", - "l = np.arange(lmax+1)\n", - "plt.plot(l, sim_cls[0], '-k', lw=2, label='simulation')\n", - "plt.plot(l, theory_cls['W1xW1']*pw**2, '-r', lw=2, label='expectation')\n", - "plt.xscale('symlog', linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", - "plt.yscale('symlog', linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", - "plt.xlabel(r'angular mode number $l$')\n", - "plt.ylabel(r'angular power spectrum $C_l^{\\kappa\\kappa}$')\n", + "l = np.arange(lmax + 1)\n", + "plt.plot(l, sim_cls[0], \"-k\", lw=2, label=\"simulation\")\n", + "plt.plot(l, theory_cls[\"W1xW1\"] * pw**2, \"-r\", lw=2, label=\"expectation\")\n", + "plt.xscale(\"symlog\", linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", + "plt.yscale(\"symlog\", linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", + "plt.xlabel(r\"angular mode number $l$\")\n", + "plt.ylabel(r\"angular power spectrum $C_l^{\\kappa\\kappa}$\")\n", "plt.legend(frameon=False)\n", "plt.show()" ] diff --git a/examples/1-basic/matter.ipynb b/examples/1-basic/matter.ipynb index 1a15d084..a7636e39 100644 --- a/examples/1-basic/matter.ipynb +++ b/examples/1-basic/matter.ipynb @@ -57,17 +57,18 @@ "lmax = 1000\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology.from_camb(pars)\n", "\n", "# shells of 200 Mpc in comoving distance spacing\n", - "zb = glass.shells.distance_grid(cosmo, 0., 1., dx=200.)\n", + "zb = glass.shells.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", + "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", @@ -118,27 +119,28 @@ "source": [ "# make a 2d grid in redshift\n", "n = 2000\n", - "zend = 1.05*zb[-1]\n", - "x, y = np.mgrid[-zend:zend:1j*n, -zend:zend:1j*n]\n", + "zend = 1.05 * zb[-1]\n", + "x, y = np.mgrid[-zend : zend : 1j * n, -zend : zend : 1j * n]\n", "z = np.hypot(x, y)\n", "grid = np.full(z.shape, np.nan)\n", "\n", "# set up the plot\n", "ax = plt.subplot(111)\n", - "ax.axis('off')\n", + "ax.axis(\"off\")\n", "\n", "# simulate and project an annulus of each matter shell onto the grid\n", "for i, delta_i in enumerate(matter):\n", - " zmin, zmax = zb[i], zb[i+1]\n", + " zmin, zmax = zb[i], zb[i + 1]\n", " g = (zmin <= z) & (z < zmax)\n", - " zg = np.sqrt(1 - (z[g]/zmax)**2)\n", - " theta, phi = hp.vec2ang(np.transpose([x[g]/zmax, y[g]/zmax, zg]))\n", + " zg = np.sqrt(1 - (z[g] / zmax) ** 2)\n", + " theta, phi = hp.vec2ang(np.transpose([x[g] / zmax, y[g] / zmax, zg]))\n", " grid[g] = hp.get_interp_val(delta_i, theta, phi)\n", - " ax.add_patch(plt.Circle((0, 0), zmax/zend, fc='none', ec='k', lw=0.5, alpha=0.5, zorder=1))\n", + " ax.add_patch(\n", + " plt.Circle((0, 0), zmax / zend, fc=\"none\", ec=\"k\", lw=0.5, alpha=0.5, zorder=1)\n", + " )\n", "\n", "# show the grid of shells\n", - "ax.imshow(grid, extent=[-1, 1, -1, 1], zorder=0,\n", - " cmap='bwr', vmin=-2, vmax=2)\n", + "ax.imshow(grid, extent=[-1, 1, -1, 1], zorder=0, cmap=\"bwr\", vmin=-2, vmax=2)\n", "\n", "# show the resulting plot\n", "plt.show()" diff --git a/examples/1-basic/photoz.ipynb b/examples/1-basic/photoz.ipynb index d4ca56a4..20d7b1df 100644 --- a/examples/1-basic/photoz.ipynb +++ b/examples/1-basic/photoz.ipynb @@ -136,9 +136,9 @@ ], "source": [ "plt.figure(figsize=(5, 5))\n", - "plt.plot(ztrue, zphot, '+k', ms=3, alpha=0.1)\n", - "plt.xlabel(r'$z_{\\rm true}$', size=12)\n", - "plt.ylabel(r'$z_{\\rm phot}$', size=12)\n", + "plt.plot(ztrue, zphot, \"+k\", ms=3, alpha=0.1)\n", + "plt.xlabel(r\"$z_{\\rm true}$\", size=12)\n", + "plt.ylabel(r\"$z_{\\rm phot}$\", size=12)\n", "plt.show()" ] }, @@ -212,14 +212,19 @@ ], "source": [ "tomo_nz = glass.observations.tomo_nz_gausserr(z, dndz, phz_sigma_0, zbins)\n", - "tomo_nz *= ARCMIN2_SPHERE*(z[-1] - z[0])/40\n", + "tomo_nz *= ARCMIN2_SPHERE * (z[-1] - z[0]) / 40\n", "\n", "for (z1, z2), nz in zip(zbins, tomo_nz):\n", - " plt.hist(ztrue[(z1 <= zphot) & (zphot < z2)], bins=40, range=(z[0], z[-1]),\n", - " histtype='stepfilled', alpha=0.5)\n", - " plt.plot(z, nz, '-k', lw=1, alpha=0.5)\n", - "plt.xlabel('true redshift $z$')\n", - "plt.ylabel('number of galaxies')\n", + " plt.hist(\n", + " ztrue[(z1 <= zphot) & (zphot < z2)],\n", + " bins=40,\n", + " range=(z[0], z[-1]),\n", + " histtype=\"stepfilled\",\n", + " alpha=0.5,\n", + " )\n", + " plt.plot(z, nz, \"-k\", lw=1, alpha=0.5)\n", + "plt.xlabel(\"true redshift $z$\")\n", + "plt.ylabel(\"number of galaxies\")\n", "plt.show()" ] } diff --git a/examples/1-basic/shells.ipynb b/examples/1-basic/shells.ipynb index 614bad7a..79691876 100644 --- a/examples/1-basic/shells.ipynb +++ b/examples/1-basic/shells.ipynb @@ -56,14 +56,15 @@ "lmax = 1000\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "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.)\n", + "zgrid = glass.shells.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", @@ -146,7 +147,7 @@ }, "outputs": [], "source": [ - "np.save('cls.npy', cls)" + "np.save(\"cls.npy\", cls)" ] } ], diff --git a/examples/2-advanced/cosmic_shear.ipynb b/examples/2-advanced/cosmic_shear.ipynb index 1d6c8c80..58aeb4a6 100644 --- a/examples/2-advanced/cosmic_shear.ipynb +++ b/examples/2-advanced/cosmic_shear.ipynb @@ -60,20 +60,21 @@ "nside = lmax = 256\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology.from_camb(pars)\n", "\n", "# shells of 200 Mpc in comoving distance spacing\n", - "zb = glass.shells.distance_grid(cosmo, 0., 1., dx=200.)\n", + "zb = glass.shells.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", "\n", "# load the angular matter power spectra previously computed with CAMB\n", - "cls = np.load('../1-basic/cls.npy')" + "cls = np.load(\"../1-basic/cls.npy\")" ] }, { @@ -160,8 +161,8 @@ "\n", "# localised redshift distribution with the given density\n", "z = np.arange(0.0, 2.0, 0.01)\n", - "dndz = np.exp(-(z - 0.5)**2/(0.1)**2)\n", - "dndz *= n_arcmin2/np.trapz(dndz, z)\n", + "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n", + "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)" @@ -197,7 +198,7 @@ "outputs": [], "source": [ "# number of HEALPix pixels in the maps\n", - "npix = 12*nside**2\n", + "npix = 12 * nside**2\n", "\n", "# map for galaxy numbers\n", "num = np.zeros(npix)\n", @@ -207,7 +208,6 @@ "\n", "# simulate the matter fields in the main loop\n", "for i, delta_i in enumerate(matter):\n", - "\n", " # compute the lensing maps for this shell\n", " convergence.add_window(delta_i, ws[i])\n", " kappa_i = convergence.kappa\n", @@ -215,13 +215,13 @@ "\n", " # generate galaxy positions uniformly over the sphere\n", " for gal_lon, gal_lat, gal_count in glass.points.uniform_positions(ngal[i]):\n", - "\n", " # generate galaxy ellipticities from the chosen distribution\n", " gal_eps = glass.shapes.ellipticity_intnorm(gal_count, sigma_e)\n", "\n", " # apply the shear fields to the ellipticities\n", - " gal_she = glass.galaxies.galaxy_shear(gal_lon, gal_lat, gal_eps,\n", - " kappa_i, gamm1_i, gamm2_i)\n", + " gal_she = glass.galaxies.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", " gal_pix = hp.ang2pix(nside, gal_lon, gal_lat, lonlat=True)\n", @@ -272,7 +272,7 @@ ], "source": [ "# compute the expected number of galaxies in each pixel\n", - "nbar = ARCMIN2_SPHERE/npix*n_arcmin2\n", + "nbar = ARCMIN2_SPHERE / npix * n_arcmin2\n", "\n", "# normalise the maps by the expected number of galaxies in each pixel\n", "she /= nbar\n", @@ -282,34 +282,36 @@ "cls = hp.anafast([num, she.real, she.imag], pol=True, lmax=lmax, use_pixel_weights=True)\n", "\n", "# get the theory cls from CAMB\n", - "pars.NonLinear = 'NonLinear_both'\n", + "pars.NonLinear = \"NonLinear_both\"\n", "pars.Want_CMB = False\n", "pars.min_l = 1\n", "pars.set_for_lmax(lmax)\n", - "pars.SourceWindows = [camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type='lensing')]\n", + "pars.SourceWindows = [\n", + " camb.sources.SplinedSourceWindow(z=z, W=dndz, source_type=\"lensing\")\n", + "]\n", "theory_cls = camb.get_results(pars).get_source_cls_dict(lmax=lmax, raw_cl=True)\n", "\n", "# factor transforming convergence to shear\n", - "l = np.arange(lmax+1)\n", - "fl = (l+2)*(l+1)*l*(l-1)/np.clip(l**2*(l+1)**2, 1, None)\n", + "l = np.arange(lmax + 1)\n", + "fl = (l + 2) * (l + 1) * l * (l - 1) / np.clip(l**2 * (l + 1) ** 2, 1, None)\n", "\n", "# the noise level from discrete observations with shape noise\n", - "nl = 4*np.pi/(nbar*npix)*sigma_e**2 * (l >= 2)\n", + "nl = 4 * np.pi / (nbar * npix) * sigma_e**2 * (l >= 2)\n", "\n", "# mixing matrix for uniform distribution of points\n", - "mm = (1 - 1/(nbar*npix))*np.eye(lmax+1, lmax+1) + (l+1/2)/(nbar*npix)\n", + "mm = (1 - 1 / (nbar * npix)) * np.eye(lmax + 1, lmax + 1) + (l + 1 / 2) / (nbar * npix)\n", "mm[:2, :] = mm[:, :2] = 0\n", "\n", "# the shear pixel window function for HEALPix\n", "_, pw = hp.pixwin(nside, lmax=lmax, pol=True)\n", "\n", "# plot the realised and expected cls\n", - "plt.plot(l, cls[1] - nl, '-k', lw=2, label='simulation')\n", - "plt.plot(l, pw**2 * mm@(fl*theory_cls['W1xW1']), '-r', lw=2, label='expectation')\n", - "plt.xscale('symlog', linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", - "plt.yscale('symlog', linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", - "plt.xlabel('angular mode number $l$')\n", - "plt.ylabel('angular power spectrum $C_l^{EE}$')\n", + "plt.plot(l, cls[1] - nl, \"-k\", lw=2, label=\"simulation\")\n", + "plt.plot(l, pw**2 * mm @ (fl * theory_cls[\"W1xW1\"]), \"-r\", lw=2, label=\"expectation\")\n", + "plt.xscale(\"symlog\", linthresh=10, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", + "plt.yscale(\"symlog\", linthresh=1e-9, linscale=0.5, subs=[2, 3, 4, 5, 6, 7, 8, 9])\n", + "plt.xlabel(\"angular mode number $l$\")\n", + "plt.ylabel(\"angular power spectrum $C_l^{EE}$\")\n", "plt.legend(frameon=False)\n", "plt.tight_layout()\n", "plt.show()" diff --git a/examples/2-advanced/stage_4_galaxies.ipynb b/examples/2-advanced/stage_4_galaxies.ipynb index 57bb16a4..536d9d73 100644 --- a/examples/2-advanced/stage_4_galaxies.ipynb +++ b/examples/2-advanced/stage_4_galaxies.ipynb @@ -64,8 +64,9 @@ "nside = lmax = 256\n", "\n", "# set up CAMB parameters for matter angular power spectrum\n", - "pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,\n", - " NonLinear=camb.model.NonLinear_both)\n", + "pars = camb.set_params(\n", + " H0=100 * h, omch2=Oc * h**2, ombh2=Ob * h**2, NonLinear=camb.model.NonLinear_both\n", + ")\n", "\n", "# get the cosmology from CAMB\n", "cosmo = Cosmology.from_camb(pars)" @@ -93,7 +94,7 @@ "outputs": [], "source": [ "# shells of 200 Mpc in comoving distance spacing\n", - "zb = glass.shells.distance_grid(cosmo, 0., 3., dx=200.)\n", + "zb = glass.shells.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", @@ -220,7 +221,7 @@ "vis = glass.observations.vmap_galactic_ecliptic(nside)\n", "\n", "# checking the mask:\n", - "hp.mollview(vis, title='Stage IV Space Survey-like Mask', unit='Visibility')\n", + "hp.mollview(vis, title=\"Stage IV Space Survey-like Mask\", unit=\"Visibility\")\n", "plt.show()" ] }, @@ -257,21 +258,30 @@ ], "source": [ "# we will store the catalogue as a structured numpy array, initially empty\n", - "catalogue = np.empty(0, dtype=[('RA', float), ('DEC', float),\n", - " ('Z_TRUE', float), ('PHZ', float), ('ZBIN', int),\n", - " ('G1', float), ('G2', float)])\n", + "catalogue = np.empty(\n", + " 0,\n", + " dtype=[\n", + " (\"RA\", float),\n", + " (\"DEC\", float),\n", + " (\"Z_TRUE\", float),\n", + " (\"PHZ\", float),\n", + " (\"ZBIN\", int),\n", + " (\"G1\", float),\n", + " (\"G2\", float),\n", + " ],\n", + ")\n", "\n", "# simulate the matter fields in the main loop, and build up the catalogue\n", "for i, delta_i in enumerate(matter):\n", - "\n", " # 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", "\n", " # generate galaxy positions from the matter density contrast\n", - " for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(ngal[i], delta_i, bias, vis):\n", - "\n", + " for gal_lon, gal_lat, gal_count in glass.points.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", "\n", @@ -285,23 +295,24 @@ " gal_eps = glass.shapes.ellipticity_intnorm(gal_count, sigma_e)\n", "\n", " # apply the shear fields to the ellipticities\n", - " gal_she = glass.galaxies.galaxy_shear(gal_lon, gal_lat, gal_eps,\n", - " kappa_i, gamm1_i, gamm2_i)\n", + " gal_she = glass.galaxies.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", " rows = np.empty(gal_count, dtype=catalogue.dtype)\n", - " rows['RA'] = gal_lon\n", - " rows['DEC'] = gal_lat\n", - " rows['Z_TRUE'] = gal_z\n", - " rows['PHZ'] = gal_phz\n", - " rows['ZBIN'] = gal_zbin\n", - " rows['G1'] = gal_she.real\n", - " rows['G2'] = gal_she.imag\n", + " rows[\"RA\"] = gal_lon\n", + " rows[\"DEC\"] = gal_lat\n", + " rows[\"Z_TRUE\"] = gal_z\n", + " rows[\"PHZ\"] = gal_phz\n", + " rows[\"ZBIN\"] = gal_zbin\n", + " rows[\"G1\"] = gal_she.real\n", + " rows[\"G2\"] = gal_she.imag\n", "\n", " # add the new rows to the catalogue\n", " catalogue = np.append(catalogue, rows)\n", "\n", - "print(f'Total number of galaxies sampled: {len(catalogue):,}')" + "print(f\"Total number of galaxies sampled: {len(catalogue):,}\")" ] }, { @@ -348,15 +359,23 @@ "\n", "# redshift distribution of tomographic bins & input distributions\n", "plt.figure()\n", - "plt.title('redshifts in catalogue')\n", - "plt.ylabel('dN/dz - normalised')\n", - "plt.xlabel('z')\n", + "plt.title(\"redshifts in catalogue\")\n", + "plt.ylabel(\"dN/dz - normalised\")\n", + "plt.xlabel(\"z\")\n", "for i in range(nbins):\n", - " in_bin = (catalogue['ZBIN'] == i)\n", - " plt.hist(catalogue['Z_TRUE'][in_bin], histtype='stepfilled', edgecolor='none', alpha=0.5, bins=50, density=1, label=f'cat. bin {i}')\n", + " in_bin = catalogue[\"ZBIN\"] == i\n", + " plt.hist(\n", + " catalogue[\"Z_TRUE\"][in_bin],\n", + " histtype=\"stepfilled\",\n", + " edgecolor=\"none\",\n", + " alpha=0.5,\n", + " bins=50,\n", + " density=1,\n", + " label=f\"cat. bin {i}\",\n", + " )\n", "for i in range(nbins):\n", - " plt.plot(z, (tomo_nz[i]/n_arcmin2)*nbins, alpha=0.5, label=f'inp. bin {i}')\n", - "plt.plot(z, dndz/n_arcmin2*nbins, ls='--', c='k')\n", + " plt.plot(z, (tomo_nz[i] / n_arcmin2) * nbins, alpha=0.5, label=f\"inp. bin {i}\")\n", + "plt.plot(z, dndz / n_arcmin2 * nbins, ls=\"--\", c=\"k\")\n", "plt.legend(ncol=2)\n", "plt.show()" ] diff --git a/glass/core/__init__.py b/glass/core/__init__.py index 7e3ee958..f91e31a6 100644 --- a/glass/core/__init__.py +++ b/glass/core/__init__.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Core functions (:mod:`glass.core`) ================================== @@ -9,4 +9,4 @@ The :mod:`glass.core` module contains core functionality for developing GLASS modules. -''' +""" diff --git a/glass/core/algorithm.py b/glass/core/algorithm.py index c5fb5a32..cb8645ff 100644 --- a/glass/core/algorithm.py +++ b/glass/core/algorithm.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -'''core module for algorithms''' +"""core module for algorithms""" from __future__ import annotations @@ -59,10 +59,10 @@ def nnls( ap = a[:, p] xp = x[p] sp = np.linalg.solve(ap.T @ ap, b @ ap) - t = (sp <= 0) + t = sp <= 0 if not np.any(t): break - alpha = -np.min(xp[t]/(xp[t] - sp[t])) + alpha = -np.min(xp[t] / (xp[t] - sp[t])) x[p] += alpha * (sp - xp) p[x <= 0] = False x[p] = sp diff --git a/glass/core/array.py b/glass/core/array.py index 09784e21..2663751a 100644 --- a/glass/core/array.py +++ b/glass/core/array.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -'''module for array utilities''' +"""module for array utilities""" import numpy as np from functools import partial @@ -15,7 +15,7 @@ def broadcast_first(*arrays): def broadcast_leading_axes(*args): - '''Broadcast all but the last N axes. + """Broadcast all but the last N axes. Returns the shape of the broadcast dimensions, and all input arrays with leading axes matching that shape. @@ -38,7 +38,7 @@ def broadcast_leading_axes(*args): >>> c.shape (3, 4, 5, 6) - ''' + """ shapes, trails = [], [] for a, n in args: @@ -52,17 +52,19 @@ def broadcast_leading_axes(*args): def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None): - '''interpolate multi-dimensional array over axis''' - return np.apply_along_axis(partial(np.interp, x, xp), axis, fp, - left=left, right=right, period=period) + """interpolate multi-dimensional array over axis""" + return np.apply_along_axis( + partial(np.interp, x, xp), axis, fp, left=left, right=right, period=period + ) def trapz_product(f, *ff, axis=-1): - '''trapezoidal rule for a product of functions''' + """trapezoidal rule for a product of functions""" x, _ = f for x_, _ in ff: - x = np.union1d(x[(x >= x_[0]) & (x <= x_[-1])], - x_[(x_ >= x[0]) & (x_ <= x[-1])]) + x = np.union1d( + x[(x >= x_[0]) & (x <= x_[-1])], x_[(x_ >= x[0]) & (x_ <= x[-1])] + ) y = np.interp(x, *f) for f_ in ff: y *= np.interp(x, *f_) @@ -70,11 +72,11 @@ def trapz_product(f, *ff, axis=-1): def cumtrapz(f, x, dtype=None, out=None): - '''cumulative trapezoidal rule along last axis''' + """cumulative trapezoidal rule along last axis""" if out is None: out = np.empty_like(f, dtype=dtype) - np.cumsum((f[..., 1:] + f[..., :-1])/2*np.diff(x), axis=-1, out=out[..., 1:]) + np.cumsum((f[..., 1:] + f[..., :-1]) / 2 * np.diff(x), axis=-1, out=out[..., 1:]) out[..., 0] = 0 return out diff --git a/glass/core/constants.py b/glass/core/constants.py index 9357dff7..b70b5b57 100644 --- a/glass/core/constants.py +++ b/glass/core/constants.py @@ -1,9 +1,9 @@ # author: Nicolas Tessore # license: MIT -'''module for constants''' +"""module for constants""" PI = 3.1415926535897932384626433832795028841971693993751 -DEGREE2_SPHERE = 60**4//100/PI -ARCMIN2_SPHERE = 60**6//100/PI -ARCSEC2_SPHERE = 60**8//100/PI +DEGREE2_SPHERE = 60**4 // 100 / PI +ARCMIN2_SPHERE = 60**6 // 100 / PI +ARCSEC2_SPHERE = 60**8 // 100 / PI diff --git a/glass/ext/__init__.py b/glass/ext/__init__.py index 695c446e..8d471cf2 100644 --- a/glass/ext/__init__.py +++ b/glass/ext/__init__.py @@ -1,20 +1,20 @@ -'''Namespace loader for GLASS extension packages. +"""Namespace loader for GLASS extension packages. Uses the pkgutil namespace mechanism to find "ext" submodules of packages that provide a "glass" module. -''' +""" def _extend_path(path, name): import os.path from pkgutil import extend_path - _pkg, _, _mod = name.partition('.') + _pkg, _, _mod = name.partition(".") - return list(filter(os.path.isdir, - (os.path.join(p, _mod) - for p in extend_path(path, _pkg)))) + return list( + filter(os.path.isdir, (os.path.join(p, _mod) for p in extend_path(path, _pkg))) + ) __path__ = _extend_path(__path__, __name__) diff --git a/glass/fields.py b/glass/fields.py index 20defdb4..2f834fea 100644 --- a/glass/fields.py +++ b/glass/fields.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Random fields (:mod:`glass.fields`) =================================== @@ -24,7 +24,7 @@ .. autofunction:: getcl -''' +""" import warnings import numpy as np @@ -32,8 +32,7 @@ from gaussiancl import gaussiancl # typing -from typing import (Any, Union, Tuple, Generator, Optional, Sequence, Callable, - Iterable) +from typing import Any, Union, Tuple, Generator, Optional, Sequence, Callable, Iterable # types Array = np.ndarray @@ -44,9 +43,10 @@ Alms = np.ndarray -def iternorm(k: int, cov: Iterable[Array], size: Size = None - ) -> Generator[Iternorm, None, None]: - '''return the vector a and variance sigma^2 for iterative normal sampling''' +def iternorm( + k: int, cov: Iterable[Array], size: Size = None +) -> Generator[Iternorm, None, None]: + """return the vector a and variance sigma^2 for iterative normal sampling""" n: Tuple[int, ...] if size is None: @@ -59,7 +59,7 @@ def iternorm(k: int, cov: Iterable[Array], size: Size = None m = np.zeros((*n, k, k)) a = np.zeros((*n, k)) s = np.zeros((*n,)) - q = (*n, k+1) + q = (*n, k + 1) j = 0 if k > 0 else None for i, x in enumerate(cov): @@ -68,48 +68,54 @@ def iternorm(k: int, cov: Iterable[Array], size: Size = None try: x = np.broadcast_to(x, q) except ValueError: - raise TypeError(f'covariance row {i}: shape {x.shape} cannot be broadcast to {q}') from None + raise TypeError( + f"covariance row {i}: shape {x.shape} cannot be broadcast to {q}" + ) from None # only need to update matrix A if there are correlations if j is not None: # compute new entries of matrix A m[..., :, j] = 0 - m[..., j:j+1, :] = np.matmul(a[..., np.newaxis, :], m) + m[..., j : j + 1, :] = np.matmul(a[..., np.newaxis, :], m) m[..., j, j] = np.where(s != 0, -1, 0) - np.divide(m[..., j, :], -s[..., np.newaxis], where=(m[..., j, :] != 0), out=m[..., j, :]) + np.divide( + m[..., j, :], + -s[..., np.newaxis], + where=(m[..., j, :] != 0), + out=m[..., j, :], + ) # compute new vector a c = x[..., 1:, np.newaxis] - a = np.matmul(m[..., :j], c[..., k-j:, :]) - a += np.matmul(m[..., j:], c[..., :k-j, :]) + a = np.matmul(m[..., :j], c[..., k - j :, :]) + a += np.matmul(m[..., j:], c[..., : k - j, :]) a = a.reshape(*n, k) # next rolling index j = (j - 1) % k # compute new standard deviation - s = x[..., 0] - np.einsum('...i,...i', a, a) + s = x[..., 0] - np.einsum("...i,...i", a, a) if np.any(s < 0): - raise ValueError('covariance matrix is not positive definite') + raise ValueError("covariance matrix is not positive definite") s = np.sqrt(s) # yield the next index, vector a, and standard deviation s yield j, a, s -def cls2cov(cls: Cls, nl: int, nf: int, nc: int - ) -> Generator[Array, None, None]: - '''Return array of cls as a covariance matrix for iterative sampling.''' - cov = np.zeros((nl, nc+1)) +def cls2cov(cls: Cls, nl: int, nf: int, nc: int) -> Generator[Array, None, None]: + """Return array of cls as a covariance matrix for iterative sampling.""" + cov = np.zeros((nl, nc + 1)) end = 0 for j in range(nf): begin, end = end, end + j + 1 - for i, cl in enumerate(cls[begin:end][:nc+1]): + for i, cl in enumerate(cls[begin:end][: nc + 1]): if cl is None: cov[:, i] = 0 else: if i == 0 and np.any(np.less(cl, 0)): - raise ValueError('negative values in cl') + raise ValueError("negative values in cl") n = len(cl) cov[:n, i] = cl cov[n:, i] = 0 @@ -118,53 +124,60 @@ def cls2cov(cls: Cls, nl: int, nf: int, nc: int def multalm(alm: Alms, bl: Array, inplace: bool = False) -> Alms: - '''multiply alm by bl''' + """multiply alm by bl""" n = len(bl) if inplace: out = np.asanyarray(alm) else: out = np.copy(alm) for m in range(n): - out[m*n-m*(m-1)//2:(m+1)*n-m*(m+1)//2] *= bl[m:] + out[m * n - m * (m - 1) // 2 : (m + 1) * n - m * (m + 1) // 2] *= bl[m:] return out -def transform_cls(cls: Cls, tfm: ClTransform, pars: Tuple[Any, ...] = () - ) -> Cls: - '''Transform Cls to Gaussian Cls.''' +def transform_cls(cls: Cls, tfm: ClTransform, pars: Tuple[Any, ...] = ()) -> Cls: + """Transform Cls to Gaussian Cls.""" gls = [] for cl in cls: if cl is not None and len(cl) > 0: if cl[0] == 0: - monopole = 0. + monopole = 0.0 else: monopole = None gl, info, err, niter = gaussiancl(cl, tfm, pars, monopole=monopole) if info == 0: - warnings.warn('Gaussian cl did not converge, inexact transform') + warnings.warn("Gaussian cl did not converge, inexact transform") else: gl = [] gls.append(gl) return gls -def gaussian_gls(cls: Cls, *, lmax: Optional[int] = None, - ncorr: Optional[int] = None, nside: Optional[int] = None - ) -> Cls: - '''Compute Gaussian Cls for a Gaussian random field. +def gaussian_gls( + cls: Cls, + *, + lmax: Optional[int] = None, + ncorr: Optional[int] = None, + nside: Optional[int] = None, +) -> Cls: + """Compute Gaussian Cls for a Gaussian random field. Depending on the given arguments, this truncates the angular power spectra to ``lmax``, removes all but ``ncorr`` correlations between fields, and applies the HEALPix pixel window function of the given ``nside``. If no arguments are given, no action is performed. - ''' + """ if ncorr is not None: - n = int((2*len(cls))**0.5) - if n*(n+1)//2 != len(cls): - raise ValueError('length of cls array is not a triangle number') - cls = [cls[i*(i+1)//2+j] if j <= ncorr else [] for i in range(n) for j in range(i+1)] + n = int((2 * len(cls)) ** 0.5) + if n * (n + 1) // 2 != len(cls): + raise ValueError("length of cls array is not a triangle number") + cls = [ + cls[i * (i + 1) // 2 + j] if j <= ncorr else [] + for i in range(n) + for j in range(i + 1) + ] if nside is not None: pw = hp.pixwin(nside, lmax=lmax) @@ -173,26 +186,35 @@ def gaussian_gls(cls: Cls, *, lmax: Optional[int] = None, for cl in cls: if cl is not None and len(cl) > 0: if lmax is not None: - cl = cl[:lmax+1] + cl = cl[: lmax + 1] if nside is not None: n = min(len(cl), len(pw)) - cl = cl[:n] * pw[:n]**2 + cl = cl[:n] * pw[:n] ** 2 gls.append(cl) return gls -def lognormal_gls(cls: Cls, shift: float = 1., *, lmax: Optional[int] = None, - ncorr: Optional[int] = None, nside: Optional[int] = None - ) -> Cls: - '''Compute Gaussian Cls for a lognormal random field.''' +def lognormal_gls( + cls: Cls, + shift: float = 1.0, + *, + lmax: Optional[int] = None, + ncorr: Optional[int] = None, + nside: Optional[int] = None, +) -> Cls: + """Compute Gaussian Cls for a lognormal random field.""" gls = gaussian_gls(cls, lmax=lmax, ncorr=ncorr, nside=nside) - return transform_cls(gls, 'lognormal', (shift,)) + return transform_cls(gls, "lognormal", (shift,)) -def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, - rng: Optional[np.random.Generator] = None - ) -> Generator[Array, None, None]: - '''Iteratively sample Gaussian random fields from Cls. +def generate_gaussian( + gls: Cls, + nside: int, + *, + ncorr: Optional[int] = None, + rng: Optional[np.random.Generator] = None, +) -> Generator[Array, None, None]: + """Iteratively sample Gaussian random fields from Cls. A generator that iteratively samples HEALPix maps of Gaussian random fields with the given angular power spectra ``gls`` and resolution parameter @@ -213,7 +235,7 @@ def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, Missing entries can be set to ``None``. - ''' + """ # get the default RNG if not given if rng is None: @@ -221,7 +243,7 @@ def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, # number of gls and number of fields ngls = len(gls) - ngrf = int((2*ngls)**0.5) + ngrf = int((2 * ngls) ** 0.5) # number of correlated fields if not specified if ncorr is None: @@ -230,14 +252,14 @@ def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, # number of modes n = max((len(gl) for gl in gls if gl is not None), default=0) if n == 0: - raise ValueError('all gls are empty') + raise ValueError("all gls are empty") # generates the covariance matrix for the iterative sampler cov = cls2cov(gls, n, ngrf, ncorr) # working arrays for the iterative sampling - z = np.zeros(n*(n+1)//2, dtype=np.complex128) - y = np.zeros((n*(n+1)//2, ncorr), dtype=np.complex128) + z = np.zeros(n * (n + 1) // 2, dtype=np.complex128) + y = np.zeros((n * (n + 1) // 2, ncorr), dtype=np.complex128) # generate the conditional normal distribution for iterative sampling conditional_dist = iternorm(ncorr, cov, size=n) @@ -246,7 +268,7 @@ def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, for j, a, s in conditional_dist: # standard normal random variates for alm # sample real and imaginary parts, then view as complex number - rng.standard_normal(n*(n+1), np.float64, z.view(np.float64)) + rng.standard_normal(n * (n + 1), np.float64, z.view(np.float64)) # scale by standard deviation of the conditional distribution # variance is distributed over real and imaginary part @@ -269,19 +291,23 @@ def generate_gaussian(gls: Cls, nside: int, *, ncorr: Optional[int] = None, yield hp.alm2map(alm, nside, pixwin=False, pol=False, inplace=True) -def generate_lognormal(gls: Cls, nside: int, shift: float = 1., *, - ncorr: Optional[int] = None, - rng: Optional[np.random.Generator] = None - ) -> Generator[Array, None, None]: - '''Iterative sample lognormal random fields from Gaussian Cls.''' +def generate_lognormal( + gls: Cls, + nside: int, + shift: float = 1.0, + *, + ncorr: Optional[int] = None, + rng: Optional[np.random.Generator] = None, +) -> Generator[Array, None, None]: + """Iterative sample lognormal random fields from Gaussian Cls.""" for i, m in enumerate(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)): # compute the variance of the auto-correlation - gl = gls[i*(i+1)//2] + gl = gls[i * (i + 1) // 2] ell = np.arange(len(gl)) - var = np.sum((2*ell + 1)*gl)/(4*np.pi) + var = np.sum((2 * ell + 1) * gl) / (4 * np.pi) # fix mean of the Gaussian random field for lognormal transformation - m -= var/2 + m -= var / 2 # exponentiate values in place and subtract 1 in one operation np.expm1(m, out=m) @@ -317,10 +343,10 @@ def getcl(cls, i, j, lmax=None): """ if j > i: i, j = j, i - cl = cls[i*(i+1)//2+i-j] + cl = cls[i * (i + 1) // 2 + i - j] if lmax is not None: if len(cl) > lmax + 1: - cl = cl[:lmax+1] + cl = cl[: lmax + 1] else: cl = np.pad(cl, (0, lmax + 1 - len(cl))) return cl @@ -356,8 +382,8 @@ def effective_cls(cls, weights1, weights2=None, *, lmax=None): from itertools import combinations_with_replacement, product # this is the number of fields - n = int((2*len(cls))**0.5) - if n*(n+1)//2 != len(cls): + n = int((2 * len(cls)) ** 0.5) + if n * (n + 1) // 2 != len(cls): raise ValueError("length of cls is not a triangle number") # find lmax if not given @@ -381,7 +407,7 @@ def effective_cls(cls, weights1, weights2=None, *, lmax=None): pairs = product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:])) # create the output array: axes for all input axes plus lmax+1 - out = np.empty(shape1[1:] + shape2[1:] + (lmax+1,)) + out = np.empty(shape1[1:] + shape2[1:] + (lmax + 1,)) # helper that will grab the entire first column (i.e. shells) c = (slice(None),) diff --git a/glass/galaxies.py b/glass/galaxies.py index 69311c60..efecd8c2 100644 --- a/glass/galaxies.py +++ b/glass/galaxies.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Galaxies (:mod:`glass.galaxies`) ================================ @@ -17,7 +17,7 @@ .. autofunction:: galaxy_shear .. autofunction:: gaussian_phz -''' +""" from __future__ import annotations @@ -30,10 +30,10 @@ from .shells import RadialWindow -def redshifts(n: int | ArrayLike, w: RadialWindow, *, - rng: np.random.Generator | None = None - ) -> np.ndarray: - '''Sample redshifts from a radial window function. +def redshifts( + n: int | ArrayLike, w: RadialWindow, *, rng: np.random.Generator | None = None +) -> np.ndarray: + """Sample redshifts from a radial window function. This function samples *n* redshifts from a distribution that follows the given radial window function *w*. @@ -53,14 +53,18 @@ def redshifts(n: int | ArrayLike, w: RadialWindow, *, z : array_like Random redshifts following the radial window function. - ''' + """ return redshifts_from_nz(n, w.za, w.wa, rng=rng) -def redshifts_from_nz(count: int | ArrayLike, z: ArrayLike, nz: ArrayLike, *, - rng: np.random.Generator | None = None - ) -> np.ndarray: - '''Generate galaxy redshifts from a source distribution. +def redshifts_from_nz( + count: int | ArrayLike, + z: ArrayLike, + nz: ArrayLike, + *, + rng: np.random.Generator | None = None, +) -> np.ndarray: + """Generate galaxy redshifts from a source distribution. The function supports sampling from multiple populations of redshifts if *count* is an array or if there are additional axes in @@ -87,7 +91,7 @@ def redshifts_from_nz(count: int | ArrayLike, z: ArrayLike, nz: ArrayLike, *, inputs with extra dimensions, returns a flattened 1-D array of samples from all populations. - ''' + """ # get default RNG if not given if rng is None: @@ -104,13 +108,14 @@ def redshifts_from_nz(count: int | ArrayLike, z: ArrayLike, nz: ArrayLike, *, # go through extra dimensions; also works if dims is empty for k in np.ndindex(dims): - # compute the CDF of each galaxy population cdf = cumtrapz(nz[k], z[k], dtype=float) cdf /= cdf[-1] # sample redshifts and store result - redshifts[total:total+count[k]] = np.interp(rng.uniform(0, 1, size=count[k]), cdf, z[k]) + redshifts[total : total + count[k]] = np.interp( + rng.uniform(0, 1, size=count[k]), cdf, z[k] + ) total += count[k] assert total == redshifts.size @@ -118,10 +123,17 @@ def redshifts_from_nz(count: int | ArrayLike, z: ArrayLike, nz: ArrayLike, *, return redshifts -def galaxy_shear(lon: np.ndarray, lat: np.ndarray, eps: np.ndarray, - kappa: np.ndarray, gamma1: np.ndarray, gamma2: np.ndarray, *, - reduced_shear: bool = True) -> np.ndarray: - '''Observed galaxy shears from weak lensing. +def galaxy_shear( + lon: np.ndarray, + lat: np.ndarray, + eps: np.ndarray, + kappa: np.ndarray, + gamma1: np.ndarray, + gamma2: np.ndarray, + *, + reduced_shear: bool = True, +) -> np.ndarray: + """Observed galaxy shears from weak lensing. Takes lensing maps for convergence and shear and produces a lensed ellipticity (shear) for each intrinsic galaxy ellipticity. @@ -143,7 +155,7 @@ def galaxy_shear(lon: np.ndarray, lat: np.ndarray, eps: np.ndarray, she : array_like Array of complex-valued observed galaxy shears (lensed ellipticities). - ''' + """ nside = healpix.npix2nside(np.broadcast(kappa, gamma1, gamma2).shape[-1]) @@ -155,7 +167,7 @@ def galaxy_shear(lon: np.ndarray, lat: np.ndarray, eps: np.ndarray, # get the lensing maps at galaxy position for i in range(0, size, 10000): - s = slice(i, i+10000) + s = slice(i, i + 10000) ipix = healpix.ang2pix(nside, lon[s], lat[s], lonlat=True) k[s] = kappa[ipix] g.real[s] = gamma1[ipix] @@ -166,7 +178,7 @@ def galaxy_shear(lon: np.ndarray, lat: np.ndarray, eps: np.ndarray, g /= 1 - k # compute lensed ellipticities - g = (eps + g)/(1 + g.conj()*eps) + g = (eps + g) / (1 + g.conj() * eps) else: # simple sum of shears g += eps @@ -174,11 +186,15 @@ def galaxy_shear(lon: np.ndarray, lat: np.ndarray, eps: np.ndarray, return g -def gaussian_phz(z: ArrayLike, sigma_0: float | ArrayLike, *, - lower: ArrayLike | None = None, - upper: ArrayLike | None = None, - rng: np.random.Generator | None = None) -> np.ndarray: - r'''Photometric redshifts assuming a Gaussian error. +def gaussian_phz( + z: ArrayLike, + sigma_0: float | ArrayLike, + *, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + rng: np.random.Generator | None = None, +) -> np.ndarray: + r"""Photometric redshifts assuming a Gaussian error. A simple toy model of photometric redshift errors that assumes a Gaussian error with redshift-dependent standard deviation @@ -222,19 +238,19 @@ def gaussian_phz(z: ArrayLike, sigma_0: float | ArrayLike, *, -------- See the :doc:`/examples/1-basic/photoz` example. - ''' + """ # get default RNG if not given if rng is None: rng = np.random.default_rng() - sigma = np.add(1, z)*sigma_0 + sigma = np.add(1, z) * sigma_0 dims = np.shape(sigma) zphot = rng.normal(z, sigma) if lower is None: - lower = 0. + lower = 0.0 if upper is None: upper = np.inf diff --git a/glass/lensing.py b/glass/lensing.py index ee01e54c..b1ab7284 100644 --- a/glass/lensing.py +++ b/glass/lensing.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Lensing (:mod:`glass.lensing`) ============================== @@ -29,7 +29,7 @@ .. autofunction:: deflect -''' +""" import numpy as np import healpy as hp @@ -37,19 +37,23 @@ # typing support from typing import Optional, Sequence, Tuple, TYPE_CHECKING from numpy.typing import NDArray, ArrayLike + if TYPE_CHECKING: # to prevent circular dependencies, only import these for type checking from cosmology import Cosmology from .shells import RadialWindow -def from_convergence(kappa: NDArray, lmax: Optional[int] = None, *, - potential: bool = False, - deflection: bool = False, - shear: bool = False, - discretized: bool = True, - ) -> Tuple[NDArray, ...]: - r'''Compute other weak lensing maps from the convergence. +def from_convergence( + kappa: NDArray, + lmax: Optional[int] = None, + *, + potential: bool = False, + deflection: bool = False, + shear: bool = False, + discretized: bool = True, +) -> Tuple[NDArray, ...]: + r"""Compute other weak lensing maps from the convergence. Takes a weak lensing convergence map and returns one or more of deflection potential, deflection, and shear maps. The maps are @@ -147,7 +151,7 @@ def from_convergence(kappa: NDArray, lmax: Optional[int] = None, *, .. [1] Tessore N., et al., OJAp, 6, 11 (2023). doi:10.21105/astro.2302.01942 - ''' + """ # no output means no computation, return empty tuple if not (potential or deflection or shear): @@ -156,19 +160,19 @@ def from_convergence(kappa: NDArray, lmax: Optional[int] = None, *, # get the NSIDE parameter nside = hp.get_nside(kappa) if lmax is None: - lmax = 3*nside - 1 + lmax = 3 * nside - 1 # compute alm alm = hp.map2alm(kappa, lmax=lmax, pol=False, use_pixel_weights=True) # mode number; all conversions are factors of this - l = np.arange(lmax+1) + l = np.arange(lmax + 1) # this tuple will be returned results = () # convert convergence to potential - fl = np.divide(-2, l*(l+1), where=(l > 0), out=np.zeros(lmax+1)) + fl = np.divide(-2, l * (l + 1), where=(l > 0), out=np.zeros(lmax + 1)) hp.almxfl(alm, fl, inplace=True) # if potential is requested, compute map and add to output @@ -184,14 +188,14 @@ def from_convergence(kappa: NDArray, lmax: Optional[int] = None, *, blm = np.zeros_like(alm) # compute deflection alms in place - fl = np.sqrt(l*(l+1)) + fl = np.sqrt(l * (l + 1)) # TODO: missing spin-1 pixel window function here hp.almxfl(alm, fl, inplace=True) # if deflection is requested, compute spin-1 maps and add to output if deflection: alpha = hp.alm2map_spin([alm, blm], nside, 1, lmax) - alpha = alpha[0] + 1j*alpha[1] + alpha = alpha[0] + 1j * alpha[1] results += (alpha,) # if no shear is requested, stop here @@ -200,25 +204,26 @@ def from_convergence(kappa: NDArray, lmax: Optional[int] = None, *, # compute shear alms in place # if discretised, factor out spin-0 kernel and apply spin-2 kernel - fl = np.sqrt((l-1)*(l+2), where=(l > 0), out=np.zeros(lmax+1)) + fl = np.sqrt((l - 1) * (l + 2), where=(l > 0), out=np.zeros(lmax + 1)) fl /= 2 if discretized: pw0, pw2 = hp.pixwin(nside, lmax=lmax, pol=True) - fl *= pw2/pw0 + fl *= pw2 / pw0 hp.almxfl(alm, fl, inplace=True) # transform to shear maps gamma = hp.alm2map_spin([alm, blm], nside, 2, lmax) - gamma = gamma[0] + 1j*gamma[1] + gamma = gamma[0] + 1j * gamma[1] results += (gamma,) # all done return results -def shear_from_convergence(kappa: np.ndarray, lmax: Optional[int] = None, *, - discretized: bool = True) -> np.ndarray: - r'''Weak lensing shear from convergence. +def shear_from_convergence( + kappa: np.ndarray, lmax: Optional[int] = None, *, discretized: bool = True +) -> np.ndarray: + r"""Weak lensing shear from convergence. .. deprecated:: 2023.6 Use the more general :func:`from_convergence` function instead. @@ -226,11 +231,11 @@ def shear_from_convergence(kappa: np.ndarray, lmax: Optional[int] = None, *, Computes the shear from the convergence using a spherical harmonic transform. - ''' + """ nside = hp.get_nside(kappa) if lmax is None: - lmax = 3*nside - 1 + lmax = 3 * nside - 1 # compute alm alm = hp.map2alm(kappa, lmax=lmax, pol=False, use_pixel_weights=True) @@ -239,15 +244,15 @@ def shear_from_convergence(kappa: np.ndarray, lmax: Optional[int] = None, *, blm = np.zeros_like(alm) # factor to convert convergence alm to shear alm - l = np.arange(lmax+1) - fl = np.sqrt((l+2)*(l+1)*l*(l-1)) - fl /= np.clip(l*(l+1), 1, None) + l = np.arange(lmax + 1) + fl = np.sqrt((l + 2) * (l + 1) * l * (l - 1)) + fl /= np.clip(l * (l + 1), 1, None) fl *= -1 # if discretised, factor out spin-0 kernel and apply spin-2 kernel if discretized: pw0, pw2 = hp.pixwin(nside, lmax=lmax, pol=True) - fl *= pw2/pw0 + fl *= pw2 / pw0 # apply correction to E-modes hp.almxfl(alm, fl, inplace=True) @@ -257,41 +262,40 @@ def shear_from_convergence(kappa: np.ndarray, lmax: Optional[int] = None, *, class MultiPlaneConvergence: - '''Compute convergence fields iteratively from multiple matter planes.''' + """Compute convergence fields iteratively from multiple matter planes.""" - def __init__(self, cosmo: 'Cosmology') -> None: - '''Create a new instance to iteratively compute the convergence.''' + def __init__(self, cosmo: "Cosmology") -> None: + """Create a new instance to iteratively compute the convergence.""" self.cosmo = cosmo # set up initial values of variables - self.z2: float = 0. - self.z3: float = 0. - self.x3: float = 0. - self.w3: float = 0. - self.r23: float = 1. - self.delta3: np.ndarray = np.array(0.) + self.z2: float = 0.0 + self.z3: float = 0.0 + self.x3: float = 0.0 + self.w3: float = 0.0 + self.r23: float = 1.0 + self.delta3: np.ndarray = np.array(0.0) self.kappa2: Optional[np.ndarray] = None self.kappa3: Optional[np.ndarray] = None - def add_window(self, delta: np.ndarray, w: 'RadialWindow') -> None: - '''Add a mass plane from a window function to the convergence. + def add_window(self, delta: np.ndarray, w: "RadialWindow") -> None: + """Add a mass plane from a window function to the convergence. The lensing weight is computed from the window function, and the source plane redshift is the effective redshift of the window. - ''' + """ zsrc = w.zeff - lens_weight = np.trapz(w.wa, w.za)/np.interp(zsrc, w.za, w.wa) + lens_weight = np.trapz(w.wa, w.za) / np.interp(zsrc, w.za, w.wa) self.add_plane(delta, zsrc, lens_weight) - def add_plane(self, delta: np.ndarray, zsrc: float, wlens: float = 1. - ) -> None: - '''Add a mass plane at redshift ``zsrc`` to the convergence.''' + def add_plane(self, delta: np.ndarray, zsrc: float, wlens: float = 1.0) -> None: + """Add a mass plane at redshift ``zsrc`` to the convergence.""" if zsrc <= self.z3: - raise ValueError('source redshift must be increasing') + raise ValueError("source redshift must be increasing") # cycle mass plane, ... delta2, self.delta3 = self.delta3, delta @@ -305,13 +309,13 @@ def add_plane(self, delta: np.ndarray, zsrc: float, wlens: float = 1. # extrapolation law x2, self.x3 = self.x3, self.cosmo.xm(self.z3) r12 = self.r23 - r13, self.r23 = self.cosmo.xm([z1, self.z2], self.z3)/self.x3 - t = r13/r12 + r13, self.r23 = self.cosmo.xm([z1, self.z2], self.z3) / self.x3 + t = r13 / r12 # lensing weight of mass plane to be added - f = 3*self.cosmo.omega_m/2 - f *= x2*self.r23 - f *= (1 + self.z2)/self.cosmo.ef(self.z2) + f = 3 * self.cosmo.omega_m / 2 + f *= x2 * self.r23 + f *= (1 + self.z2) / self.cosmo.ef(self.z2) f *= w2 # create kappa planes on first iteration @@ -327,27 +331,27 @@ def add_plane(self, delta: np.ndarray, zsrc: float, wlens: float = 1. # compute next convergence plane in place of last self.kappa3 *= 1 - t - self.kappa3 += t*self.kappa2 - self.kappa3 += f*delta2 + self.kappa3 += t * self.kappa2 + self.kappa3 += f * delta2 @property def zsrc(self) -> float: - '''The redshift of the current convergence plane.''' + """The redshift of the current convergence plane.""" return self.z3 @property def kappa(self) -> Optional[np.ndarray]: - '''The current convergence plane.''' + """The current convergence plane.""" return self.kappa3 @property def delta(self) -> np.ndarray: - '''The current matter plane.''' + """The current matter plane.""" return self.delta3 @property def wlens(self) -> float: - '''The weight of the current matter plane.''' + """The weight of the current matter plane.""" return self.w3 @@ -448,16 +452,16 @@ def deflect(lon: ArrayLike, lat: ArrayLike, alpha: ArrayLike) -> NDArray: # δ = arctan2(sin|α| sinγ, sinθ cos|α| - cosθ sin|α| cosγ) t = np.radians(lat) - ct, st = np.sin(t), np.cos(t) # sin and cos flipped: lat not co-lat + ct, st = np.sin(t), np.cos(t) # sin and cos flipped: lat not co-lat - a = np.hypot(alpha1, alpha2) # abs(alpha) + a = np.hypot(alpha1, alpha2) # abs(alpha) g = np.arctan2(alpha2, alpha1) # arg(alpha) ca, sa = np.cos(a), np.sin(a) cg, sg = np.cos(g), np.sin(g) # flipped atan2 arguments for lat instead of co-lat - tp = np.arctan2(ct*ca + st*sa*cg, np.hypot(ct*sa - st*ca*cg, st*sg)) + tp = np.arctan2(ct * ca + st * sa * cg, np.hypot(ct * sa - st * ca * cg, st * sg)) - d = np.arctan2(sa*sg, st*ca - ct*sa*cg) + d = np.arctan2(sa * sg, st * ca - ct * sa * cg) return lon - np.degrees(d), np.degrees(tp) diff --git a/glass/observations.py b/glass/observations.py index 44883e02..acd8b5f9 100644 --- a/glass/observations.py +++ b/glass/observations.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Observations (:mod:`glass.observations`) ======================================== @@ -26,7 +26,7 @@ .. autofunction:: vmap_galactic_ecliptic -''' +""" import numpy as np import healpy as hp @@ -38,10 +38,12 @@ from .core.array import cumtrapz -def vmap_galactic_ecliptic(nside: int, galactic: Tuple[float, float] = (30, 90), - ecliptic: Tuple[float, float] = (20, 80) - ) -> np.ndarray: - '''visibility map masking galactic and ecliptic plane +def vmap_galactic_ecliptic( + nside: int, + galactic: Tuple[float, float] = (30, 90), + ecliptic: Tuple[float, float] = (20, 80), +) -> np.ndarray: + """visibility map masking galactic and ecliptic plane This function returns a :term:`visibility map` that blocks out stripes for the galactic and ecliptic planes. The location of the stripes is set with @@ -65,24 +67,29 @@ def vmap_galactic_ecliptic(nside: int, galactic: Tuple[float, float] = (30, 90), TypeError If the ``galactic`` or ``ecliptic`` arguments are not pairs of numbers. - ''' + """ if np.ndim(galactic) != 1 or len(galactic) != 2: - raise TypeError('galactic stripe must be a pair of numbers') + raise TypeError("galactic stripe must be a pair of numbers") if np.ndim(ecliptic) != 1 or len(ecliptic) != 2: - raise TypeError('ecliptic stripe must be a pair of numbers') + raise TypeError("ecliptic stripe must be a pair of numbers") m = np.ones(hp.nside2npix(nside)) m[hp.query_strip(nside, *galactic)] = 0 - m = hp.Rotator(coord='GC').rotate_map_pixel(m) + m = hp.Rotator(coord="GC").rotate_map_pixel(m) m[hp.query_strip(nside, *ecliptic)] = 0 - m = hp.Rotator(coord='CE').rotate_map_pixel(m) + m = hp.Rotator(coord="CE").rotate_map_pixel(m) return m -def gaussian_nz(z: np.ndarray, mean: ArrayLike, sigma: ArrayLike, *, - norm: Optional[ArrayLike] = None) -> np.ndarray: - r'''Gaussian redshift distribution. +def gaussian_nz( + z: np.ndarray, + mean: ArrayLike, + sigma: ArrayLike, + *, + norm: Optional[ArrayLike] = None, +) -> np.ndarray: + r"""Gaussian redshift distribution. The redshift follows a Gaussian distribution with the given mean and standard deviation. @@ -106,11 +113,11 @@ def gaussian_nz(z: np.ndarray, mean: ArrayLike, sigma: ArrayLike, *, nz : array_like Redshift distribution at the given ``z`` values. - ''' - mean = np.reshape(mean, np.shape(mean) + (1,)*np.ndim(z)) - sigma = np.reshape(sigma, np.shape(sigma) + (1,)*np.ndim(z)) + """ + mean = np.reshape(mean, np.shape(mean) + (1,) * np.ndim(z)) + sigma = np.reshape(sigma, np.shape(sigma) + (1,) * np.ndim(z)) - nz = np.exp(-((z - mean)/sigma)**2/2) + nz = np.exp(-(((z - mean) / sigma) ** 2) / 2) nz /= np.trapz(nz, z, axis=-1)[..., np.newaxis] if norm is not None: @@ -119,10 +126,15 @@ def gaussian_nz(z: np.ndarray, mean: ArrayLike, sigma: ArrayLike, *, return nz -def smail_nz(z: np.ndarray, z_mode: ArrayLike, alpha: ArrayLike, - beta: ArrayLike, *, norm: Optional[ArrayLike] = None - ) -> np.ndarray: - r'''Redshift distribution following Smail et al. (1994). +def smail_nz( + z: np.ndarray, + z_mode: ArrayLike, + alpha: ArrayLike, + beta: ArrayLike, + *, + norm: Optional[ArrayLike] = None, +) -> np.ndarray: + r"""Redshift distribution following Smail et al. (1994). The redshift follows the Smail et al. [1]_ redshift distribution. @@ -161,12 +173,12 @@ def smail_nz(z: np.ndarray, z_mode: ArrayLike, alpha: ArrayLike, .. [1] Smail I., Ellis R. S., Fitchett M. J., 1994, MNRAS, 270, 245 .. [2] Amara A., Refregier A., 2007, MNRAS, 381, 1018 - ''' + """ z_mode = np.asanyarray(z_mode)[..., np.newaxis] alpha = np.asanyarray(alpha)[..., np.newaxis] beta = np.asanyarray(beta)[..., np.newaxis] - pz = z**alpha*np.exp(-alpha/beta*(z/z_mode)**beta) + pz = z**alpha * np.exp(-alpha / beta * (z / z_mode) ** beta) pz /= np.trapz(pz, z, axis=-1)[..., np.newaxis] if norm is not None: @@ -175,10 +187,10 @@ def smail_nz(z: np.ndarray, z_mode: ArrayLike, alpha: ArrayLike, return pz -def fixed_zbins(zmin: float, zmax: float, *, - nbins: Optional[int] = None, dz: Optional[float] = None - ) -> List[Tuple[float, float]]: - '''tomographic redshift bins of fixed size +def fixed_zbins( + zmin: float, zmax: float, *, nbins: Optional[int] = None, dz: Optional[float] = None +) -> List[Tuple[float, float]]: + """tomographic redshift bins of fixed size This function creates contiguous tomographic redshift bins of fixed size. It takes either the number or size of the bins. @@ -196,21 +208,22 @@ def fixed_zbins(zmin: float, zmax: float, *, ------- zbins : list of tuple of float List of redshift bin edges. - ''' + """ if nbins is not None and dz is None: - zbinedges = np.linspace(zmin, zmax, nbins+1) + zbinedges = np.linspace(zmin, zmax, nbins + 1) if nbins is None and dz is not None: zbinedges = np.arange(zmin, zmax, dz) else: - raise ValueError('exactly one of nbins and dz must be given') + raise ValueError("exactly one of nbins and dz must be given") return list(zip(zbinedges, zbinedges[1:])) -def equal_dens_zbins(z: np.ndarray, nz: np.ndarray, nbins: int - ) -> List[Tuple[float, float]]: - '''equal density tomographic redshift bins +def equal_dens_zbins( + z: np.ndarray, nz: np.ndarray, nbins: int +) -> List[Tuple[float, float]]: + """equal density tomographic redshift bins This function subdivides a source redshift distribution into ``nbins`` tomographic redshift bins with equal density. @@ -227,21 +240,22 @@ def equal_dens_zbins(z: np.ndarray, nz: np.ndarray, nbins: int zbins : list of tuple of float List of redshift bin edges. - ''' + """ # compute the normalised cumulative distribution function # 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 = cumtrapz(nz, z) cuml_nz /= cuml_nz[[-1]] - zbinedges = np.interp(np.linspace(0, 1, nbins+1), cuml_nz, z) + zbinedges = np.interp(np.linspace(0, 1, nbins + 1), cuml_nz, z) return list(zip(zbinedges, zbinedges[1:])) -def tomo_nz_gausserr(z: np.ndarray, nz: np.ndarray, sigma_0: float, - zbins: List[Tuple[float, float]]) -> np.ndarray: - '''tomographic redshift bins with a Gaussian redshift error +def tomo_nz_gausserr( + z: np.ndarray, nz: np.ndarray, sigma_0: float, zbins: List[Tuple[float, float]] +) -> np.ndarray: + """tomographic redshift bins with a Gaussian redshift error This function takes a _true_ overall source redshift distribution ``z``, ``nz`` and returns tomographic source redshift distributions for the @@ -277,7 +291,7 @@ def tomo_nz_gausserr(z: np.ndarray, nz: np.ndarray, sigma_0: float, .. [1] Amara A., Réfrégier A., 2007, MNRAS, 381, 1018. doi:10.1111/j.1365-2966.2007.12271.x - ''' + """ # converting zbins into an array: zbins_arr = np.asanyarray(zbins) # type: ignore[no-redef] @@ -291,10 +305,10 @@ def tomo_nz_gausserr(z: np.ndarray, nz: np.ndarray, sigma_0: float, # compute the probabilities that redshifts z end up in each bin # then apply probability as weights to given nz # leading axis corresponds to the different bins - sz = 2**0.5*sigma_0*(1 + z) - binned_nz = erf((z - z_lower)/sz) - binned_nz -= erf((z - z_upper)/sz) - binned_nz /= 1 + erf(z/sz) + sz = 2**0.5 * sigma_0 * (1 + z) + binned_nz = erf((z - z_lower) / sz) + binned_nz -= erf((z - z_upper) / sz) + binned_nz /= 1 + erf(z / sz) binned_nz *= nz return binned_nz diff --git a/glass/points.py b/glass/points.py index c2b4ddae..0c35ebb9 100644 --- a/glass/points.py +++ b/glass/points.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Random points (:mod:`glass.points`) =================================== @@ -29,7 +29,7 @@ .. autofunction:: linear_bias .. autofunction:: loglinear_bias -''' +""" import numpy as np import healpix @@ -39,7 +39,7 @@ def effective_bias(z, bz, w): - '''Effective bias parameter from a redshift-dependent bias function. + """Effective bias parameter from a redshift-dependent bias function. This function takes a redshift-dependent bias function :math:`b(z)` and computes an effective bias parameter :math:`\\bar{b}` for a @@ -67,28 +67,36 @@ def effective_bias(z, bz, w): \\bar{b} = \\frac{\\int b(z) \\, w(z) \\, dz}{\\int w(z) \\, dz} \\;. - ''' + """ norm = np.trapz(w.wa, w.za) - return trapz_product((z, bz), (w.za, w.wa))/norm + return trapz_product((z, bz), (w.za, w.wa)) / norm def linear_bias(delta, b): - '''linear bias model :math:`\\delta_g = b \\, \\delta`''' - return b*delta + """linear bias model :math:`\\delta_g = b \\, \\delta`""" + return b * delta def loglinear_bias(delta, b): - '''log-linear bias model :math:`\\ln(1 + \\delta_g) = b \\ln(1 + \\delta)`''' + """log-linear bias model :math:`\\ln(1 + \\delta_g) = b \\ln(1 + \\delta)`""" delta_g = np.log1p(delta) delta_g *= b np.expm1(delta_g, out=delta_g) return delta_g -def positions_from_delta(ngal, delta, bias=None, vis=None, *, - bias_model='linear', remove_monopole=False, - batch=1_000_000, rng=None): - '''Generate positions tracing a density contrast. +def positions_from_delta( + ngal, + delta, + bias=None, + vis=None, + *, + bias_model="linear", + remove_monopole=False, + batch=1_000_000, + rng=None, +): + """Generate positions tracing a density contrast. The map of expected number counts is constructed from the number density, density contrast, an optional bias model, and an optional @@ -143,7 +151,7 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, sampled, an array of counts in the shape of the extra dimensions is returned. - ''' + """ # get default RNG if not given if rng is None: @@ -151,9 +159,9 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, # get the bias model if isinstance(bias_model, str): - bias_model = globals()[f'{bias_model}_bias'] + bias_model = globals()[f"{bias_model}_bias"] elif not callable(bias_model): - raise ValueError('bias_model must be string or callable') + raise ValueError("bias_model must be string or callable") # broadcast inputs to common shape of extra dimensions inputs = [(ngal, 0), (delta, 1)] @@ -169,7 +177,6 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, # iterate the leading dimensions for k in np.ndindex(dims): - # compute density contrast from bias model, or copy if bias is None: n = np.copy(delta[k]) @@ -182,7 +189,7 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, # turn into number count, modifying the array in place n += 1 - n *= ARCMIN2_SPHERE/n.size*ngal[k] + n *= ARCMIN2_SPHERE / n.size * ngal[k] # apply visibility if given if vis is not None: @@ -216,7 +223,7 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, start, stop, size = 0, 0, 0 while count: # tally this group of pixels - q = np.cumsum(n[stop:stop+step]) + q = np.cumsum(n[stop : stop + step]) # does this group of pixels fill the batch? if size + q[-1] < min(batch, count): # no, we need the next group of pixels to fill the batch @@ -224,7 +231,7 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, size += q[-1] else: # how many pixels from this group do we need? - stop += np.searchsorted(q, batch - size, side='right') + stop += np.searchsorted(q, batch - size, side="right") # if the first pixel alone is too much, use it anyway if stop == start: stop += 1 @@ -236,14 +243,14 @@ def positions_from_delta(ngal, delta, bias=None, vis=None, *, # keep track of remaining number of points count -= ipix.size # yield the batch - yield lon, lat, ipix.size*cmask + yield lon, lat, ipix.size * cmask # make sure that the correct number of pixels was sampled assert np.sum(n[stop:]) == 0 def uniform_positions(ngal, *, rng=None): - '''Generate positions uniformly over the sphere. + """Generate positions uniformly over the sphere. The function supports array input for the ``ngal`` parameter. @@ -262,7 +269,7 @@ def uniform_positions(ngal, *, rng=None): The number of sampled points. For array inputs, an array of counts with the same shape is returned. - ''' + """ # get default RNG if not given if rng is None: @@ -279,7 +286,6 @@ def uniform_positions(ngal, *, rng=None): # sample each set of points for k in np.ndindex(dims): - # sample uniformly over the sphere lon = rng.uniform(-180, 180, size=ngal[k]) lat = np.rad2deg(np.arcsin(rng.uniform(-1, 1, size=ngal[k]))) diff --git a/glass/shapes.py b/glass/shapes.py index 2d2d9acd..3a71f3a5 100644 --- a/glass/shapes.py +++ b/glass/shapes.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Observed shapes (:mod:`glass.shapes`) ===================================== @@ -22,7 +22,7 @@ .. autofunction:: triaxial_axis_ratio -''' +""" from __future__ import annotations @@ -31,7 +31,7 @@ def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None): - r'''axis ratio of a randomly projected triaxial ellipsoid + r"""axis ratio of a randomly projected triaxial ellipsoid Given the two axis ratios `1 >= zeta >= xi` of a randomly oriented triaxial ellipsoid, computes the axis ratio `q` of the projection. @@ -61,7 +61,7 @@ def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None): ---------- .. [1] Binney J., 1985, MNRAS, 212, 767. doi:10.1093/mnras/212.4.767 - ''' + """ # default RNG if not provided if rng is None: @@ -72,10 +72,10 @@ def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None): size = np.broadcast(zeta, xi).shape # draw random viewing angle (theta, phi) - cos2_theta = rng.uniform(low=-1., high=1., size=size) + cos2_theta = rng.uniform(low=-1.0, high=1.0, size=size) cos2_theta *= cos2_theta sin2_theta = 1 - cos2_theta - cos2_phi = np.cos(rng.uniform(low=0., high=2*np.pi, size=size)) + cos2_phi = np.cos(rng.uniform(low=0.0, high=2 * np.pi, size=size)) cos2_phi *= cos2_phi sin2_phi = 1 - cos2_phi @@ -85,18 +85,20 @@ def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None): x2 = np.square(xi) # eq. (11) multiplied by xi^2 zeta^2 - A = (1 + z2m1*sin2_phi)*cos2_theta + x2*sin2_theta - B2 = 4*z2m1**2*cos2_theta*sin2_phi*cos2_phi - C = 1 + z2m1*cos2_phi + A = (1 + z2m1 * sin2_phi) * cos2_theta + x2 * sin2_theta + B2 = 4 * z2m1**2 * cos2_theta * sin2_phi * cos2_phi + C = 1 + z2m1 * cos2_phi # eq. (12) - q = np.sqrt((A+C-np.sqrt((A-C)**2+B2))/(A+C+np.sqrt((A-C)**2+B2))) + q = np.sqrt( + (A + C - np.sqrt((A - C) ** 2 + B2)) / (A + C + np.sqrt((A - C) ** 2 + B2)) + ) return q def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None): - r'''ellipticity distribution following Ryden (2004) + r"""ellipticity distribution following Ryden (2004) The ellipticities are sampled by randomly projecting a 3D ellipsoid with principal axes :math:`A > B > C` [1]_. The distribution of :math:`\log(1 - @@ -130,7 +132,7 @@ def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None): .. [1] Ryden B. S., 2004, ApJ, 601, 214. .. [2] Padilla N. D., Strauss M. A., 2008, MNRAS, 388, 1321. - ''' + """ # default RNG if not provided if rng is None: @@ -139,10 +141,10 @@ def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None): # draw gamma and epsilon from truncated normal -- eq.s (10)-(11) # first sample unbounded normal, then rejection sample truncation eps = rng.normal(mu, sigma, size=size) - bad = (eps > 0) + bad = eps > 0 while np.any(bad): eps[bad] = rng.normal(mu, sigma, size=eps[bad].shape) - bad = (eps > 0) + bad = eps > 0 gam = rng.normal(gamma, sigma_gamma, size=size) bad = (gam < 0) | (gam > 1) while np.any(bad): @@ -151,23 +153,23 @@ def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None): # compute triaxial axis ratios zeta = B/A, xi = C/A zeta = -np.expm1(eps) - xi = (1 - gam)*zeta + xi = (1 - gam) * zeta # random projection of random triaxial ellipsoid q = triaxial_axis_ratio(zeta, xi, rng=rng) # assemble ellipticity with random complex phase - e = np.exp(1j*rng.uniform(0, 2*np.pi, size=np.shape(q))) - e *= (1-q)/(1+q) + e = np.exp(1j * rng.uniform(0, 2 * np.pi, size=np.shape(q))) + e *= (1 - q) / (1 + q) # return the ellipticity return e -def ellipticity_gaussian(count: int | ArrayLike, sigma: ArrayLike, *, - rng: np.random.Generator | None = None - ) -> NDArray: - r'''Sample Gaussian galaxy ellipticities. +def ellipticity_gaussian( + count: int | ArrayLike, sigma: ArrayLike, *, rng: np.random.Generator | None = None +) -> NDArray: + r"""Sample Gaussian galaxy ellipticities. The ellipticities are sampled from a normal distribution with standard deviation ``sigma`` for each component. Samples where the @@ -189,7 +191,7 @@ def ellipticity_gaussian(count: int | ArrayLike, sigma: ArrayLike, *, eps : array_like Array of galaxy :term:`ellipticity`. - ''' + """ # default RNG if not provided if rng is None: @@ -205,23 +207,23 @@ def ellipticity_gaussian(count: int | ArrayLike, sigma: ArrayLike, *, # reject those where abs(e) > 0 i = 0 for k in np.ndindex(count.shape): - e = rng.standard_normal(2*count[k], np.float64).view(np.complex128) + e = rng.standard_normal(2 * count[k], np.float64).view(np.complex128) e *= sigma[k] r = np.where(np.abs(e) > 1)[0] while len(r) > 0: - rng.standard_normal(2*len(r), np.float64, e[r].view(np.float64)) + rng.standard_normal(2 * len(r), np.float64, e[r].view(np.float64)) e[r] *= sigma[k] r = r[np.abs(e[r]) > 1] - eps[i:i+count[k]] = e + eps[i : i + count[k]] = e i += count[k] return eps -def ellipticity_intnorm(count: int | ArrayLike, sigma: ArrayLike, *, - rng: np.random.Generator | None = None - ) -> NDArray: - r'''Sample galaxy ellipticities with intrinsic normal distribution. +def ellipticity_intnorm( + count: int | ArrayLike, sigma: ArrayLike, *, rng: np.random.Generator | None = None +) -> NDArray: + r"""Sample galaxy ellipticities with intrinsic normal distribution. The ellipticities are sampled from an intrinsic normal distribution with standard deviation ``sigma`` for each component. @@ -240,7 +242,7 @@ def ellipticity_intnorm(count: int | ArrayLike, sigma: ArrayLike, *, eps : array_like Array of galaxy :term:`ellipticity`. - ''' + """ # default RNG if not provided if rng is None: @@ -251,10 +253,10 @@ def ellipticity_intnorm(count: int | ArrayLike, sigma: ArrayLike, *, # make sure sigma is admissible if not np.all((0 <= sigma) & (sigma < 0.5**0.5)): - raise ValueError('sigma must be between 0 and sqrt(0.5)') + raise ValueError("sigma must be between 0 and sqrt(0.5)") # convert to sigma_eta using fit - sigma_eta = sigma*((8 + 5*sigma**2)/(2 - 4*sigma**2))**0.5 + sigma_eta = sigma * ((8 + 5 * sigma**2) / (2 - 4 * sigma**2)) ** 0.5 # allocate flattened output array eps = np.empty(count.sum(), dtype=np.complex128) @@ -262,11 +264,11 @@ def ellipticity_intnorm(count: int | ArrayLike, sigma: ArrayLike, *, # sample complex ellipticities i = 0 for k in np.ndindex(count.shape): - e = rng.standard_normal(2*count[k], np.float64).view(np.complex128) + e = rng.standard_normal(2 * count[k], np.float64).view(np.complex128) e *= sigma_eta[k] r = np.hypot(e.real, e.imag) - e *= np.divide(np.tanh(r/2), r, where=(r > 0), out=r) - eps[i:i+count[k]] = e + e *= np.divide(np.tanh(r / 2), r, where=(r > 0), out=r) + eps[i : i + count[k]] = e i += count[k] return eps diff --git a/glass/shells.py b/glass/shells.py index b1667d6d..000e9d9e 100644 --- a/glass/shells.py +++ b/glass/shells.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" Shells (:mod:`glass.shells`) ============================ @@ -43,7 +43,7 @@ .. autofunction:: volume_weight .. autofunction:: density_weight -''' +""" import warnings from collections import namedtuple @@ -52,9 +52,9 @@ from .core.array import ndinterp # type checking -from typing import (Union, Sequence, List, Tuple, Optional, Callable, - TYPE_CHECKING) +from typing import Union, Sequence, List, Tuple, Optional, Callable, TYPE_CHECKING from numpy.typing import ArrayLike + if TYPE_CHECKING: from cosmology import Cosmology @@ -63,23 +63,23 @@ WeightFunc = Callable[[ArrayLike1D], np.ndarray] -def distance_weight(z: ArrayLike, cosmo: 'Cosmology') -> np.ndarray: - '''Uniform weight in comoving distance.''' - return 1/cosmo.ef(z) +def distance_weight(z: ArrayLike, cosmo: "Cosmology") -> np.ndarray: + """Uniform weight in comoving distance.""" + return 1 / cosmo.ef(z) -def volume_weight(z: ArrayLike, cosmo: 'Cosmology') -> np.ndarray: - '''Uniform weight in comoving volume.''' - return cosmo.xm(z)**2/cosmo.ef(z) +def volume_weight(z: ArrayLike, cosmo: "Cosmology") -> np.ndarray: + """Uniform weight in comoving volume.""" + return cosmo.xm(z) ** 2 / cosmo.ef(z) -def density_weight(z: ArrayLike, cosmo: 'Cosmology') -> np.ndarray: - '''Uniform weight in matter density.''' - return cosmo.rho_m_z(z)*cosmo.xm(z)**2/cosmo.ef(z) +def density_weight(z: ArrayLike, cosmo: "Cosmology") -> np.ndarray: + """Uniform weight in matter density.""" + return cosmo.rho_m_z(z) * cosmo.xm(z) ** 2 / cosmo.ef(z) -RadialWindow = namedtuple('RadialWindow', 'za, wa, zeff') -RadialWindow.__doc__ = '''A radial window, defined by a window function. +RadialWindow = namedtuple("RadialWindow", "za, wa, zeff") +RadialWindow.__doc__ = """A radial window, defined by a window function. The radial window is defined by a window function in redshift, which is given by a pair of arrays ``za``, ``wa``. @@ -119,16 +119,18 @@ def density_weight(z: ArrayLike, cosmo: 'Cosmology') -> np.ndarray: ------- _replace - ''' -RadialWindow.za.__doc__ = '''Redshift array; the abscissae of the window function.''' -RadialWindow.wa.__doc__ = '''Weight array; the values (ordinates) of the window function.''' -RadialWindow.zeff.__doc__ = '''Effective redshift of the window.''' + """ +RadialWindow.za.__doc__ = """Redshift array; the abscissae of the window function.""" +RadialWindow.wa.__doc__ = ( + """Weight array; the values (ordinates) of the window function.""" +) +RadialWindow.zeff.__doc__ = """Effective redshift of the window.""" -def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3, - weight: Optional[WeightFunc] = None - ) -> List[RadialWindow]: - '''Tophat window functions from the given redshift bin edges. +def tophat_windows( + zbins: ArrayLike1D, dz: float = 1e-3, weight: Optional[WeightFunc] = None +) -> List[RadialWindow]: + """Tophat window functions from the given redshift bin edges. Uses the *N+1* given redshifts as bin edges to construct *N* tophat window functions. The redshifts of the windows have linear spacing @@ -160,11 +162,11 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3, -------- :ref:`user-window-functions` - ''' + """ if len(zbins) < 2: - raise ValueError('zbins must have at least two entries') + raise ValueError("zbins must have at least two entries") if zbins[0] != 0: - warnings.warn('first tophat window does not start at redshift zero') + warnings.warn("first tophat window does not start at redshift zero") wht: WeightFunc if weight is not None: @@ -174,18 +176,20 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3, ws = [] for zmin, zmax in zip(zbins, zbins[1:]): - n = max(round((zmax - zmin)/dz), 2) + n = max(round((zmax - zmin) / dz), 2) z = np.linspace(zmin, zmax, n) w = wht(z) - zeff = np.trapz(w*z, z)/np.trapz(w, z) + zeff = np.trapz(w * z, z) / np.trapz(w, z) ws.append(RadialWindow(z, w, zeff)) return ws -def linear_windows(zgrid: ArrayLike1D, dz: float = 1e-3, - weight: Optional[WeightFunc] = None, - ) -> List[RadialWindow]: - '''Linear interpolation window functions. +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 @@ -217,30 +221,34 @@ def linear_windows(zgrid: ArrayLike1D, dz: float = 1e-3, -------- :ref:`user-window-functions` - ''' + """ if len(zgrid) < 3: - raise ValueError('nodes must have at least 3 entries') + raise ValueError("nodes must have at least 3 entries") if zgrid[0] != 0: - warnings.warn('first triangular window does not start at z=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)]) + 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.0, 1.0, n, endpoint=False), np.linspace(1.0, 0.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. +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. @@ -273,30 +281,32 @@ def cubic_windows(zgrid: ArrayLike1D, dz: float = 1e-3, -------- :ref:`user-window-functions` - ''' + """ if len(zgrid) < 3: - raise ValueError('nodes must have at least 3 entries') + 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') + 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)]) + 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.0, 1.0, n, endpoint=False) + v = np.linspace(1.0, 0.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. +def restrict( + z: ArrayLike1D, f: ArrayLike1D, w: RadialWindow +) -> Tuple[np.ndarray, np.ndarray]: + """Restrict a function to a redshift window. Multiply the function :math:`f(z)` by a window function :math:`w(z)` to produce :math:`w(z) f(z)` over the support of :math:`w`. @@ -324,20 +334,21 @@ def restrict(z: ArrayLike1D, f: ArrayLike1D, w: RadialWindow zr, fr : array The restricted function. - ''' + """ z_ = np.compress(np.greater(z, w.za[0]) & np.less(z, w.za[-1]), z) zr = np.union1d(w.za, z_) - fr = ndinterp(zr, z, f, left=0., right=0.) * ndinterp(zr, w.za, w.wa) + fr = ndinterp(zr, z, f, left=0.0, right=0.0) * ndinterp(zr, w.za, w.wa) return zr, fr -def partition(z: ArrayLike, - fz: ArrayLike, - shells: Sequence[RadialWindow], - *, - method: str = "nnls", - ) -> ArrayLike: +def partition( + z: ArrayLike, + fz: ArrayLike, + shells: Sequence[RadialWindow], + *, + method: str = "nnls", +) -> ArrayLike: """Partition a function by a sequence of windows. Returns a vector of weights :math:`x_1, x_2, \\ldots` such that the @@ -464,16 +475,16 @@ def partition_lstsq( dz = np.gradient(zp) # create the window function matrix - a = [np.interp(zp, za, wa, left=0., right=0.) for za, wa, _ in shells] - a = a/np.trapz(a, zp, axis=-1)[..., None] - a = a*dz + a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells] + a = a / np.trapz(a, zp, axis=-1)[..., None] + a = a * dz # create the target vector of distribution values - b = ndinterp(zp, z, fz, left=0., right=0.) - b = b*dz + b = ndinterp(zp, z, fz, left=0.0, right=0.0) + b = b * dz # append a constraint for the integral - mult = 1/sumtol + mult = 1 / sumtol a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1) b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1) @@ -521,16 +532,16 @@ def partition_nnls( dz = np.gradient(zp) # create the window function matrix - a = [np.interp(zp, za, wa, left=0., right=0.) for za, wa, _ in shells] - a = a/np.trapz(a, zp, axis=-1)[..., None] - a = a*dz + a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells] + a = a / np.trapz(a, zp, axis=-1)[..., None] + a = a * dz # create the target vector of distribution values - b = ndinterp(zp, z, fz, left=0., right=0.) - b = b*dz + b = ndinterp(zp, z, fz, left=0.0, right=0.0) + b = b * dz # append a constraint for the integral - mult = 1/sumtol + mult = 1 / sumtol a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1) b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1) @@ -540,7 +551,7 @@ def partition_nnls( # reduce the dimensionality of the problem using a thin QR decomposition q, r = np.linalg.qr(a.T) - y = np.einsum('ji,...j', q, b) + y = np.einsum("ji,...j", q, b) # for each dim, find non-negative weights x such that y == r @ x x = np.empty([len(shells)] + dims) @@ -566,23 +577,23 @@ def partition_restrict( def redshift_grid(zmin, zmax, *, dz=None, num=None): - '''Redshift grid with uniform spacing in redshift.''' + """Redshift grid with uniform spacing in redshift.""" if dz is not None and num is None: - z = np.arange(zmin, np.nextafter(zmax+dz, zmax), dz) + z = np.arange(zmin, np.nextafter(zmax + dz, zmax), dz) elif dz is None and num is not None: - z = np.linspace(zmin, zmax, num+1) + z = np.linspace(zmin, zmax, num + 1) else: raise ValueError('exactly one of "dz" or "num" must be given') return z def distance_grid(cosmo, zmin, zmax, *, dx=None, num=None): - '''Redshift grid with uniform spacing in comoving distance.''' + """Redshift grid with uniform spacing in comoving distance.""" xmin, xmax = cosmo.dc(zmin), cosmo.dc(zmax) if dx is not None and num is None: - x = np.arange(xmin, np.nextafter(xmax+dx, xmax), dx) + x = np.arange(xmin, np.nextafter(xmax + dx, xmax), dx) elif dx is None and num is not None: - x = np.linspace(xmin, xmax, num+1) + x = np.linspace(xmin, xmax, num + 1) else: raise ValueError('exactly one of "dx" or "num" must be given') return cosmo.dc_inv(x) @@ -624,12 +635,13 @@ def combine( """ return sum( - np.expand_dims(weight, -1) * np.interp( + np.expand_dims(weight, -1) + * np.interp( z, shell.za, shell.wa / np.trapz(shell.wa, shell.za), - left=0., - right=0., + left=0.0, + right=0.0, ) for shell, weight in zip(shells, weights) ) diff --git a/glass/user.py b/glass/user.py index 7ff63b74..138a746e 100644 --- a/glass/user.py +++ b/glass/user.py @@ -1,6 +1,6 @@ # author: Nicolas Tessore # license: MIT -''' +""" User utilities (:mod:`glass.user`) ================================== @@ -17,19 +17,19 @@ .. autofunction:: load_cls .. autofunction:: write_catalog -''' +""" import numpy as np from contextlib import contextmanager def save_cls(filename, cls): - '''Save a list of Cls to file. + """Save a list of Cls to file. Uses :func:`numpy.savez` internally. The filename should therefore have a ``.npz`` suffix, or it will be given one. - ''' + """ split = np.cumsum([len(cl) if cl is not None else 0 for cl in cls[:-1]]) values = np.concatenate([cl for cl in cls if cl is not None]) @@ -37,28 +37,28 @@ def save_cls(filename, cls): def load_cls(filename): - '''Load a list of Cls from file. + """Load a list of Cls from file. Uses :func:`numpy.load` internally. - ''' + """ with np.load(filename) as npz: - values = npz['values'] - split = npz['split'] + values = npz["values"] + split = npz["split"] return np.split(values, split) class _FitsWriter: - '''Writer that creates a FITS file. Initialised with the fits object and extention name.''' + """Writer that creates a FITS file. Initialised with the fits object and extention name.""" def __init__(self, fits, ext=None): - '''Create a new, uninitialised writer.''' + """Create a new, uninitialised writer.""" self.fits = fits self.ext = ext def _append(self, data, names=None): - '''Internal method where the FITS writing is done''' + """Internal method where the FITS writing is done""" if self.ext is None or self.ext not in self.fits: self.fits.write_table(data, names=names, extname=self.ext) @@ -70,9 +70,9 @@ def _append(self, data, names=None): hdu.write(data, names=names, firstrow=hdu.get_nrows()) def write(self, data=None, /, **columns): - '''Writes to FITS by calling the internal _append method. + """Writes to FITS by calling the internal _append method. Pass either a positional variable (data) - or multiple named arguments (**columns)''' + or multiple named arguments (**columns)""" # if data is given, write it as it is if data is not None: @@ -101,6 +101,7 @@ def write_catalog(filename, *, ext=None): """ import fitsio + with fitsio.FITS(filename, "rw", clobber=True) as fits: fits.write(None) writer = _FitsWriter(fits, ext) diff --git a/pyproject.toml b/pyproject.toml index ddf2bd8d..96d6f4ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -88,5 +88,12 @@ testpaths = [ ] xfail_strict = true +[tool.ruff] +fix = true +force-exclude = true +src = [ + "glass", +] + [tool.tomlsort] overrides."project.classifiers".inline_arrays = false diff --git a/tests/core/test_array.py b/tests/core/test_array.py index 348e0048..e4065318 100644 --- a/tests/core/test_array.py +++ b/tests/core/test_array.py @@ -38,13 +38,11 @@ def test_ndinterp(): x = [[0.5, 1.5], [2.5, 3.5]] y = ndinterp(x, xp, yp) assert np.shape(y) == (2, 2) - npt.assert_allclose(y, [[1.15, 1.25], - [1.35, 1.45]], atol=1e-15) + npt.assert_allclose(y, [[1.15, 1.25], [1.35, 1.45]], atol=1e-15) # test nd interpolation in final axis - yp = [[1.1, 1.2, 1.3, 1.4, 1.5], - [2.1, 2.2, 2.3, 2.4, 2.5]] + yp = [[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]] x = 0.5 y = ndinterp(x, xp, yp) @@ -54,21 +52,18 @@ def test_ndinterp(): x = [0.5, 1.5, 2.5] y = ndinterp(x, xp, yp) assert np.shape(y) == (2, 3) - npt.assert_allclose(y, [[1.15, 1.25, 1.35], - [2.15, 2.25, 2.35]], atol=1e-15) + npt.assert_allclose(y, [[1.15, 1.25, 1.35], [2.15, 2.25, 2.35]], atol=1e-15) x = [[0.5, 1.5], [2.5, 3.5]] y = ndinterp(x, xp, yp) assert np.shape(y) == (2, 2, 2) - npt.assert_allclose(y, [[[1.15, 1.25], - [1.35, 1.45]], - [[2.15, 2.25], - [2.35, 2.45]]], atol=1e-15) + npt.assert_allclose( + y, [[[1.15, 1.25], [1.35, 1.45]], [[2.15, 2.25], [2.35, 2.45]]], atol=1e-15 + ) # test nd interpolation in middle axis - yp = [[[1.1], [1.2], [1.3], [1.4], [1.5]], - [[2.1], [2.2], [2.3], [2.4], [2.5]]] + yp = [[[1.1], [1.2], [1.3], [1.4], [1.5]], [[2.1], [2.2], [2.3], [2.4], [2.5]]] x = 0.5 y = ndinterp(x, xp, yp, axis=1) @@ -78,18 +73,29 @@ def test_ndinterp(): x = [0.5, 1.5, 2.5] y = ndinterp(x, xp, yp, axis=1) assert np.shape(y) == (2, 3, 1) - npt.assert_allclose(y, [[[1.15], [1.25], [1.35]], - [[2.15], [2.25], [2.35]]], atol=1e-15) + npt.assert_allclose( + y, [[[1.15], [1.25], [1.35]], [[2.15], [2.25], [2.35]]], atol=1e-15 + ) x = [[0.5, 1.5, 2.5, 3.5], [3.5, 2.5, 1.5, 0.5], [0.5, 3.5, 1.5, 2.5]] y = ndinterp(x, xp, yp, axis=1) assert np.shape(y) == (2, 3, 4, 1) - npt.assert_allclose(y, [[[[1.15], [1.25], [1.35], [1.45]], - [[1.45], [1.35], [1.25], [1.15]], - [[1.15], [1.45], [1.25], [1.35]]], - [[[2.15], [2.25], [2.35], [2.45]], - [[2.45], [2.35], [2.25], [2.15]], - [[2.15], [2.45], [2.25], [2.35]]]], atol=1e-15) + npt.assert_allclose( + y, + [ + [ + [[1.15], [1.25], [1.35], [1.45]], + [[1.45], [1.35], [1.25], [1.15]], + [[1.15], [1.45], [1.25], [1.35]], + ], + [ + [[2.15], [2.25], [2.35], [2.45]], + [[2.45], [2.35], [2.25], [2.15]], + [[2.15], [2.45], [2.25], [2.35]], + ], + ], + atol=1e-15, + ) def test_trapz_product(): diff --git a/tests/test_fields.py b/tests/test_fields.py index d273a02a..18bbad2a 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -1,5 +1,6 @@ def test_getcl(): from glass.fields import getcl + # make a mock Cls array with the index pairs as entries cls = [{i, j} for i in range(10) for j in range(i, -1, -1)] # make sure indices are retrieved correctly diff --git a/tests/test_fits.py b/tests/test_fits.py index 25549388..055c0b33 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -2,6 +2,7 @@ # check if fitsio is available for testing import importlib.util + if importlib.util.find_spec("fitsio") is not None: HAVE_FITSIO = True else: @@ -12,8 +13,8 @@ def _test_append(fits, data, names): - '''Write routine for FITS test cases''' - cat_name = 'CATALOG' + """Write routine for FITS test cases""" + cat_name = "CATALOG" if cat_name not in fits: fits.write_table(data, names=names, extname=cat_name) else: @@ -30,25 +31,30 @@ def _test_append(fits, data, names): @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio") def test_basic_write(tmp_path): import fitsio + d = tmp_path / "sub" d.mkdir() filename_gfits = "gfits.fits" # what GLASS creates filename_tfits = "tfits.fits" # file created on the fly to test against - with user.write_catalog(d / filename_gfits, ext="CATALOG") as out, fitsio.FITS(d / filename_tfits, "rw", clobber=True) as myFits: + with user.write_catalog(d / filename_gfits, ext="CATALOG") as out, fitsio.FITS( + d / filename_tfits, "rw", clobber=True + ) as myFits: for i in range(0, my_max): - array = np.arange(i, i+1, delta) # array of size 1/delta - array2 = np.arange(i+1, i+2, delta) # array of size 1/delta + array = np.arange(i, i + 1, delta) # array of size 1/delta + array2 = np.arange(i + 1, i + 2, delta) # array of size 1/delta out.write(RA=array, RB=array2) arrays = [array, array2] - names = ['RA', 'RB'] + names = ["RA", "RB"] _test_append(myFits, arrays, names) - with fitsio.FITS(d / filename_gfits) as g_fits, fitsio.FITS(d / filename_tfits) as t_fits: + with fitsio.FITS(d / filename_gfits) as g_fits, fitsio.FITS( + d / filename_tfits + ) as t_fits: glass_data = g_fits[1].read() test_data = t_fits[1].read() - assert glass_data['RA'].size == test_data['RA'].size - assert glass_data['RB'].size == test_data['RA'].size + assert glass_data["RA"].size == test_data["RA"].size + assert glass_data["RB"].size == test_data["RA"].size @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio") @@ -61,22 +67,25 @@ def test_write_exception(tmp_path): for i in range(0, my_max): if i == except_int: raise Exception("Unhandled exception") - array = np.arange(i, i+1, delta) # array of size 1/delta - array2 = np.arange(i+1, i+2, delta) # array of size 1/delta + array = np.arange(i, i + 1, delta) # array of size 1/delta + array2 = np.arange(i + 1, i + 2, delta) # array of size 1/delta out.write(RA=array, RB=array2) except Exception: import fitsio + with fitsio.FITS(d / filename) as hdul: data = hdul[1].read() - assert data['RA'].size == except_int/delta - assert data['RB'].size == except_int/delta + assert data["RA"].size == except_int / delta + assert data["RB"].size == except_int / delta - fitsMat = data['RA'].reshape(except_int, int(1/delta)) - fitsMat2 = data['RB'].reshape(except_int, int(1/delta)) + fitsMat = data["RA"].reshape(except_int, int(1 / delta)) + fitsMat2 = data["RB"].reshape(except_int, int(1 / delta)) for i in range(0, except_int): - array = np.arange(i, i+1, delta) # re-create array to compare to read data - array2 = np.arange(i+1, i+2, delta) + array = np.arange( + i, i + 1, delta + ) # re-create array to compare to read data + array2 = np.arange(i + 1, i + 2, delta) assert array.tolist() == fitsMat[i].tolist() assert array2.tolist() == fitsMat2[i].tolist() @@ -84,6 +93,7 @@ def test_write_exception(tmp_path): @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio") def test_out_filename(tmp_path): import fitsio + fits = fitsio.FITS(filename, "rw", clobber=True) writer = user._FitsWriter(fits) assert writer.fits._filename == filename diff --git a/tests/test_galaxies.py b/tests/test_galaxies.py index c852b7d0..3520df85 100644 --- a/tests/test_galaxies.py +++ b/tests/test_galaxies.py @@ -8,13 +8,13 @@ def test_redshifts(): # create a mock radial window function w = Mock() - w.za = np.linspace(0., 1., 20) - w.wa = np.exp(-0.5*(w.za - 0.5)**2/0.1**2) + w.za = np.linspace(0.0, 1.0, 20) + w.wa = np.exp(-0.5 * (w.za - 0.5) ** 2 / 0.1**2) # sample redshifts (scalar) z = redshifts(13, w) assert z.shape == (13,) - assert z.min() >= 0. and z.max() <= 1. + assert z.min() >= 0.0 and z.max() <= 1.0 # sample redshifts (array) z = redshifts([[1, 2], [3, 4]], w) @@ -45,7 +45,7 @@ def test_redshifts_from_nz(): count = 10 z = np.linspace(0, 1, 100) - nz = z*(1-z) + nz = z * (1 - z) redshifts = redshifts_from_nz(count, z, nz) @@ -56,7 +56,7 @@ def test_redshifts_from_nz(): count = [10, 20, 30] z = np.linspace(0, 1, 100) - nz = z*(1-z) + nz = z * (1 - z) redshifts = redshifts_from_nz(count, z, nz) @@ -66,7 +66,7 @@ def test_redshifts_from_nz(): count = 10 z = np.linspace(0, 1, 100) - nz = [z*(1-z), (z-0.5)**2] + nz = [z * (1 - z), (z - 0.5) ** 2] redshifts = redshifts_from_nz(count, z, nz) @@ -76,7 +76,7 @@ def test_redshifts_from_nz(): count = [[10], [20], [30]] z = np.linspace(0, 1, 100) - nz = [z*(1-z), (z-0.5)**2] + nz = [z * (1 - z), (z - 0.5) ** 2] redshifts = redshifts_from_nz(count, z, nz) @@ -86,7 +86,7 @@ def test_redshifts_from_nz(): count = [10, 20, 30] z = np.linspace(0, 1, 100) - nz = [z*(1-z), (z-0.5)**2] + nz = [z * (1 - z), (z - 0.5) ** 2] with pytest.raises(ValueError): redshifts_from_nz(count, z, nz) @@ -101,7 +101,7 @@ def test_gaussian_phz(): # case: zero variance z = np.linspace(0, 1, 100) - sigma_0 = 0. + sigma_0 = 0.0 phz = gaussian_phz(z, sigma_0) @@ -109,7 +109,7 @@ def test_gaussian_phz(): # case: truncated normal - z = 0. + z = 0.0 sigma_0 = np.ones(100) phz = gaussian_phz(z, sigma_0) @@ -119,7 +119,7 @@ def test_gaussian_phz(): # case: upper and lower bound - z = 1. + z = 1.0 sigma_0 = np.ones(100) phz = gaussian_phz(z, sigma_0, lower=0.5, upper=1.5) @@ -132,8 +132,8 @@ def test_gaussian_phz(): # case: scalar redshift, scalar sigma_0 - z = 1. - sigma_0 = 0. + z = 1.0 + sigma_0 = 0.0 phz = gaussian_phz(z, sigma_0) @@ -143,7 +143,7 @@ def test_gaussian_phz(): # case: array redshift, scalar sigma_0 z = np.linspace(0, 1, 10) - sigma_0 = 0. + sigma_0 = 0.0 phz = gaussian_phz(z, sigma_0) @@ -152,7 +152,7 @@ def test_gaussian_phz(): # case: scalar redshift, array sigma_0 - z = 1. + z = 1.0 sigma_0 = np.zeros(10) phz = gaussian_phz(z, sigma_0) diff --git a/tests/test_lensing.py b/tests/test_lensing.py index da0780de..d1c60270 100644 --- a/tests/test_lensing.py +++ b/tests/test_lensing.py @@ -8,11 +8,11 @@ def shells(): from glass.shells import RadialWindow shells = [ - RadialWindow([0., 1., 2.], [0., 1., 0.], 1.), - RadialWindow([1., 2., 3.], [0., 1., 0.], 2.), - RadialWindow([2., 3., 4.], [0., 1., 0.], 3.), - RadialWindow([3., 4., 5.], [0., 1., 0.], 4.), - RadialWindow([4., 5., 6.], [0., 1., 0.], 5.), + RadialWindow([0.0, 1.0, 2.0], [0.0, 1.0, 0.0], 1.0), + RadialWindow([1.0, 2.0, 3.0], [0.0, 1.0, 0.0], 2.0), + RadialWindow([2.0, 3.0, 4.0], [0.0, 1.0, 0.0], 3.0), + RadialWindow([3.0, 4.0, 5.0], [0.0, 1.0, 0.0], 4.0), + RadialWindow([4.0, 5.0, 6.0], [0.0, 1.0, 0.0], 5.0), ] return shells @@ -21,13 +21,12 @@ def shells(): @pytest.fixture def cosmo(): class MockCosmology: - @property def omega_m(self): return 0.3 def ef(self, z): - return (self.omega_m * (1 + z)**3 + 1 - self.omega_m) ** 0.5 + return (self.omega_m * (1 + z) ** 3 + 1 - self.omega_m) ** 0.5 def xm(self, z, z2=None): if z2 is None: @@ -42,31 +41,33 @@ def xm(self, z, z2=None): def test_deflect_nsew(usecomplex): from glass.lensing import deflect - d = 5. + d = 5.0 r = np.radians(d) if usecomplex: + def alpha(re, im): - return re + 1j*im + return re + 1j * im else: + def alpha(re, im): return [re, im] # north - lon, lat = deflect(0., 0., alpha(r, 0)) - assert np.allclose([lon, lat], [0., d]) + lon, lat = deflect(0.0, 0.0, alpha(r, 0)) + assert np.allclose([lon, lat], [0.0, d]) # south - lon, lat = deflect(0., 0., alpha(-r, 0)) - assert np.allclose([lon, lat], [0., -d]) + lon, lat = deflect(0.0, 0.0, alpha(-r, 0)) + assert np.allclose([lon, lat], [0.0, -d]) # east - lon, lat = deflect(0., 0., alpha(0, r)) - assert np.allclose([lon, lat], [-d, 0.]) + lon, lat = deflect(0.0, 0.0, alpha(0, r)) + assert np.allclose([lon, lat], [-d, 0.0]) # west - lon, lat = deflect(0., 0., alpha(0, -r)) - assert np.allclose([lon, lat], [d, 0.]) + lon, lat = deflect(0.0, 0.0, alpha(0, -r)) + assert np.allclose([lon, lat], [d, 0.0]) def test_deflect_many(): @@ -74,18 +75,18 @@ def test_deflect_many(): from glass.lensing import deflect n = 1000 - abs_alpha = np.random.uniform(0, 2*np.pi, size=n) + abs_alpha = np.random.uniform(0, 2 * np.pi, size=n) arg_alpha = np.random.uniform(-np.pi, np.pi, size=n) lon_ = np.degrees(np.random.uniform(-np.pi, np.pi, size=n)) lat_ = np.degrees(np.arcsin(np.random.uniform(-1, 1, size=n))) - lon, lat = deflect(lon_, lat_, abs_alpha*np.exp(1j*arg_alpha)) + lon, lat = deflect(lon_, lat_, abs_alpha * np.exp(1j * arg_alpha)) x_, y_, z_ = healpix.ang2vec(lon_, lat_, lonlat=True) x, y, z = healpix.ang2vec(lon, lat, lonlat=True) - dotp = x*x_ + y*y_ + z*z_ + dotp = x * x_ + y * y_ + z * z_ npt.assert_allclose(dotp, np.cos(abs_alpha)) @@ -130,4 +131,4 @@ def test_multi_plane_weights(shells, cosmo): wmat = multi_plane_weights(weights, shells, cosmo) - npt.assert_allclose(np.einsum('ij,ik', wmat, deltas), kappa) + npt.assert_allclose(np.einsum("ij,ik", wmat, deltas), kappa) diff --git a/tests/test_points.py b/tests/test_points.py index fd4f9bf6..9247ec02 100644 --- a/tests/test_points.py +++ b/tests/test_points.py @@ -102,7 +102,6 @@ def test_position_weights(): for bshape in None, (), (100,), (100, 1): for cshape in (100,), (100, 50), (100, 3, 2): - counts = np.random.rand(*cshape) bias = None if bshape is None else np.random.rand(*bshape) @@ -111,9 +110,13 @@ def test_position_weights(): expected = counts / counts.sum(axis=0, keepdims=True) if bias is not None: if np.ndim(bias) > np.ndim(expected): - expected = np.expand_dims(expected, tuple(range(np.ndim(expected), np.ndim(bias)))) + expected = np.expand_dims( + expected, tuple(range(np.ndim(expected), np.ndim(bias))) + ) else: - bias = np.expand_dims(bias, tuple(range(np.ndim(bias), np.ndim(expected)))) + bias = np.expand_dims( + bias, tuple(range(np.ndim(bias), np.ndim(expected))) + ) expected = bias * expected npt.assert_allclose(weights, expected) diff --git a/tests/test_shapes.py b/tests/test_shapes.py index 3a4100c0..1d3d36af 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -3,7 +3,6 @@ def test_ellipticity_gaussian(): - from glass.shapes import ellipticity_gaussian n = 1_000_000 @@ -19,7 +18,7 @@ def test_ellipticity_gaussian(): eps = ellipticity_gaussian([n, n], [0.128, 0.256]) - assert eps.shape == (2*n,) + assert eps.shape == (2 * n,) assert np.all(np.abs(eps) < 1) @@ -30,7 +29,6 @@ def test_ellipticity_gaussian(): def test_ellipticity_intnorm(): - from glass.shapes import ellipticity_intnorm n = 1_000_000 @@ -46,7 +44,7 @@ def test_ellipticity_intnorm(): eps = ellipticity_intnorm([n, n], [0.128, 0.256]) - assert eps.shape == (2*n,) + assert eps.shape == (2 * n,) assert np.all(np.abs(eps) < 1) diff --git a/tests/test_shells.py b/tests/test_shells.py index 10f76318..ff92497d 100644 --- a/tests/test_shells.py +++ b/tests/test_shells.py @@ -10,13 +10,13 @@ def test_tophat_windows(): ws = tophat_windows(zb, dz) - assert len(ws) == len(zb)-1 + assert len(ws) == len(zb) - 1 - assert all(z0 == w.za[0] and zn == w.za[-1] - for w, z0, zn in zip(ws, zb, zb[1:])) + assert all(z0 == w.za[0] and zn == w.za[-1] for w, z0, zn in zip(ws, zb, zb[1:])) - assert all(zn <= z0+len(w.za)*dz <= zn+dz - for w, z0, zn in zip(ws, zb, zb[1:])) + assert all( + zn <= z0 + len(w.za) * dz <= zn + dz for w, z0, zn in zip(ws, zb, zb[1:]) + ) assert all(np.all(w.wa == 1) for w in ws) @@ -25,28 +25,28 @@ def test_restrict(): from glass.shells import restrict, RadialWindow # Gaussian test function - z = np.linspace(0., 5., 1000) - f = np.exp(-((z - 2.)/0.5)**2/2) + z = np.linspace(0.0, 5.0, 1000) + f = np.exp(-(((z - 2.0) / 0.5) ** 2) / 2) # window for restriction - w = RadialWindow(za=[1., 2., 3., 4.], wa=[0., .5, .5, 0.], zeff=None) + w = RadialWindow(za=[1.0, 2.0, 3.0, 4.0], wa=[0.0, 0.5, 0.5, 0.0], zeff=None) zr, fr = restrict(z, f, w) assert zr[0] == w.za[0] and zr[-1] == w.za[-1] - assert fr[0] == fr[-1] == 0. + assert fr[0] == fr[-1] == 0.0 for zi, wi in zip(w.za, w.wa): i = np.searchsorted(zr, zi) assert zr[i] == zi - assert fr[i] == wi*np.interp(zi, z, f) + assert fr[i] == wi * np.interp(zi, z, f) for zi, fi in zip(z, f): if w.za[0] <= zi <= w.za[-1]: i = np.searchsorted(zr, zi) assert zr[i] == zi - assert fr[i] == fi*np.interp(zi, w.za, w.wa) + assert fr[i] == fi * np.interp(zi, w.za, w.wa) @pytest.mark.parametrize("method", ["lstsq", "nnls", "restrict"]) @@ -55,15 +55,15 @@ def test_partition(method): from glass.shells import RadialWindow, partition shells = [ - RadialWindow(np.array([0., 1.]), np.array([1., 0.]), 0.0), - RadialWindow(np.array([0., 1., 2.]), np.array([0., 1., 0.]), 0.5), - RadialWindow(np.array([1., 2., 3.]), np.array([0., 1., 0.]), 1.5), - RadialWindow(np.array([2., 3., 4.]), np.array([0., 1., 0.]), 2.5), - RadialWindow(np.array([3., 4., 5.]), np.array([0., 1., 0.]), 3.5), - RadialWindow(np.array([4., 5.]), np.array([0., 1.]), 5.0), + RadialWindow(np.array([0.0, 1.0]), np.array([1.0, 0.0]), 0.0), + RadialWindow(np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 0.0]), 0.5), + RadialWindow(np.array([1.0, 2.0, 3.0]), np.array([0.0, 1.0, 0.0]), 1.5), + RadialWindow(np.array([2.0, 3.0, 4.0]), np.array([0.0, 1.0, 0.0]), 2.5), + RadialWindow(np.array([3.0, 4.0, 5.0]), np.array([0.0, 1.0, 0.0]), 3.5), + RadialWindow(np.array([4.0, 5.0]), np.array([0.0, 1.0]), 5.0), ] - z = np.linspace(0., 5., 1000) + z = np.linspace(0.0, 5.0, 1000) k = 1 + np.arange(6).reshape(3, 2, 1) fz = np.exp(-z / k)