Repository: skypyproject/skypy
Pull request #581 (2022-11-21 16:43:54)

MAINT: merge main into module/halos
## Description

## Checklist
- [ ] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
- [ ] Write unit tests
- [ ] Write documentation strings
- [ ] Assign someone from your working team to review this pull request
- [ ] Assign someone from the infrastructure team to review this pull request



Patch:
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
index 05b4e83..c6aac7e 100644
--- a/.github/PULL_REQUEST_TEMPLATE.md
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -1,7 +1,7 @@
 ## Description
 
 ## Checklist
-- [ ] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/master/CONTRIBUTING.md)
+- [ ] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
 - [ ] Write unit tests
 - [ ] Write documentation strings
 - [ ] Assign someone from your working team to review this pull request
diff --git a/.github/workflows/codestyle.yaml b/.github/workflows/codestyle.yaml
index cdc8583..2e14de4 100644
--- a/.github/workflows/codestyle.yaml
+++ b/.github/workflows/codestyle.yaml
@@ -2,11 +2,11 @@ name: Code Style
 on:
   push:
     branches:
-    - master
+    - main
     - module/*
   pull_request:
     branches:
-    - master
+    - main
     - module/*
 jobs:
   flake8:
diff --git a/.github/workflows/compatibility.yaml b/.github/workflows/compatibility.yaml
index 4a4c4d9..2fd34cd 100644
--- a/.github/workflows/compatibility.yaml
+++ b/.github/workflows/compatibility.yaml
@@ -9,15 +9,15 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest]
-        python: [3.9]
-        toxenv: [py39-test, py39-test-dev]
-        release: [master, latest]
+        python: ['3.10']
+        toxenv: [py310-test, py310-test-dev]
+        release: [main, latest]
     steps:
     - name: Checkout Repository
       uses: actions/checkout@v2
       with:
         fetch-depth: 0
-    - if: matrix.release != 'master'
+    - if: matrix.release != 'main'
       name: Checkout Release
       run: |
         git checkout tags/$(curl -s https://api.github.com/repos/skypyproject/skypy/releases/${{ matrix.release }} | python -c "import sys, json; print(json.load(sys.stdin)['tag_name'])")
diff --git a/.mailmap b/.mailmap
index 21a2dc8..3dc76d1 100644
--- a/.mailmap
+++ b/.mailmap
@@ -13,5 +13,6 @@
 Adam Amara <adam.amara@phys.ethz.ch> adamamara <adamamara@penguin>
 Coleman Krawczyk <coleman.krawczyk@gmail.com> CKrawczyk
 Ian Harry <ian.harry@ligo.org> <iwharry@googlemail.com>
+Lucia F. de la Bella <55983939+Lucia-Fonseca@users.noreply.github.com>
 Nicolas Tessore <n.tessore@ucl.ac.uk> <nicolas.tessore@manchester.ac.uk>
 Richard R <58728519+rrjbca@users.noreply.github.com> <rrjbca@users.noreply.github.com>
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 46f0150..3369dad 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -1,17 +1,24 @@
 version: 2
 
 build:
-  image: latest
+  os: ubuntu-20.04
+  apt_packages:
+    - graphviz
+  tools:
+    python: "3.7"
 
 sphinx:
+  builder: html
+  configuration: docs/conf.py
   fail_on_warning: true
 
 python:
-  version: 3.7
+  system_packages: false
   install:
     - method: pip
       path: .
       extra_requirements:
         - docs
+        - all
 
 formats: []
diff --git a/.zenodo.json b/.zenodo.json
index 3a0998b..ef09137 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -25,7 +25,7 @@
         },
         {
             "orcid": "0000-0002-1064-3400",
-            "affiliation": "University of Manchester",
+            "affiliation": "University of Portsmouth",
             "name": "Lucia F. de la Bella"
         },
         {
@@ -43,6 +43,11 @@
             "affiliation": "University of Manchester",
             "name": "Juan Pablo Cordero"
         },
+        {
+            "orcid": "0000-0003-1560-7959",
+            "affiliation": "University of Portsmouth",
+            "name": "Fox Davidson"
+        },
         {
             "orcid": "0000-0002-8218-563X",
             "affiliation": "University of Portsmouth",
@@ -50,7 +55,7 @@
         },
         {
             "orcid": "0000-0002-4437-0770",
-            "affiliation": "University of Manchester",
+            "affiliation": "University of Oxford, University of Manchester",
             "name": "Ian Harrison"
         },
         {
@@ -58,6 +63,11 @@
             "affiliation": "University of Portsmouth",
             "name": "Ian Harry"
         },
+        {
+            "orcid": "0000-0001-9994-1115",
+            "affiliation": "University of Geneva",
+            "name": "William G. Hartley"
+        },
         {
             "orcid": "0000-0001-9233-2341",
             "affiliation": "University of Portsmouth",
@@ -80,14 +90,19 @@
         },
         {
             "orcid": "0000-0003-1291-1023",
-            "affiliation": "University of Manchester",
-            "name": "Richard Rollins"
+            "affiliation": "University of Edinburgh",
+            "name": "Richard P. Rollins"
         },
         {
             "orcid": "0000-0001-8685-2308",
             "affiliation": "University of Portsmouth",
             "name": "Philipp Sudek"
         },
+        {
+            "orcid": "0000-0002-6724-833X",
+            "affiliation": "Institute of Astronomy and Astrophysics Academia Sinica",
+            "name": "Sut-Ieng Tam"
+        },
         {
             "orcid": "0000-0002-9696-7931",
             "affiliation": "UCL",
@@ -98,10 +113,20 @@
             "affiliation": "Institute of Astronomy and Astrophysics Academia Sinica",
             "name": "Keiichi Umetsu"
         },
+        {
+            "orcid": "0000-0001-9841-943X",
+            "affiliation": "University of Portsmouth",
+            "name": "Arthur E. Tolley"
+        },
         {
             "orcid": "0000-0002-7627-8688",
             "affiliation": "University of Portsmouth",
             "name": "Andrew R. Williamson"
+        },
+        {
+            "orcid": "0000-0003-3334-3037",
+            "affiliation": "University of Manchester",
+            "name": "Laura Wolz"
         }
     ],
     "contributors": [
diff --git a/CITATION.rst b/CITATION.rst
index 04e06d0..f5a4565 100644
--- a/CITATION.rst
+++ b/CITATION.rst
@@ -1,13 +1,27 @@
+Citation Guidelines
+===================
+
+|JOSS| |Zenodo|
+
+
 If you use SkyPy for work or research presented in a publication (whether
-directly, or as a dependency to another package) we recommend and encourage
+directly, or as a dependency of another package) we recommend and encourage
 the following acknowledgment:
 
   This research made use of SkyPy, a Python package for forward modeling
-  astronomical surveys (SkyPy Collaboration, 2020).
+  astronomical surveys (Amara et. al., 2021, SkyPy Collaboration, 202x).
+
+where the citations are to our publication in the `Journal of Open Source
+Software`_ and the `Zenodo DOI`_ for the specific version of the software that
+you used. We also encourage citations within the main text wherever
+appropriate. DOIs and BibTeX keys are available through the links above.
+
+.. _Journal of Open Source Software: https://joss.theoj.org/papers/10.21105/joss.03056
+.. _Zenodo DOI: https://zenodo.org/record/3755531
+
 
-where (SkyPy Collaboration, 2020) is a citation to the software DOI for the
-specific version that you used. We also encourage you to cite the software DOI
-in the main text wherever appropriate. DOIs and BibTeX keys for each release
-can be found on Zenodo_.
+.. |JOSS| image:: https://joss.theoj.org/papers/10.21105/joss.03056/status.svg
+    :target: https://doi.org/10.21105/joss.03056
 
-.. _Zenodo: https://zenodo.org/record/3755531
+.. |Zenodo| image:: https://zenodo.org/badge/doi/10.5281/zenodo.4475347.svg
+    :target: https://doi.org/10.5281/zenodo.3755531
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 767689b..758a3d0 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -32,12 +32,12 @@ Finally add the ``skypyproject`` repository as a *remote*. This will allow you t
 Create a branch for your new feature
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-Create a *branch* off the ``skypyproject`` master branch. Working on unique branches for each new feature simplifies the development, review and merge processes by maintining logical separation. To create a feature branch:
+Create a *branch* off the ``skypyproject`` main branch. Working on unique branches for each new feature simplifies the development, review and merge processes by maintining logical separation. To create a feature branch:
 
 ::
 
   git fetch skypyproject
-  git checkout -b <your-branch-name> skypyproject/master
+  git checkout -b <your-branch-name> skypyproject/main
 
 
 Hack away!
@@ -73,7 +73,7 @@ When you feel that work on your new feature is complete, you should create a *Pu
   1. Go to `SkyPy Pull Requests <https://github.com/skypyproject/skypy/pulls>`_
   2. Click the green **New pull request** button
   3. Click **compare across forks**
-  4. Confirm that the base fork is ``skypyproject/skypy`` and the base branch is ``master``
+  4. Confirm that the base fork is ``skypyproject/skypy`` and the base branch is ``main``
   5. Confirm the head fork is ``<your-account>/skypy`` and the compare branch is ``<your-branch-name>``
   6. Give your pull request a title and fill out the the template for the description
   7. Click the green **Create pull request** button
@@ -91,22 +91,22 @@ A series of automated checks will be run on your pull request, some of which wil
 Updating your branch
 ^^^^^^^^^^^^^^^^^^^^
 
-As you work on your feature, new commits might be made to the ``skypyproject`` master branch. You will need to update your branch with these new commits before your pull request can be accepted. You can achieve this in a few different ways:
+As you work on your feature, new commits might be made to the ``skypyproject`` main branch. You will need to update your branch with these new commits before your pull request can be accepted. You can achieve this in a few different ways:
 
   - If your pull request has no conflicts, click **Update branch**
   - If your pull request has conflicts, click **Resolve conflicts**, manually resolve the conflicts and click **Mark as resolved**
-  - *merge* the ``skypyproject`` master branch from the command line:
+  - *merge* the ``skypyproject`` main branch from the command line:
 
     ::
 
         git fetch skypyproject
-        git merge skypyproject/master
+        git merge skypyproject/main
 
-  - *rebase* your feature branch onto the ``skypyproject`` master branch from the command line:
+  - *rebase* your feature branch onto the ``skypyproject`` main branch from the command line:
     ::
 
         git fetch skypyproject
-        git rebase skypyproject/master
+        git rebase skypyproject/main
 
 
 **Warning**: It is bad practice to *rebase* commits that have already been pushed to a remote such as your fork. Rebasing creates new copies of your commits that can cause the local and remote branches to diverge. ``git push --force`` will **overwrite** the remote branch with your newly rebased local branch. This is strongly discouraged, particularly when working on a shared branch where you could erase a collaborators commits.
@@ -129,7 +129,7 @@ Before your pull request can be merged into the codebase, it will be reviewed by
 General Guidelines
 ^^^^^^^^^^^^^^^^^^
 
-- SkyPy is compatible with Python>=3.6 (see `setup.cfg <https://github.com/skypyproject/skypy/blob/master/setup.cfg>`_). SkyPy *does not* support backwards compatibility with Python 2.x; `six`, `__future__` and `2to3` should not be used.
+- SkyPy is compatible with Python>=3.7 (see `setup.cfg <https://github.com/skypyproject/skypy/blob/main/setup.cfg>`_). SkyPy *does not* support backwards compatibility with Python 2.x; `six`, `__future__` and `2to3` should not be used.
 - All contributions should follow the `PEP8 Style Guide for Python Code <https://www.python.org/dev/peps/pep-0008/>`_. We recommend using `flake8 <https://flake8.pycqa.org/>`__ to check your code for PEP8 compliance.
 - Importing SkyPy should only depend on having `NumPy <https://www.numpy.org>`_, `SciPy <https://www.scipy.org/>`_ and `Astropy <https://www.astropy.org/>`__ installed.
 - Code is grouped into submodules based on broad science areas e.g. `galaxies <https://skypy.readthedocs.io/en/stable/galaxies.html>`_. There is also a `utils <https://skypy.readthedocs.io/en/stable/utils/index.html>`_ submodule for general utility functions.
@@ -150,7 +150,6 @@ All public classes, methods and functions require docstrings. You can build docu
   - Description
   - Parameters
   - Notes
-  - Examples
   - References
 
 For more information see the Astropy guide to `Writing Documentation <https://docs.astropy.org/en/stable/development/docguide.html>`_.
diff --git a/README.rst b/README.rst
index 1869bd9..9e02852 100644
--- a/README.rst
+++ b/README.rst
@@ -2,32 +2,56 @@
 SkyPy: A package for modelling the Universe
 ===========================================
 
-|PyPI| |conda-forge| |Read the Docs| |GitHub| |Codecov| |Zenodo|
+|Read the Docs| |GitHub| |Codecov| |Compatibility|
 
-This package contains methods for modelling the Universe, galaxies and the Milky
-Way. Also included are methods for generating observed data.
+This package contains methods for modelling the Universe, galaxies and the
+Milky Way. SkyPy simulates populations of astronomical objects, generating
+random realisations of intrinsic and observed properties, with the
+intention the simulations can then be compared to data as part of an inference
+pipeline.
+
+Currently, SkyPy implements the following modules:
 
 * Galaxies_: morphology, luminosity and redshift distributions
 * Halo_ and subhalo mass distributions
 * `Power Spectra`_ using CAMB and Halofit
 * Pipelines_ to generate populations of astronomical objects
 
-The full list of features can be found in the `SkyPy Documentation`_.
+The `full list of features`_ can be found in the `SkyPy Documentation`_.
+
+For more information on the people involved and how SkyPy is developed, please
+visit the SkyPy Collaboration website: `http://skypyproject.org`_
+
+.. _Galaxies: https://skypy.readthedocs.io/en/latest/galaxies.html
+.. _Pipelines: https://skypy.readthedocs.io/en/latest/pipeline/index.html
+.. _full list of features: https://skypy.readthedocs.io/en/latest/feature_list.html
+.. _SkyPy Documentation: https://skypy.readthedocs.io/en/latest/
+.. _http://skypyproject.org: http://skypyproject.org
+
+Citation
+--------
+
+|JOSS| |Zenodo|
 
 If you use SkyPy for work or research presented in a publication please follow
 our `Citation Guidelines`_.
 
+<<<<<<< HEAD
 .. _Galaxies: https://skypy.readthedocs.io/en/latest/galaxies.html
 .. _Halo: https://skypy.readthedocs.io/en/latest/halos/index.html
 .. _Power Spectra: https://skypy.readthedocs.io/en/latest/power_spectrum/index.html
 .. _Pipelines: https://skypy.readthedocs.io/en/latest/pipeline/index.html
 .. _SkyPy Documentation: https://skypy.readthedocs.io/en/latest/
+=======
+>>>>>>> main
 .. _Citation Guidelines: CITATION.rst
 
 
 Installation
 ------------
 
+|PyPI| |conda-forge|
+
 SkyPy releases are distributed through PyPI_ and conda-forge_. Instructions for
 installing SkyPy and its dependencies can be found in the Installation_
 section of the documentation.
@@ -47,6 +71,16 @@ Examples_ that demonstrate how to use it.
 .. _Examples: https://skypy.readthedocs.io/en/stable/examples/index.html
 
 
+Get in Touch
+------------
+
+You are welcome to talk about the SkyPy package and code using our
+`Discussions Page`_. For any other questions about the project in general,
+please get in touch with the `SkyPy Co-ordinators`_.
+
+ .. _Discussions Page: https://github.com/skypyproject/skypy/discussions
+ .. _SkyPy Co-ordinators: mailto:skypy-coordinators@googlegroups.com
+
 Contributing
 ------------
 
@@ -56,8 +90,8 @@ For information on how to contribute see our `Contributor Guidelines`_.
 All communication relating to The SkyPy Project must meet the standards set out
 in the `Code of Conduct`_.
 
-.. _Contributor Guidelines: https://skypy.readthedocs.io/en/stable/developer/contributing.html
-.. _Code of Conduct: https://skypy.readthedocs.io/en/stable/project/CODE_OF_CONDUCT.html
+.. _Contributor Guidelines: https://skypy.readthedocs.io/en/latest/developer/contributing.html
+.. _Code of Conduct: https://skypy.readthedocs.io/en/stable/project/code_of_conduct.html
 
 .. |PyPI| image:: https://img.shields.io/pypi/v/skypy?label=PyPI&logo=pypi
     :target: https://pypi.python.org/pypi/skypy
@@ -71,9 +105,15 @@ in the `Code of Conduct`_.
 .. |GitHub| image:: https://github.com/skypyproject/skypy/workflows/Tests/badge.svg
     :target: https://github.com/skypyproject/skypy/actions
 
-.. |Codecov| image:: https://codecov.io/gh/skypyproject/skypy/branch/master/graph/badge.svg
+.. |Compatibility| image:: https://github.com/skypyproject/skypy/actions/workflows/compatibility.yaml/badge.svg
+    :target: https://github.com/skypyproject/skypy/actions/workflows/compatibility.yaml
+
+.. |Codecov| image:: https://codecov.io/gh/skypyproject/skypy/branch/main/graph/badge.svg
     :target: https://codecov.io/gh/skypyproject/skypy
 
 .. |Zenodo| image:: https://zenodo.org/badge/221432358.svg
     :target: https://zenodo.org/badge/latestdoi/221432358
     :alt: SkyPy Concept DOI
+
+.. |JOSS| image:: https://joss.theoj.org/papers/d4fac0604318190d6627ab29b568a48d/status.svg
+    :target: https://joss.theoj.org/papers/d4fac0604318190d6627ab29b568a48d
diff --git a/docs/arch/adr-03.md b/docs/arch/adr-03.md
new file mode 100644
index 0000000..809d115
--- /dev/null
+++ b/docs/arch/adr-03.md
@@ -0,0 +1,120 @@
+
+ADR 3: Position sampling and patches of the sky
+===============================================
+
+Author: Nicolas Tessore  
+Date: 5 February 2021  
+Status: Accepted
+
+
+Context
+-------
+
+This ADR addresses two related open problems for SkyPy simulations:
+
+- How can we describe the regions in which to sample the positions of e.g.
+  galaxies or supernovae?
+- How can we break up the sampling over large fractions of the sky into smaller
+  chunks?
+
+The second point is particularly relevant when the simulation becomes so large
+that it no longer fits into memory, as well as for parallelisation.
+
+
+Extension to pipelines
+----------------------
+
+This ADR proposes to introduce a new top-level keyword `regions` [alt: `sky`,
+`geometry`, `patches`] into pipelines that describes the geometry of the
+simulation. For example, to include a "rectangular" and a "circular" region:
+
+```yaml
+regions:
+- !rectangular [ ... ]
+- !circular [ ... ]
+```
+
+The `regions` list can be traversed by the pipeline runner to create what are
+effectively independent parallel simulations. The list items are objects with
+the following interface.
+
+
+Region interface
+----------------
+
+The regions need to support two sets of operations:
+
+- Information queries: For example, there should be a `.area` [alt:
+  `.solid_angle`] attribute that returns the solid angle of the region.
+- Random sampling: There needs to be at least a `random_point()` [alt:
+  `random_position()`, `uniform_in()`] function that can uniformly sample a
+  random position from within the region.
+
+When the pipeline runner traverses the list of regions, it can keep track of
+the current region in a `$region` reference that can be used where necessary.
+For example, to sample from a luminosity function with positions:
+
+```yaml
+tables:
+  galaxies:
+    z, M: !schechter_lf
+      ...
+      sky_area: $region.area
+    ra, dec: !random_point [ $region ]
+```
+
+
+Support for HEALPix maps
+------------------------
+
+The above proposal is powerful enough to support advanced features such as
+regions that are described by HEALPix maps. There may be a `healpix()` function
+that generates a list of regions from HEALPix pixels:
+
+```yaml
+regions: !healpix
+  nside: 8
+```
+
+The resulting list would contain `12 * nside**2 = 768` regions corresponding
+to the HEALPix pixels of a map with `nside = 8`.
+
+The function is easily extensible. For example, instead of using all HEALPix
+pixels, there might be a footprint that describes a specific survey:
+
+```yaml
+regions: !healpix
+  mask: survey-footprint-n512.fits
+```
+
+The `mask` keyword can be combined with the `nside` parameter to change the
+resolution of the mask if requested.
+
+If the HEALPix maps become finely resolved, it may be desirable to combine
+several pixels into a single region. There may be a `batch` [alt: `combine`]
+keyword for this purpose:
+
+```yaml
+regions: !healpix
+  nside: 8
+  batch: 16
+```
+
+The resulting list of regions will contain `768/16 = 48` regions. The `batch`
+keyword may also take a quantity of solid angle and automatically choose the
+number of pixels to combine accordingly.
+
+
+Map making
+----------
+
+This ADR does not address the problem of how maps will be generated from the
+list of regions. For example, a very real use case would be to generate
+populations of galaxies and simply count the total number in each HEALPix pixel
+to generate a density map. This will be addressed in a separate ADR.
+
+
+Consequences
+------------
+The existing `Pipeline` class must be extended to support iterating regions. No
+existing interfaces are affected.
diff --git a/docs/conf.py b/docs/conf.py
index 761d5f7..86922fb 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -165,7 +165,7 @@ if setup_cfg.get('edit_on_github').lower() == 'true':
     extensions += ['sphinx_astropy.ext.edit_on_github']
 
     edit_on_github_project = setup_cfg['github_project']
-    edit_on_github_branch = "master"
+    edit_on_github_branch = "main"
 
     edit_on_github_source_root = ""
     edit_on_github_doc_root = "docs"
diff --git a/docs/configuration_files.rst b/docs/configuration_files.rst
new file mode 100644
index 0000000..7e3d710
--- /dev/null
+++ b/docs/configuration_files.rst
@@ -0,0 +1,349 @@
+###################
+Configuration Files
+###################
+
+This page outlines how to construct configuration files to run your own routines
+with `~skypy.pipeline.Pipeline`.
+
+`SkyPy` is an astrophysical simulation pipeline tool that allows to define any
+arbitrary workflow and store data in table format. You may use `SkyPy` `~skypy.pipeline.Pipeline`
+to call any function --your own implementation, from any compatible external software or from the `SkyPy library`.
+Then `SkyPy` deals with the data dependencies and provides a library of functions to be used with it.
+
+These guidelines start with an example using one of the `SkyPy` functions, and it follows
+the concrete YAML syntax necessary for you to write your own configuration files, beyond using `SkyPy`
+functions.
+
+SkyPy example
+-------------
+
+In this section, we exemplify how you can write a configuration file and use some of the `SkyPy` functions.
+In this example, we sample redshifts and magnitudes from the SkyPy luminosity function, `~skypy.galaxies.schechter_lf`.
+
+- `Define variables`:
+
+From the documentation, the parameters for the `~skypy.galaxies.schechter_lf` function are: ``redshift``, the characteristic absolute magnitude ``M_star``, the amplitude ``phi_star``, faint-end slope parameter ``alpha``,
+the magnitude limit ``magnitude_limit``, the fraction of sky ``sky_area``, ``cosmology`` and ``noise``.
+If you are planning to reuse some of these parameters, you can define them at the top-level of your configuration file.
+In our example, we are using ``Astropy`` linear and exponential models for the characteristic absolute magnitude and the amplitude, respectively.
+Also, ``noise`` is an optional parameter and you could use its default value by omitting its definition.
+
+  .. code:: yaml
+
+    cosmology: !astropy.cosmology.default_cosmology.get []
+    z_range: !numpy.linspace [0, 2, 21]
+    M_star: !astropy.modeling.models.Linear1D [-0.9, -20.4]
+    phi_star: !astropy.modeling.models.Exponential1D [3e-3, -9.7]
+    magnitude_limit: 23
+    sky_area: 0.1 deg2
+
+- `Tables and columns`:
+
+You can create a table ``blue_galaxies`` and for now add the columns for redshift and magnitude (note here the ``schechter_lf`` returns a 2D object)
+
+  .. code:: yaml
+
+    tables:
+      blue_galaxies:
+        redshift, magnitude: !skypy.galaxies.schechter_lf
+          redshift: $z_range
+      	  M_star: $M_star
+      	  phi_star: $phi_star
+      	  alpha: -1.3
+      	  m_lim: $magnitude_limit
+      	  sky_area: $sky_area
+
+`Important:` if cosmology is detected as a parameter but is not set, it automatically uses the cosmology variable defined at the top-level of the file.
+
+This is how the entire configuration file looks like!
+
+.. literalinclude:: luminosity.yml
+   :language: yaml
+
+You may now save it as ``luminosity.yml`` and run it using the `SkyPy` `~skypy.pipeline.Pipeline`:
+
+.. plot::
+   :include-source: true
+   :context: close-figs
+
+    import matplotlib.pyplot as plt
+    from skypy.pipeline import Pipeline
+
+    # Execute SkyPy luminosity pipeline
+    pipeline = Pipeline.read("luminosity.yml")
+    pipeline.execute()
+
+    # Blue population
+    skypy_galaxies = pipeline['blue_galaxies']
+
+    # Plot histograms
+    fig, axs = plt.subplots(1, 2, figsize=(9, 3))
+
+    axs[0].hist(skypy_galaxies['redshift'], bins=50, histtype='step', color='purple')
+    axs[0].set_xlabel(r'$Redshift$')
+    axs[0].set_ylabel(r'$\mathrm{N}$')
+    axs[0].set_yscale('log')
+
+    axs[1].hist(skypy_galaxies['magnitude'], bins=50, histtype='step', color='green')
+    axs[1].set_xlabel(r'$Magnitude$')
+    axs[1].set_yscale('log')
+
+    plt.tight_layout()
+    plt.show()
+
+You can also run the pipeline directly from the command line and write the outputs to a fits file:
+
+.. code-block:: bash
+
+    $ skypy luminosity.yml luminosity.fits
+
+
+
+Don’t forget to check out for more complete examples_!
+
+.. _examples: https://skypy.readthedocs.io/en/stable/examples/index.html
+
+
+YAML syntax
+-----------
+YAML_ is a file format designed to be readable by both computers and humans.
+Fundamentally, a file written in YAML consists of a set of key-value pairs.
+Each pair is written as ``key: value``, where whitespace after the ``:`` is optional.
+The hash character ``#`` denotes the start of a comment and all further text on that
+line is ignored by the parser.
+
+
+This guide introduces the main syntax of YAML relevant when writing
+a configuration file to use with ``SkyPy``. Essentially, it begins with
+definitions of individual variables at the top level, followed by the tables,
+and, within the table entries, the features of objects to simulate are included.
+Main keywords: ``parameters``, ``cosmology``, ``tables``.
+
+
+Variables
+^^^^^^^^^
+* `Variable definition`: a variable is defined as a key-value pair at the top of the file.
+  YAML is able to interpret any numeric data with the appropriate type: integer, float, boolean.
+  Similarly for lists and dictionaries.
+  In addition, SkyPy has added extra functionality to interpret and store Astropy Quantities_.
+  Everything else is stored as a string (with or without explicitly using quotes)
+
+  .. code:: yaml
+
+      # YAML interprets
+      counter: 100 # An integer
+      miles: 1000.0 # A floating point
+      name: "Joy" # A string
+      planet: Earth # Another string
+      mylist: [ 'abc', 789, 2.0e3 ] # A list
+      mydict: { 'fruit': 'orange', 'year': 2020 } # A dictionary
+
+      # SkyPy extra functionality
+      angle: 10 deg
+      distance: 300 kpc
+
+
+* `Import objects`:
+  the SkyPy configuration syntax allows objects to be imported directly from external
+  (sub)modules using the ``!`` tag and providing neither a list of arguments or a
+  dictionary of keywords. For example, this enables the import and usage of any Astropy cosmology:
+
+  .. code:: yaml
+
+      cosmology: !astropy.cosmology.Planck13 # import the Planck13 object and bind it to the variable named "cosmology"
+
+
+Parameters
+^^^^^^^^^^
+
+* `Parameters definition`: parameters are variables that can be modified at execution.
+
+  For example,
+
+  .. code:: yaml
+
+      parameters:
+        hubble_constant: 70
+        omega_matter: 0.3
+
+
+Functions
+^^^^^^^^^
+* `Function call`: functions are defined as tuples where the first entry is the fully qualified function name tagged with and exclamation mark ``!`` and the second entry is either a list of positional arguments or a dictionary of keyword arguments.
+
+  For example, if you need to call the ``log10()`` and ``linspace()`` NumPy_ functions, then you define the following key-value pairs:
+
+  .. code:: yaml
+
+      log_of_2: !numpy.log10 [2]
+      myarray: !numpy.linspace [0, 2.5, 10]
+
+  You can also define parameters of functions with a dictionary of keyword arguments.
+  Imagine you want to compute the comoving distance for a range of redshifts and an `Astropy` Planck 2015 cosmology.
+  To run it with the `SkyPy` pipeline, call the function and define the parameters as an indented dictionary.
+
+  .. code:: yaml
+
+      comoving_distance: !astropy.cosmology.Planck15.comoving_distance
+        z: !numpy.linspace [ 0, 1.3, 10 ]
+
+  Similarly, you can specify the functions arguments as a dictionary:
+
+  .. code:: yaml
+
+      comoving_distance: !astropy.cosmology.Planck15.comoving_distance
+        z: !numpy.linspace {start: 0, stop: 1.3, num: 10}
+
+  `N.B.` To call a function with no arguments, you should pass an empty list of
+  ``args`` or an empty dictionary of ``kwargs``. For example:
+
+  .. code:: yaml
+
+      cosmo: !astropy.cosmology.default_cosmology.get []
+
+
+* `Variable reference`: variables can be referenced by their full name tagged with a dollar sign ``$``.
+  In the previous example you could also define the variables at the top-level of the file and then reference them:
+
+  .. code:: yaml
+
+      redshift: !numpy.linspace [ 0, 1.3, 10 ]
+      comoving_distance: !astropy.cosmology.Planck15.comoving_distance
+        z: $redshift
+
+* The `cosmology` to be used by functions within the pipeline only needs to be set up at the top. If a function needs ``cosmology`` as an input, you need not define it again, it is automatically detected.
+
+  For example, calculate the angular size of a galaxy with a given physical size, at a fixed redshift and for a given cosmology:
+
+  .. code:: yaml
+
+      cosmology: !astropy.cosmology.FlatLambdaCDM
+        H0: 70
+        Om0: 0.3
+      size: !skypy.galaxies.morphology.angular_size
+        physical_size: 10 kpc
+        redshift: 0.2
+
+* `Job completion`: ``.depends`` can be used to force any function call to wait for completion
+  of any other job.
+
+  A simple example where, for some reason, the comoving distance needs to be called after
+  completion of the angular size function:
+
+  .. code:: yaml
+
+    cosmology: !astropy.cosmology.Planck15
+    size: !skypy.galaxies.morphology.angular_size
+      physical_size: 10 kpc
+      redshift: 0.2
+    comoving_distance: !astropy.cosmology.Planck15.comoving_distance
+      z: !numpy.linspace [ 0, 1.3, 10 ]
+      .depends: size
+
+  By doing so, you force the function call ``redshift`` to be completed before is used to compute the comoving distance.
+
+
+Tables
+^^^^^^
+
+* `Table creation`: a dictionary of table names, each resolving to a dictionary of column names for that table.
+
+  Let us create a table called ``telescope`` with a column to store the width of spectral lines that follow a normal distribution
+
+  .. code:: yaml
+
+      tables:
+        telescope:
+          spectral_lines: !scipy.stats.norm.rvs
+            loc: 550
+            scale: 1.6
+            size: 100
+
+* `Column addition`: you can add as many columns to a table as you need.
+    Imagine you want to add a column for the telescope collecting surface
+
+  .. code:: yaml
+
+      tables:
+        telescope:
+          spectral_lines: !scipy.stats.norm.rvs
+            loc: 550
+            scale: 1.6
+            size: 100
+          collecting_surface: !numpy.random.uniform
+            low: 6.9
+            high: 7.1
+            size: 100
+
+* `Column reference`: columns in the pipeline can be referenced by their full name tagged with a dollar sign ``$``.
+  Example: the galaxy mass that follows a lognormal distribution. You can create a table ``galaxies``
+  with a column ``mass`` where you sample 10000 object and a second column, ``radius`` which also follows a lognormal distribution
+  but the mean depends on how massive the galaxies are:
+
+  .. code:: yaml
+
+    tables:
+      galaxies:
+        mass: !numpy.random.lognormal
+          mean: 5.
+          size: 10000
+        radius: !numpy.random.lognormal
+          mean: $galaxies.mass
+
+
+* `Multi-column assignment`: multi-column assignment is performed with any 2d-array, where one of the dimensions is interpreted
+  as the rows of the table and the second dimension, as separate columns. Or you can do it from a function that returns a tuple.
+
+  We use multi-column assignment in the following example where we sample a two-dimensional array of values from a lognormal distribution and then store them as three columns in a table:
+
+  .. code:: yaml
+
+    tables:
+      halos:
+        mass, radius, concentration: !numpy.random.lognormal
+          size: [10000, 3]
+
+
+* `Table initialisation`: by default tables are initialised using ``astropy.table.Table()`` however this can be overridden using the ``.init`` keyword to initialise the table with any function call.
+
+  For example, you can stack galaxy properties such as radii and mass:
+
+  .. code:: yaml
+
+    radii: !numpy.logspace [ 1, 2, 100 ]
+    mass: !numpy.logspace [ 9, 12, 100 ]
+    tables:
+      galaxies:
+        .init: !astropy.table.vstack [[ $radii, $mass ]]
+
+
+* `Table reference`: when a function call depends on tables, you need to ensure the referenced table has the necessary content and is not empty.
+  You can do that with ``.complete``.
+
+  Example: you want to perform a very simple abundance matching, i.e. painting galaxies within your halos.
+  You can create two tables ``halos`` and ``galaxies`` storing the halo mass and galaxy luminosities.
+  Then you can stack these two tables and store it in a third table called ``matching``.
+
+  .. code:: yaml
+
+    tables:
+      halos:
+        halo_mass: !numpy.random.uniform
+          low: 1.0e8
+          high: 1.0e14
+          size: 20
+      galaxies:
+        luminosity: !numpy.random.uniform
+          low: 0.05
+          high: 10.0
+          size: 20
+      matching:
+        .init: !astropy.table.hstack
+          tables: [ $halos, $galaxies ]
+          .depends: [ halos.complete, galaxies.complete ]
+
+
+.. _YAML: https://yaml.org
+.. _NumPy: https://numpy.org
+.. _Quantities: https://docs.astropy.org/en/stable/units/
+.. _clone(): https://docs.astropy.org/en/stable/api/astropy.cosmology.FLRW.html?highlight=clone#astropy.cosmology.FLRW.clone
diff --git a/docs/feature_list.rst b/docs/feature_list.rst
new file mode 100644
index 0000000..6df41db
--- /dev/null
+++ b/docs/feature_list.rst
@@ -0,0 +1,20 @@
+#############
+Feature List
+#############
+
+This page outlines the main features in the `SkyPy` library.
+
+Galaxies
+--------
+
+- `Luminosity Distributions`_
+- `Morphological Distributions`_
+- `Redshift Distributions`_
+- `Spectral Energy Distribution Modelling`_
+- `Stellar Mass Distributions`_
+
+.. _Luminosity Distributions: https://skypy.readthedocs.io/en/latest/galaxies.html#module-skypy.galaxies.luminosity
+.. _Morphological Distributions: https://skypy.readthedocs.io/en/latest/galaxies.html#module-skypy.galaxies.morphology
+.. _Redshift Distributions: https://skypy.readthedocs.io/en/latest/galaxies.html#module-skypy.galaxies.redshift
+.. _Spectral Energy Distribution Modelling: https://skypy.readthedocs.io/en/latest/galaxies.html#module-skypy.galaxies.spectrum
+.. _Stellar Mass Distributions: https://skypy.readthedocs.io/en/latest/galaxies.html#module-skypy.galaxies.stellar_mass
diff --git a/docs/galaxies.rst b/docs/galaxies.rst
index b805a21..bfddeee 100644
--- a/docs/galaxies.rst
+++ b/docs/galaxies.rst
@@ -63,6 +63,13 @@ Spectrum
 
 The following models are found in the `skypy.galaxies.spectrum` package.
 
+SkyPy uses the `speclite <https://speclite.readthedocs.io/>`_ package for
+photometric calculations. Some of the following functions take the names of
+photometric filters as an input parameter. Users can choose from the available
+`Speclite Filters <https://speclite.readthedocs.io/en/latest/filters.html>`_
+following the naming syntax described in `speclite.filters.load_filters`, or
+create their own named `speclite.filters.FilterResponse`.
+
 .. currentmodule:: skypy.galaxies.spectrum
 .. autosummary::
    :nosignatures:
diff --git a/docs/halos/index.rst b/docs/halos/index.rst
index d0d82ef..3aed5f3 100644
--- a/docs/halos/index.rst
+++ b/docs/halos/index.rst
@@ -22,7 +22,7 @@ a python script, for example.
     k = np.logspace(-3, 1, num=1000, base=10.0)
     mass = 10**np.arange(9.0, 15.0, 0.1)
 
-    pk0 = eisenstein_hu(k, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)
+    pk0 = eisenstein_hu(k, A_s, n_s, Planck15, kwmap=0.02, wiggle=True) # doctest: +SKIP
     sigma = np.sqrt(_sigma_squared(mass, k, pk0, growth_0, Planck15))
 
     # Collapse functions
@@ -59,36 +59,6 @@ a python script, for example.
     plt.show()
 
 
-You can also sample halos using their mass function. For this, you can use a config
-file and run the pipeline, for example.
-
-.. literalinclude:: examples/halos.yml
-   :language: yaml
-
-.. plot::
-   :include-source: false
-   :context: close-figs
-
-    from skypy.pipeline import Pipeline
-    pipeline = Pipeline.read('examples/halos.yml')
-    pipeline.execute()
-
-    # Draw from different halo mass samplers
-    halo_massST = pipeline['sheth-tormen']
-    halo_massPS = pipeline['press-schechter']
-
-    plt.hist(np.log10(halo_massST), histtype='step', label='Sheth-Tormen')
-    plt.hist(np.log10(halo_massPS), histtype='step', label='Press-Schechter')
-
-    # axis label and title
-    plt.xlabel(r'$log(mass)$')
-    plt.title('Halo sampler')
-
-    # show plot labels
-    plt.legend()
-    plt.show()
-
-
 Halos (`skypy.halos`)
 =====================
 
diff --git a/docs/index.rst b/docs/index.rst
index 11458ce..2d916e3 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -4,7 +4,15 @@
 SkyPy Documentation
 ###################
 
-This package contains methods for modelling the Universe, galaxies and the Milky Way. Also included are methods for generating observed data.
+This package contains methods for modelling the Universe, galaxies and the
+Milky Way. SkyPy simulates populations of astronomical objects, generating
+random realisations of intrinsic and observed properties, with the
+intention the simulations can then be compared to data as part of an inference
+pipeline.
+
+.. Important:: If you use SkyPy for work presented in a publication or talk
+   please follow our :doc:`project/citation`.
+
 
 
 .. _getting-started:
@@ -17,6 +25,8 @@ Getting Started
    :maxdepth: 1
 
    install
+   feature_list
+   configuration_files
    examples/index
 
 .. _user-docs:
@@ -67,6 +77,7 @@ Project details
    :maxdepth: 1
 
    project/code_of_conduct
+   project/citation
 
 
 *****
diff --git a/docs/install.rst b/docs/install.rst
index e52b8bb..986f8d9 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -29,18 +29,18 @@ channel. To install the latest version for your active conda environment:
 From GitHub
 -----------
 
-The latest development version of SkyPy can be found on the master branch of
+The latest development version of SkyPy can be found on the main branch of
 the `skypyproject/skypy`_ GitHub repository. This and any other branch or tag
 can be installed directly from GitHub using a recent version of pip:
 
 .. code:: console
 
-    $ pip install skypy@git+https://github.com/skypyproject/skypy.git@master
+    $ pip install skypy@git+https://github.com/skypyproject/skypy.git@main
 
 Dependencies
 ------------
 
-SkyPy is compatble with Python versions 3.6 or later on Ubuntu, macOS and
+SkyPy is compatble with Python versions 3.7 or later on Ubuntu, macOS and
 Windows operating systems. It has the following core dependencies:
 
 - `astropy <https://www.astropy.org/>`__
@@ -53,6 +53,7 @@ Installing using pip or conda will automatically install or update these core
 dependencies if necessary. SkyPy also has a number of optional dependencies
 that enable additional features:
 
+- `h5py <https://www.h5py.org/>`_
 - `speclite <https://speclite.readthedocs.io/>`_
 
 To install SkyPy with all optional dependencies using pip:
diff --git a/docs/luminosity.yml b/docs/luminosity.yml
new file mode 100644
index 0000000..0ca5ed0
--- /dev/null
+++ b/docs/luminosity.yml
@@ -0,0 +1,15 @@
+cosmology: !astropy.cosmology.default_cosmology.get []
+z_range: !numpy.linspace [0, 2, 21]
+M_star: !astropy.modeling.models.Linear1D [-0.9, -20.4]
+phi_star: !astropy.modeling.models.Exponential1D [3e-3, -9.7]
+magnitude_limit: 23
+sky_area: 0.1 deg2
+tables:
+  blue_galaxies:
+    redshift, magnitude: !skypy.galaxies.schechter_lf
+      redshift: $z_range
+      M_star: $M_star
+      phi_star: $phi_star
+      alpha: -1.3
+      m_lim: $magnitude_limit
+      sky_area: $sky_area
diff --git a/docs/pipeline/index.rst b/docs/pipeline/index.rst
index 51c65df..adb1d8e 100644
--- a/docs/pipeline/index.rst
+++ b/docs/pipeline/index.rst
@@ -19,7 +19,10 @@ fits files:
 
 .. code-block:: bash
 
-    $ skypy examples/galaxies/sdss_photometry.yml --format fits
+    $ skypy examples/galaxies/sdss_photometry.yml sdss_photometry.fits
+
+To view the progress of the pipeline as it runs you can enable logging using the
+`--verbose` flag.
 
 Config files are written in YAML format and read using the
 `~skypy.pipeline.load_skypy_yaml` funciton. Each entry in the config specifices
diff --git a/docs/power_spectrum/index.rst b/docs/power_spectrum/index.rst
index a138981..1f00ee5 100644
--- a/docs/power_spectrum/index.rst
+++ b/docs/power_spectrum/index.rst
@@ -9,29 +9,7 @@ matter power spectrum, including `~skypy.power_spectrum.camb` and
 matter power spectrum using `~skypy.power_spectrum.eisenstein_hu` and the
 non-linear corrections using `~skypy.power_spectrum.halofit_smith`:
 
-.. literalinclude:: examples/power_spectrum.yml
-   :language: yaml
 
-.. plot::
-   :include-source: false
-   :context: close-figs
-
-    import matplotlib.pyplot as plt
-    from skypy.pipeline import Pipeline
-
-    pipeline = Pipeline.read('examples/power_spectrum.yml')
-    pipeline.execute()
-
-    # Eisenstein and Hu power spectrum and Halofit matter power spectra
-    k = pipeline['wavenumber']
-    power_EH_w = pipeline['eisenstein_hu_wiggle']
-    hf_Smith = pipeline['halofit']
-
-    plt.loglog(k, power_EH_w, label='Eisenstein & Hu')
-    plt.loglog(k, hf_Smith, '--', label='Halofit')
-    plt.xlabel(r'Wavenumber $(1/Mpc)$')
-    plt.ylabel(r'Power spectrum $(Mpc^3)$')
-    plt.legend(frameon=False, loc='lower left');
 
 
 
diff --git a/docs/project/citation.rst b/docs/project/citation.rst
new file mode 100644
index 0000000..ff43a05
--- /dev/null
+++ b/docs/project/citation.rst
@@ -0,0 +1,1 @@
+.. include:: ../../CITATION.rst
diff --git a/docs/utils/index.rst b/docs/utils/index.rst
index a99ca10..090453e 100644
--- a/docs/utils/index.rst
+++ b/docs/utils/index.rst
@@ -126,6 +126,14 @@ Photometry (`skypy.utils.photometry`)
 This module contains methods that model spectral energy distributions and
 calculate photometric properties.
 
+SkyPy uses the `speclite <https://speclite.readthedocs.io/>`_ package for
+photometric calculations. Some of the following functions take the names of
+photometric filters as an input parameter. Users can choose from the available
+`Speclite Filters <https://speclite.readthedocs.io/en/latest/filters.html>`_
+following the naming syntax described in `speclite.filters.load_filters`, or
+create their own named `speclite.filters.FilterResponse`.
+
+
 .. currentmodule:: skypy.utils.photometry
 .. autosummary::
    :nosignatures:
@@ -135,6 +143,8 @@ calculate photometric properties.
    luminosity_in_band
    mag_ab
    SpectrumTemplates
+   magnitude_error_rykoff
+   logistic_completeness_function
 
 
 Random sampling (`skypy.utils.random`)
diff --git a/examples/galaxies/plot_photometry.py b/examples/galaxies/plot_photometry.py
index 9df4b46..c41988d 100644
--- a/examples/galaxies/plot_photometry.py
+++ b/examples/galaxies/plot_photometry.py
@@ -40,7 +40,7 @@ kcorrect spectral energy distribution templates as implemented in SkyPy.
 #
 # .. code-block:: bash
 #
-#    $ skypy examples/galaxies/sdss_photometry.yml --format fits
+#    $ skypy examples/galaxies/sdss_photometry.yml sdss_photometry.fits
 #
 # or in a python script using the :class:`Pipeline <skypy.pipeline.Pipeline>`
 # class as demonstrated in the `SDSS Photometry`_ section below. For more
diff --git a/examples/galaxies/plot_schechter.py b/examples/galaxies/plot_schechter.py
index 48e14f5..11b060b 100644
--- a/examples/galaxies/plot_schechter.py
+++ b/examples/galaxies/plot_schechter.py
@@ -100,7 +100,7 @@ for ax, (z_min, z_max) in zip([a1, a2, a3, a4], z_slices):
     ax.plot(bins, phi_model, label='Model', color='b')
     ax.text(-14.3, 5e-7, r'${:.2f} \leq z < {:.2f}$'.format(z_min, z_max))
     ax.set_xlabel(r'$M_B$')
-    ax.set_ylabel(r'$\log_{10} \Phi \, [\mathrm{Mpc}^{-3} \, \mathrm{mag}^{-1}]$')
+    ax.set_ylabel(r'$\Phi \, [\mathrm{Mpc}^{-3} \, \mathrm{mag}^{-1}]$')
     ax.set_yscale('log')
     ax.set_xlim([-14, -24])
     ax.set_ylim([3e-7, 3e-2])
diff --git a/examples/galaxies/sdss_dered_10deg2.ecsv b/examples/galaxies/sdss_dered_10deg2.ecsv
index f827d3c..5ae82dd 100644
--- a/examples/galaxies/sdss_dered_10deg2.ecsv
+++ b/examples/galaxies/sdss_dered_10deg2.ecsv
@@ -9,15 +9,15 @@
 # - name: magnitude
 #   datatype: float32
 # - name: dered_u
-#   datatype: int
+#   datatype: uint16
 # - name: dered_g
-#   datatype: int
+#   datatype: uint16
 # - name: dered_r
-#   datatype: int
+#   datatype: uint16
 # - name: dered_i
-#   datatype: int
+#   datatype: uint16
 # - name: dered_z
-#   datatype: int
+#   datatype: uint16
 magnitude dered_u dered_g dered_r dered_i dered_z
 15.0 5 6 5 9 9
 15.1 4 5 11 8 17
diff --git a/examples/galaxies/sdss_photometry.yml b/examples/galaxies/sdss_photometry.yml
index 5d370bf..651e2e5 100644
--- a/examples/galaxies/sdss_photometry.yml
+++ b/examples/galaxies/sdss_photometry.yml
@@ -1,4 +1,4 @@
-cosmology: !astropy.cosmology.default_cosmology.get
+cosmology: !astropy.cosmology.default_cosmology.get []
 z_range: !numpy.linspace [0, 2, 21]
 M_star: !astropy.modeling.models.Linear1D [-0.9408582, -20.40492365]
 phi_star: !astropy.modeling.models.Exponential1D [0.00370253, -9.73858]
diff --git a/setup.cfg b/setup.cfg
index 44eb052..467facf 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -14,7 +14,7 @@ github_project = skypyproject/skypy
 [options]
 zip_safe = False
 packages = find:
-python_requires = >=3.6
+python_requires = >=3.7
 setup_requires = setuptools_scm
 install_requires =
     astropy>=4
@@ -31,15 +31,19 @@ console_scripts =
 test =
     pytest-astropy
     pytest-rerunfailures
-    speclite>=0.11
+    speclite>=0.14
 all =
-    speclite>=0.11
+    h5py
+    speclite>=0.14
+    six
     colossus
 docs =
     sphinx-astropy
+    graphviz
     matplotlib
+    six
     colossus
-    speclite>=0.11
+    speclite>=0.14
 
 [options.package_data]
 skypy = data/*,data/*/*,*/tests/data/*
diff --git a/skypy/galaxies/__init__.py b/skypy/galaxies/__init__.py
index 12630ea..f6307c3 100644
--- a/skypy/galaxies/__init__.py
+++ b/skypy/galaxies/__init__.py
@@ -5,6 +5,7 @@ populations.
 
 __all__ = [
     'schechter_lf',
+    'schechter_smf',
 ]
 
 from . import luminosity  # noqa F401,F403
diff --git a/skypy/galaxies/luminosity.py b/skypy/galaxies/luminosity.py
index 97717ef..57ce1ec 100644
--- a/skypy/galaxies/luminosity.py
+++ b/skypy/galaxies/luminosity.py
@@ -46,18 +46,6 @@ def schechter_lf_magnitude(redshift, M_star, alpha, m_lim, cosmology, size=None,
         Absolute magnitude sampled from a Schechter luminosity function for
         each input galaxy redshift.
 
-    Examples
-    --------
-
-    Sample a number of blue (alpha = -1.3, M_star = -20.5) galaxy magnitudes
-    brighter than m = 22.0 around redshift 0.5.
-
-    >>> import numpy as np
-    >>> from skypy.galaxies.luminosity import schechter_lf_magnitude
-    >>> from astropy.cosmology import Planck15
-    >>> z = np.random.uniform(4.9, 5.1, size=20)
-    >>> M = schechter_lf_magnitude(z, -20.5, -1.3, 22.0, Planck15)
-
     '''
 
     # only alpha scalars supported at the moment
diff --git a/skypy/galaxies/morphology.py b/skypy/galaxies/morphology.py
index 5272c48..debfc8d 100644
--- a/skypy/galaxies/morphology.py
+++ b/skypy/galaxies/morphology.py
@@ -40,13 +40,6 @@ def angular_size(physical_size, redshift, cosmology):
     angular_size : astropy.Quantity
         Angular distances in units of [rad] for a given radius.
 
-    Examples
-    --------
-    >>> from astropy import units
-    >>> from skypy.galaxies.morphology import angular_size
-    >>> from astropy.cosmology import Planck15
-    >>> r = angular_size(10*units.kpc, 1, Planck15)
-
     References
     ----------
     .. [1] D. W. Hogg, (1999), astro-ph/9905116.
@@ -89,14 +82,6 @@ def beta_ellipticity(e_ratio, e_sum, size=None):
     ----------
     .. [1] Kacprzak T., Herbel J., Nicola A. et al., arXiv:1906.01018
 
-    Examples
-    --------
-    Sample 10000 random variates from the Kacprzak model with
-    :math:`e_{\rm ratio} = 0.5` and :math:`e_{\rm sum} = 1.0`:
-
-    >>> from skypy.galaxies.morphology import beta_ellipticity
-    >>> ellipticity = beta_ellipticity(0.5, 1.0, size=10000)
-
     '''
 
     # convert to beta distribution parameters
@@ -139,17 +124,6 @@ def late_type_lognormal_size(magnitude, alpha, beta, gamma, M0, sigma1, sigma2,
         shape is (ns,). If magnitude has shape (nm,) and size=None,
         shape is (nm,).
 
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from skypy.galaxies import morphology
-    >>> magnitude = -16.0
-    >>> alpha, beta, gamma, M0 = 0.21, 0.53, -1.31, -20.52
-    >>> sigma1, sigma2 = 0.48, 0.25
-    >>> s = morphology.late_type_lognormal_size(magnitude, alpha, beta, gamma,
-    ...                                         M0, sigma1, sigma2)
-
-
     References
     ----------
     .. [1] S. Shen, H.J. Mo, S.D.M. White, M.R. Blanton, G. Kauffmann, W. Voges,
@@ -200,17 +174,6 @@ def early_type_lognormal_size(magnitude, a, b, M0, sigma1, sigma2, size=None):
         shape is (ns,). If magnitude has shape (nm,) and size=None,
         shape is (nm,).
 
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from skypy.galaxies import morphology
-    >>> magnitude = -20.0
-    >>> a, b, M0 = 0.6, -4.63, -20.52
-    >>> sigma1, sigma2 = 0.48, 0.25
-    >>> s = morphology.early_type_lognormal_size(magnitude, a, b, M0, sigma1,
-    ...                                          sigma2)
-
-
     References
     ----------
     .. [1] S. Shen, H.J. Mo, S.D.M. White, M.R. Blanton, G. Kauffmann, W. Voges,
@@ -253,14 +216,6 @@ def linear_lognormal_size(magnitude, a_mu, b_mu, sigma, size=None):
         shape is (ns,). If magnitude has shape (nm,) and size=None,
         shape is (nm,).
 
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from skypy.galaxies import morphology
-    >>> magnitude = -20.0
-    >>> a_mu, b_mu, sigma =-0.24, -4.63, 0.4
-    >>> s = morphology.linear_lognormal_size(magnitude, a_mu, b_mu, sigma)
-
     References
     ----------
     .. [1] J. Herbel, T. Kacprzak, A. Amara, A. Refregier, C.Bruderer and
@@ -296,19 +251,15 @@ def ryden04_ellipticity(mu_gamma, sigma_gamma, mu, sigma, size=None):
     size : int or tuple of ints or None
         Size of the sample. If `None` the size is inferred from the parameters.
 
+    Returns
+    -------
+    ellipticity: (size,) array_like
+        Ellipticities sampled from the Ryden 2004 model.
+
     References
     ----------
     .. [1] Ryden B. S., 2004, ApJ, 601, 214
 
-    Examples
-    --------
-    Sample 10000 random variates from the Ryden (2004) model with parameters
-    :math:`\mu_\gamma = 0.222`, :math:`\sigma_\gamma = 0.056`,
-    :math:`\mu = -1.85`, and :math:`\sigma = 0.89`.
-
-    >>> from skypy.galaxies.morphology import ryden04_ellipticity
-    >>> ellipticity = ryden04_ellipticity(0.222, 0.056, -1.85, 0.89, size=10000)
-
     '''
 
     # get size if not given
diff --git a/skypy/galaxies/redshift.py b/skypy/galaxies/redshift.py
index fbe672a..f048a35 100644
--- a/skypy/galaxies/redshift.py
+++ b/skypy/galaxies/redshift.py
@@ -58,14 +58,6 @@ def smail(z_median, alpha, beta, size=None):
     .. [1] Smail I., Ellis R. S., Fitchett M. J., 1994, MNRAS, 270, 245
     .. [2] Amara A., Refregier A., 2007, MNRAS, 381, 1018
 
-    Examples
-    --------
-    Sample 10 random variates from the Smail model with `alpha = 1.5` and
-    `beta = 2` and median redshift `z_median = 1.2`.
-
-    >>> from skypy.galaxies.redshift import smail
-    >>> redshift = smail(1.2, 1.5, 2.0, size=10)
-
     '''
 
     k = (alpha+1)/beta
@@ -117,23 +109,6 @@ def schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, sky_area,
         Redshifts of the galaxy sample described by the Schechter luminosity
         function.
 
-    Examples
-    --------
-    Compute the number density of galaxies with redshifts between 0 and 5
-    for typical values of the "blue" galaxy luminosity function above an
-    apparent magnitude cut of 22 for a survey of 1 square degree = 1/41253 of
-    the sky.
-
-    >>> from skypy.galaxies.redshift import schechter_lf_redshift
-    >>> from astropy import units
-    >>> from astropy.cosmology import Planck15
-    >>> z = [0., 5.]
-    >>> M_star = -20.5
-    >>> phi_star = 3.5e-3
-    >>> alpha = -1.3
-    >>> sky_area = 1*units.deg**2
-    >>> z_gal = schechter_lf_redshift(z, M_star, phi_star, alpha, 22, sky_area, Planck15)
-
     '''
 
     # compute lower truncation of scaled Schechter random variable
@@ -260,17 +235,12 @@ def redshifts_from_comoving_density(redshift, density, sky_area, cosmology, nois
         Sampled redshifts such that the comoving number density of galaxies
         corresponds to the input distribution.
 
-    Examples
+    Warnings
     --------
-    Sample redshifts with a constant comoving number density 1e-3/Mpc3 up to
-    redshift 1 for a survey of 1 square degree = 1/41253 of the sky.
-
-    >>> from skypy.galaxies.redshift import redshifts_from_comoving_density
-    >>> from astropy import units
-    >>> from astropy.cosmology import Planck15
-    >>> z_range = np.arange(0, 1.01, 0.1)
-    >>> sky_area = 1*units.deg**2
-    >>> z_gal = redshifts_from_comoving_density(z_range, 1e-3, sky_area, Planck15)
+    The inverse cumulative distribution function is approximated from the
+    number density and comoving volume calculated at the given `redshift`
+    values. The user must choose suitable `redshift` values to satisfy their
+    desired numerical accuracy.
 
     '''
 
diff --git a/skypy/galaxies/spectrum.py b/skypy/galaxies/spectrum.py
index 0095d24..f549276 100644
--- a/skypy/galaxies/spectrum.py
+++ b/skypy/galaxies/spectrum.py
@@ -73,20 +73,6 @@ def dirichlet_coefficients(redshift, alpha0, alpha1, z1=1., weight=None):
     .. [2] Blanton M. R., Roweis S., 2007, The Astronomical Journal,
            Volume 133, Page 734
 
-    Examples
-    --------
-    >>> from skypy.galaxies.spectrum import dirichlet_coefficients
-    >>> import numpy as np
-
-    Sample the coefficients according to [1]_ for n blue galaxies with
-    redshifts between 0 and 1.
-
-    >>> n = 100000
-    >>> alpha0 = np.array([2.079, 3.524, 1.917, 1.992, 2.536])
-    >>> alpha1 = np.array([2.265, 3.862, 1.921, 1.685, 2.480])
-    >>> redshift = np.linspace(0,2, n)
-    >>> coefficients = dirichlet_coefficients(redshift, alpha0, alpha1)
-
     """
 
     if np.ndim(alpha0) != 1 or np.ndim(alpha1) != 1:
@@ -157,6 +143,28 @@ class KCorrectTemplates(SpectrumTemplates):
         Mt = self.absolute_magnitudes(coefficients, filter)
         return np.power(10, 0.4*(Mt-magnitudes))
 
+    def stellar_mass_remain(self, coefficients, magnitudes, filter):
+        r'''Compute total surviving stellar mass from absolute magnitudes.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of template coefficients.
+        magnitudes : (ng,) array_like
+            The magnitudes to match in the reference bandpass.
+        filter : str
+            A single reference bandpass filter specification for
+            `~speclite.filters.load_filters`.
+
+        Returns
+        -------
+        stellar_mass : (ng,) array_like
+            Total surviving stellar mass of each galaxy in template units.
+        '''
+        m = self.stellar_mass(coefficients, magnitudes, filter)
+        m *= np.dot(coefficients, self.mremain)
+        return m
+
     def metallicity(self, coefficients):
         r'''Galaxy metallicities from kcorrect templates.
 
diff --git a/skypy/galaxies/stellar_mass.py b/skypy/galaxies/stellar_mass.py
index 2b1a33e..0fd3cd8 100644
--- a/skypy/galaxies/stellar_mass.py
+++ b/skypy/galaxies/stellar_mass.py
@@ -4,6 +4,7 @@
 import numpy as np
 
 from ..utils.random import schechter
+from ..utils import dependent_argument
 
 
 __all__ = [
@@ -11,6 +12,8 @@ __all__ = [
 ]
 
 
+@dependent_argument('m_star', 'redshift')
+@dependent_argument('alpha', 'redshift')
 def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
                        resolution=1000):
     r""" Stellar masses following the Schechter mass function [1]_.
@@ -19,9 +22,10 @@ def schechter_smf_mass(redshift, alpha, m_star, m_min, m_max, size=None,
     ----------
     redshift : array_like
         Galaxy redshifts for which to sample magnitudes.
-    alpha : float
-        The alpha parameter in the Schechter stellar mass function.
-    m_star : (nm,) array-like
+    alpha : float or function
+        The alpha parameter in the Schechter stellar mass function. If function,
+        it must return a scalar value.
+    m_star : (nm,) array-like or function
         Characteristic stellar mass m_*.
     size: int, optional
          Output shape of stellar mass samples. If size is None and m_star
diff --git a/skypy/pipeline/__init__.py b/skypy/pipeline/__init__.py
index 48f02f8..9adc050 100644
--- a/skypy/pipeline/__init__.py
+++ b/skypy/pipeline/__init__.py
@@ -3,5 +3,9 @@ This module provides methods to pipeline together multiple models with
 dependencies and handle their outputs.
 """
 
+import logging
+
+log = logging.getLogger(__name__)
+
 from ._config import *  # noqa
 from ._pipeline import *  # noqa
diff --git a/skypy/pipeline/_config.py b/skypy/pipeline/_config.py
index 2b839d4..656402b 100644
--- a/skypy/pipeline/_config.py
+++ b/skypy/pipeline/_config.py
@@ -52,23 +52,23 @@ class SkyPyLoader(yaml.SafeLoader):
         return Ref(ref)
 
     def construct_call(self, name, node):
-        if isinstance(node, yaml.ScalarNode):
-            arg = self.construct_scalar(node)
-            args = [arg] if arg != '' else []
-            kwargs = {}
-        elif isinstance(node, yaml.SequenceNode):
-            args = self.construct_sequence(node)
-            kwargs = {}
-        elif isinstance(node, yaml.MappingNode):
-            args = []
-            kwargs = self.construct_mapping(node)
-
         try:
-            function = import_function(name)
+            object = import_function(name)
         except (ModuleNotFoundError, AttributeError) as e:
             raise ImportError(f'{e}\n{node.start_mark}') from e
 
-        return Call(function, args, kwargs)
+        if isinstance(node, yaml.ScalarNode):
+            if node.value:
+                raise ValueError(f'{node.value}: ScalarNode should be empty to import an object')
+            return object
+        else:
+            if isinstance(node, yaml.SequenceNode):
+                args = self.construct_sequence(node)
+                kwargs = {}
+            if isinstance(node, yaml.MappingNode):
+                args = []
+                kwargs = self.construct_mapping(node)
+            return Call(object, args, kwargs)
 
     def construct_quantity(self, node):
         value = self.construct_scalar(node)
diff --git a/skypy/pipeline/_items.py b/skypy/pipeline/_items.py
index 80481db..9886d22 100644
--- a/skypy/pipeline/_items.py
+++ b/skypy/pipeline/_items.py
@@ -1,6 +1,7 @@
 '''item types in the pipeline'''
 
 from collections.abc import Sequence, Mapping
+from . import log
 import inspect
 
 
@@ -88,4 +89,5 @@ class Call(Item):
         '''execute the call in the given pipeline'''
         args = pipeline.evaluate(self.args)
         kwargs = pipeline.evaluate(self.kwargs)
+        log.info(f"Calling {self.function.__name__}")
         return self.function(*args, **kwargs)
diff --git a/skypy/pipeline/_pipeline.py b/skypy/pipeline/_pipeline.py
index 1ba6641..84be0a6 100644
--- a/skypy/pipeline/_pipeline.py
+++ b/skypy/pipeline/_pipeline.py
@@ -9,7 +9,9 @@ from copy import copy, deepcopy
 from collections.abc import Sequence, Mapping
 from ._config import load_skypy_yaml
 from ._items import Item, Call, Ref
+from . import log
 import networkx
+import pathlib
 
 
 __all__ = [
@@ -65,7 +67,7 @@ class Pipeline:
 
         References
         ----------
-        .. [1] https://github.com/skypyproject/skypy/tree/master/examples
+        .. [1] https://github.com/skypyproject/skypy/tree/main/examples
 
         '''
 
@@ -100,8 +102,6 @@ class Pipeline:
             self.dag.add_node(job, skip=False)
             if isinstance(settings, Item):
                 items[job] = settings
-                # infer additional item properties from context
-                settings.infer(context)
         for table, columns in self.table_config.items():
             table_complete = '.'.join((table, 'complete'))
             self.dag.add_node(table_complete)
@@ -113,8 +113,6 @@ class Pipeline:
                 self.dag.add_edge(job, table_complete)
                 if isinstance(settings, Item):
                     items[job] = settings
-                    # infer additional item properties from context
-                    settings.infer(context)
                 # DAG nodes for individual columns in multi-column assignment
                 names = [n.strip() for n in column.split(',')]
                 if len(names) > 1:
@@ -136,6 +134,8 @@ class Pipeline:
                 while c:
                     self.dag.add_edge(c, d)
                     c, d = c.rpartition('.')[0], c
+            # infer additional item properties from context
+            settings.infer(context)
 
     def execute(self, parameters={}):
         r'''Run a pipeline.
@@ -163,6 +163,7 @@ class Pipeline:
 
         # Initialise cosmology from config parameters
         if self.cosmology is not None:
+            log.info("Setting cosmology")
             self.state['cosmology'] = self.evaluate(self.cosmology)
 
         # go through the jobs in dependency order
@@ -171,7 +172,8 @@ class Pipeline:
             skip = node.get('skip', True)
             if skip:
                 continue
-            elif job in self.config:
+            log.info(f"Generating {job}")
+            if job in self.config:
                 settings = self.config.get(job)
                 self.state[job] = self.evaluate(settings)
             else:
@@ -186,27 +188,65 @@ class Pipeline:
                     # Single column assignment
                     self.state[table][column] = self.evaluate(settings)
 
-    def write(self, file_format=None, overwrite=False):
+    def write(self, filename, overwrite=False):
         r'''Write pipeline results to disk.
 
         Parameters
         ----------
-        file_format : str
-            File format used to write tables. Files are written using the
-            Astropy unified file read/write interface; see [1]_ for supported
-            file formats. If None (default) tables are not written to file.
+        filename : str
+            Name of output file to be written. It must have one of the
+            supported file extensions for FITS (.fit .fits .fts) or HDF5
+            (.hdf5 .hd5 .he5 .h5).
         overwrite : bool
-            Whether to overwrite any existing files without warning.
+            If filename already exists, this flag indicates whether or not to
+            overwrite it (without warning).
+        '''
 
-        References
+        suffix = pathlib.Path(filename).suffix.lower()
+        _fits_suffixes = ('.fit', '.fits', '.fts')
+        _hdf5_suffixes = ('.hdf5', '.hd5', '.he5', '.h5')
+
+        if suffix in _fits_suffixes:
+            self.write_fits(filename, overwrite)
+        elif suffix in _hdf5_suffixes:
+            self.write_hdf5(filename, overwrite)
+        else:
+            raise ValueError(f'{suffix} is an unsupported file format. SkyPy supports '
+                             'FITS (' + ' '.join(_fits_suffixes) + ') and '
+                             'HDF5 (' + ' '.join(_hdf5_suffixes) + ').')
+
+    def write_fits(self, filename, overwrite=False):
+        r'''Write pipeline results to a FITS file.
+
+        Parameters
         ----------
-        .. [1] https://docs.astropy.org/en/stable/io/unified.html
+        filename : str
+            Name of output file to be written.
+        overwrite : bool
+            If filename already exists, this flag indicates whether or not to
+            overwrite it (without warning).
+        '''
+        from astropy.io.fits import HDUList, PrimaryHDU, table_to_hdu
+        hdul = [PrimaryHDU()]
+        for t in self.table_config:
+            hdu = table_to_hdu(self[t])
+            hdu.header['EXTNAME'] = t
+            hdul.append(hdu)
+        HDUList(hdul).writeto(filename, overwrite=overwrite)
+
+    def write_hdf5(self, filename, overwrite=False):
+        r'''Write pipeline results to a HDF5 file.
 
+        Parameters
+        ----------
+        filename : str
+            Name of output file to be written.
+        overwrite : bool
+            If filename already exists, this flag indicates whether or not to
+            overwrite it (without warning).
         '''
-        if file_format:
-            for table in self.table_config.keys():
-                filename = '.'.join((table, file_format))
-                self.state[table].write(filename, overwrite=overwrite)
+        for t in self.table_config:
+            self[t].write(filename, path=f'tables/{t}', append=True, overwrite=overwrite)
 
     def evaluate(self, value):
         '''evaluate an item in the pipeline'''
diff --git a/skypy/pipeline/scripts/skypy.py b/skypy/pipeline/scripts/skypy.py
index 0c6e825..6d723bc 100644
--- a/skypy/pipeline/scripts/skypy.py
+++ b/skypy/pipeline/scripts/skypy.py
@@ -1,6 +1,7 @@
 '''skypy Command Line Script'''
 
 import argparse
+import logging
 from skypy import __version__ as skypy_version
 from skypy.pipeline import Pipeline, load_skypy_yaml
 import sys
@@ -11,19 +12,40 @@ def main(args=None):
     parser = argparse.ArgumentParser(description="SkyPy pipeline driver")
     parser.add_argument('--version', action='version', version=skypy_version)
     parser.add_argument('config', help='Config file name')
-    parser.add_argument('-f', '--format', required=False,
-                        choices=['fits', 'hdf5'], help='Table file format')
+    parser.add_argument('output', help='Output file name')
     parser.add_argument('-o', '--overwrite', action='store_true',
                         help='Whether to overwrite existing files')
+    parser.add_argument("-v", "--verbose", action="count", default=0,
+                        help="Increase logging verbosity")
+    parser.add_argument("-q", "--quiet", action="count", default=0,
+                        help="Decrease logging verbosity")
 
     # get system args if none passed
     if args is None:
         args = sys.argv[1:]
 
     args = parser.parse_args(args or ['--help'])
-    config = load_skypy_yaml(args.config)
 
-    pipeline = Pipeline(config)
-    pipeline.execute()
-    pipeline.write(file_format=args.format, overwrite=args.overwrite)
-    return(0)
+    # Setup skypy logger
+    default_level = logging._nameToLevel['WARNING']
+    logging_level = default_level + 10 * (args.quiet - args.verbose)
+    formatter = logging.Formatter('%(asctime)s [%(levelname)s] %(name)s: %(message)s')
+    stream_handler = logging.StreamHandler()
+    stream_handler.setLevel(logging_level)
+    stream_handler.setFormatter(formatter)
+    logger = logging.getLogger('skypy')
+    logger.setLevel(logging_level)
+    logger.addHandler(stream_handler)
+
+    try:
+        config = load_skypy_yaml(args.config)
+        pipeline = Pipeline(config)
+        pipeline.execute()
+        if args.output:
+            logger.info(f"Writing {args.output}")
+            pipeline.write(args.output, overwrite=args.overwrite)
+    except Exception as e:
+        logger.exception(e)
+        raise SystemExit(2) from e
+
+    return 0
diff --git a/skypy/power_spectrum/_eisenstein_hu.py b/skypy/power_spectrum/_eisenstein_hu.py
index a1c8cea..136a7b5 100644
--- a/skypy/power_spectrum/_eisenstein_hu.py
+++ b/skypy/power_spectrum/_eisenstein_hu.py
@@ -53,7 +53,7 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
     >>> from astropy.cosmology import Planck15
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
-    >>> transfer_with_wiggles(wavenumber, A_s, n_s, Planck15, kwmap=0.02)
+    >>> transfer_with_wiggles(wavenumber, A_s, n_s, Planck15, kwmap=0.02)  # doctest: +SKIP
     array([9.92144790e-01, 7.78548704e-01, 1.29998169e-01, 4.63863054e-03,
        8.87918075e-05])
 
@@ -189,7 +189,7 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
     >>> from astropy.cosmology import Planck15
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
-    >>> transfer_no_wiggles(wavenumber, A_s, n_s, Planck15)
+    >>> transfer_no_wiggles(wavenumber, A_s, n_s, Planck15)  # doctest: +SKIP
     array([9.91959695e-01, 7.84518347e-01, 1.32327555e-01, 4.60773671e-03,
        8.78447096e-05])
 
@@ -271,7 +271,7 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
     >>> from astropy.cosmology import Planck15
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
-    >>> eisenstein_hu(wavenumber, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)
+    >>> eisenstein_hu(wavenumber, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)  # doctest: +SKIP
     array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03, 1.14604456e+02,
        3.91399918e-01])
 
diff --git a/skypy/power_spectrum/_growth.py b/skypy/power_spectrum/_growth.py
index 22ca95f..8ee7bd1 100644
--- a/skypy/power_spectrum/_growth.py
+++ b/skypy/power_spectrum/_growth.py
@@ -98,7 +98,7 @@ def growth_factor(redshift, cosmology, gamma=6.0/11.0):
     >>> import numpy as np
     >>> from astropy.cosmology import FlatLambdaCDM
     >>> cosmology = FlatLambdaCDM(H0=67.04, Om0=0.3183, Ob0=0.047745)
-    >>> growth_factor(0, cosmology)
+    >>> growth_factor(0, cosmology)  # doctest: +SKIP
     0.5355746155304598
 
     References
@@ -147,7 +147,7 @@ def growth_function(redshift, cosmology, gamma=6.0/11.0, z_upper=1100):
     >>> from scipy import integrate
     >>> from astropy.cosmology import FlatLambdaCDM
     >>> cosmology = FlatLambdaCDM(H0=67.04, Om0=0.3183, Ob0=0.047745)
-    >>> growth_function(0, cosmology)
+    >>> growth_function(0, cosmology)  # doctest: +SKIP
     0.7909271056297236
 
     References
diff --git a/skypy/power_spectrum/_halofit.py b/skypy/power_spectrum/_halofit.py
index 3c2d394..f91fcb2 100644
--- a/skypy/power_spectrum/_halofit.py
+++ b/skypy/power_spectrum/_halofit.py
@@ -104,8 +104,8 @@ def halofit(wavenumber, redshift, linear_power_spectrum,
     >>> z, A_s, n_s = 0, 2.2e-09, 0.97
     >>> cosmology = default_cosmology.get()
     >>> dz = growth_function(z, cosmology)
-    >>> linear_power = eisenstein_hu(k, A_s, n_s, cosmology) * np.square(dz)
-    >>> nonlinear_power = halofit_smith(k, z, linear_power, cosmology)
+    >>> linear_power = eisenstein_hu(k, A_s, n_s, cosmology) * np.square(dz)  # doctest: +SKIP
+    >>> nonlinear_power = halofit_smith(k, z, linear_power, cosmology)  # doctest: +SKIP
     '''
 
     # Manage shapes of input arrays
diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
index 6efa272..2b72cfd 100644
--- a/skypy/utils/photometry.py
+++ b/skypy/utils/photometry.py
@@ -5,6 +5,7 @@
 
 from abc import ABCMeta, abstractmethod
 import numpy as np
+import scipy.special
 
 
 __all__ = [
@@ -13,10 +14,12 @@ __all__ = [
     'luminosity_in_band',
     'mag_ab',
     'SpectrumTemplates',
+    'magnitude_error_rykoff',
+    'logistic_completeness_function',
 ]
 
 try:
-    import speclite.filters
+    import speclite.filters  # noqa F401
 except ImportError:
     HAS_SPECLITE = False
 else:
@@ -84,11 +87,12 @@ def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
     .. [1] M. R. Blanton et al., 2003, AJ, 125, 2348
 
     '''
+    from speclite.filters import load_filters
 
     # load the filters
     if np.ndim(filters) == 0:
         filters = (filters,)
-    filters = speclite.filters.load_filters(*filters)
+    filters = load_filters(*filters)
     if np.shape(filters) == (1,):
         filters = filters[0]
 
@@ -127,7 +131,7 @@ def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
         u = u.reshape(u.shape + (1,)*(nd_s+nd_f))
         m = np.ascontiguousarray(m[n])
         m += u*dm[n]
-        del(dm, n, u)
+        del dm, n, u
 
     # combine spectra if asked to
     if coefficients is not None:
@@ -303,3 +307,141 @@ def absolute_magnitude_from_luminosity(luminosity, zeropoint=None):
         zeropoint = -luminosity_in_band[zeropoint]
 
     return -2.5*np.log10(luminosity) - zeropoint
+
+
+def magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit=np.inf):
+    r"""Magnitude error according to the model from Rykoff et al. (2015).
+
+    Given an apparent magnitude calculate the magnitude error that is introduced
+    by the survey specifications and follows the model described in Rykoff et al. (2015).
+
+    Parameters
+    ----------
+    magnitude: array_like
+        Apparent magnitude. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    magnitude_limit: array_like
+        :math:`10\sigma` limiting magnitude of the survey. This and the other
+        array_like parameters must be broadcastable to the same shape.
+    magnitude_zp: array_like
+        Zero-point magnitude of the survey. This and the other array_like parameters must
+        be broadcastable to the same shape.
+    a,b: array_like
+        Model parameters: a is the intercept and
+        b is the slope of the logarithmic effective time.
+        These and the other array_like parameters must be broadcastable to the same shape.
+    error_limit: float, optional
+        Upper limit of the returned error. If given, all values larger than this value
+        will be set to error_limit. Default is None.
+
+    Returns
+    -------
+    error: ndarray
+        The apparent magnitude error in the Rykoff et al. (2015) model. This is a scalar
+        if magnitude, magnitude_limit, magnitude_zp, a and b are scalars.
+
+    Notes
+    -----
+    Rykoff et al. (2015) (see [1]_) describe the error of the apparent magnitude :math:`m` as
+
+    .. math::
+
+        \sigma_m(m;m_{\mathrm{lim}}, t_{\mathrm{eff}}) &=
+            \sigma_m(F(m);F_{\mathrm{noise}}(m_{\mathrm{lim}}), t_{\mathrm{eff}}) \\
+        &= \frac{2.5}{\ln(10)} \left[ \frac{1}{Ft_{\mathrm{eff}}}
+            \left( 1 + \frac{F_{\mathrm{noise}}}{F} \right) \right]^{1/2} \;,
+
+    where
+
+    .. math::
+
+        F=10^{-0.4(m - m_{\mathrm{ZP}})}
+
+    is the source's flux,
+
+    .. math::
+
+        F_\mathrm{noise} = \frac{F_{\mathrm{lim}}^2 t_{\mathrm{eff}}}{10^2} - F_{\mathrm{lim}}
+
+    is the effective noise flux and :math:`t_\mathrm{eff}` is the effective exposure time
+    (we absorbed the normalisation constant :math:`k` in the definition of
+    :math:`t_\mathrm{eff}`).
+    Furthermore, :math:`m_\mathrm{ZP}` is the zero-point magnitude of the survey and
+    :math:`F_\mathrm{lim}` is the :math:`10\sigma` limiting flux.
+    Accordingly, :math:`m_\mathrm{lim}` is the :math:`10\sigma` limiting magnitud
+    associated with :math:`F_\mathrm{lim}`.
+
+    The effective exposure time is described by
+
+    .. math::
+
+        \ln{t_\mathrm{eff}} = a + b(m_\mathrm{lim} - 21)\;,
+
+    where :math:`a` and :math:`b` are free parameters.
+
+    Further note that the model was originally used for SDSS galaxy photometry.
+
+    References
+    ----------
+    .. [1] Rykoff E. S., Rozo E., Keisler R., 2015, eprint arXiv:1509.00870
+
+    """
+
+    flux = luminosity_from_absolute_magnitude(magnitude, -magnitude_zp)
+    flux_limit = luminosity_from_absolute_magnitude(magnitude_limit, -magnitude_zp)
+    t_eff = np.exp(a + b * np.subtract(magnitude_limit, 21.0))
+    flux_noise = np.square(flux_limit / 10) * t_eff - flux_limit
+    error = 2.5 / np.log(10) * np.sqrt((1 + flux_noise / flux) / (flux * t_eff))
+
+    return np.minimum(error, error_limit)
+
+
+def logistic_completeness_function(magnitude, magnitude_95, magnitude_50):
+    r'''Logistic completeness function.
+
+    This function calculates the logistic completeness function (based on eq. (7) in
+    López-Sanjuan, C. et al. (2017) [1]_.
+
+    .. math::
+
+        p(m) = \frac{1}{1 + \exp[\kappa (m - m_{50})]}\;,
+
+    which describes the probability :math:`p(m)` that an object of magnitude :math:`m` is detected
+    in a specific band and with
+
+    .. math::
+
+        \kappa = \frac{\ln(\frac{1}{19})}{m_{95} - m_{50}}\;.
+
+    Here, :math:`m_{95}` and :math:`m_{50}` are the 95% and 50% completeness
+    magnitudes, respectively.
+
+    Parameters
+    ----------
+    magnitude : array_like
+        Magnitudes. Can be multidimensional for computing with multiple filter bands.
+    magnitude_95 : scalar or 1-D array_like
+        95% completeness magnitude.
+        If `magnitude_50` is 1-D array it has to be scalar or 1-D array of the same shape.
+    magnitude_50 : scalar or 1-D array_like
+        50% completeness magnitude.
+        If `magnitude_95` is 1-D array it has to be scalar or 1-D array of the same shape.
+
+    Returns
+    -------
+    probability : scalar or array_like
+        Probability of detecting an object with magnitude :math:`m`.
+        Returns array_like of the same shape as magnitude.
+        Exemption: If magnitude is scalar and `magnitude_95` or `magnitude_50`
+        is array_like of shape (nb, ) it returns array_like of shape (nb, ).
+
+    References
+    -----------
+    .. [1] López-Sanjuan, C. et al., `2017A&A…599A..62L`_
+    .. _2017A&A…599A..62L: https://ui.adsabs.harvard.edu/abs/2017A%26A...599A..62L
+
+    '''
+
+    kappa = np.log(1. / 19) / np.subtract(magnitude_95, magnitude_50)
+    arg = kappa * np.subtract(magnitude, magnitude_50)
+    return scipy.special.expit(-arg)
diff --git a/skypy/utils/random.py b/skypy/utils/random.py
index 30500a0..c301eee 100644
--- a/skypy/utils/random.py
+++ b/skypy/utils/random.py
@@ -41,13 +41,12 @@ def schechter(alpha, x_min, x_max, resolution=100, size=None, scale=1.):
     x_sample : array_like
         Samples drawn from the Schechter function.
 
-    Examples
+    Warnings
     --------
-    >>> import skypy.utils.random as random
-    >>> alpha = -1.3
-    >>> sample = random.schechter(alpha, x_min=1e-10, x_max=1e2,
-    ...                           resolution=100, size=1000)
-
+    The inverse cumulative distribution function is approximated from the
+    Schechter function evaluated on a logarithmically-spaced grid. The user
+    must choose the `resolution` of this grid to satisfy their desired
+    numerical accuracy.
 
     References
     ----------
diff --git a/skypy/utils/special.py b/skypy/utils/special.py
index dd51b4a..a2fae69 100644
--- a/skypy/utils/special.py
+++ b/skypy/utils/special.py
@@ -70,15 +70,6 @@ def gammaincc(a, x):
     --------
     scipy.special.gammaincc : Computes the start of the recurrence.
 
-    Examples
-    --------
-    This implementation of `gammaincc` supports positive and negative values
-    of `a`.
-
-    >>> from skypy.utils.special import gammaincc
-    >>> gammaincc([-1.5, -0.5, 0.5, 1.5], 1.2)
-    array([ 0.03084582, -0.03378949,  0.12133525,  0.49363462])
-
     '''
     if np.broadcast(a, x).ndim == 0:
         return _gammaincc(a, x)
diff --git a/tox.ini b/tox.ini
index 905f321..f0a9531 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,9 +1,9 @@
 [tox]
 envlist =
-    py{36,37,38,39}-test{,-all}{,-dev,-latest,-oldest}{,-cov}
-    py{36,37,38,39}-test-numpy{116,117,118,119}
-    py{36,37,38,39}-test-scipy{12,13,14,15,16}
-    py{36,37,38,39}-test-astropy{40,41,42}
+    py{37,38,39,310}-test{,-all}{,-dev,-latest,-oldest}{,-cov}
+    py{37,38,39,310}-test-numpy{116,117,118,119,120,121,122,123}
+    py{37,38,39,310}-test-scipy{12,13,14,15,16,17,18}
+    py{37,38,39,310}-test-astropy{40,41,42,43,50,51}
     build_docs
     linkcheck
     codestyle
@@ -17,7 +17,7 @@ indexserver =
 [testenv]
 
 # Pass through the following environment variables which may be needed for the CI
-passenv = HOME WINDIR LC_ALL LC_CTYPE CC CI TRAVIS
+passenv = HOME,WINDIR,LC_ALL,LC_CTYPE,CC,CI,TRAVIS
 
 # Run the tests in a temporary directory to make sure that we don't import
 # this package from the source tree
@@ -42,14 +42,23 @@ description =
     numpy117: with numpy 1.17.*
     numpy118: with numpy 1.18.*
     numpy119: with numpy 1.19.*
+    numpy120: with numpy 1.20.*
+    numpy121: with numpy 1.21.*
+    numpy122: with numpy 1.22.*
+    numpy123: with numpy 1.23.*
     scipy12: with scipy 1.2.*
     scipy13: with scipy 1.3.*
     scipy14: with scipy 1.4.*
     scipy15: with scipy 1.5.*
     scipy16: with scipy 1.6.*
+    scipy17: with scipy 1.7.*
+    scipy18: with scipy 1.8.*
     astropy40: with astropy 4.0.*
     astropy41: with astropy 4.1.*
     astropy42: with astropy 4.2.*
+    astropy43: with astropy 4.3.*
+    astropy50: with astropy 5.0.*
+    astropy51: with astropy 5.1.*
 
 # The following provides some specific pinnings for key packages
 deps =
@@ -58,24 +67,33 @@ deps =
     numpy117: numpy==1.17.*
     numpy118: numpy==1.18.*
     numpy119: numpy==1.19.*
+    numpy120: numpy==1.20.*
+    numpy121: numpy==1.21.*
+    numpy122: numpy==1.22.*
+    numpy123: numpy==1.23.*
 
     scipy12: scipy==1.2.*
     scipy13: scipy==1.3.*
     scipy14: scipy==1.4.*
     scipy15: scipy==1.5.*
     scipy16: scipy==1.6.*
+    scipy17: scipy==1.7.*
+    scipy18: scipy==1.8.*
 
     astropy40: astropy==4.0.*
     astropy41: astropy==4.1.*
     astropy42: astropy==4.2.*
+    astropy43: astropy==4.3.*
+    astropy50: astropy==5.0.*
+    astropy51: astropy==5.1.*
 
     dev: :NIGHTLY:numpy
     dev: :NIGHTLY:scipy
     dev: git+https://github.com/astropy/astropy.git#egg=astropy
 
-    latest: astropy==4.2.*
-    latest: numpy==1.19.*
-    latest: scipy==1.6.*
+    latest: astropy==5.1.*
+    latest: numpy==1.23.*
+    latest: scipy==1.8.*
 
     oldest: astropy==4.0.*
     oldest: numpy==1.16.*
@@ -86,6 +104,10 @@ extras =
     test
     all: all
 
+# Enable pip to install pre-releases in the dev environment
+pip_pre =
+    dev: True
+
 commands =
     pip freeze
     !cov: pytest --pyargs skypy {toxinidir}/docs {posargs}
@@ -97,12 +119,13 @@ description = invoke sphinx-build to build the HTML docs
 extras = docs
 commands =
     pip freeze
-    sphinx-build -W -b html . _build/html
+    sphinx-build -W -b html . _build/html {posargs}
 
 [testenv:linkcheck]
 changedir = docs
 description = check the links in the HTML docs
 extras = docs
+deps = pydot
 commands =
     pip freeze
     sphinx-build -W -b linkcheck . _build/html


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 1051df7..a1a452f 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -2,11 +2,11 @@ name: Tests
 on:
   push:
     branches:
-    - master
+    - main
     - module/*
   pull_request:
     branches:
-    - master
+    - main
     - module/*
 jobs:
   test:
@@ -19,24 +19,24 @@ jobs:
 
           - name: latest supported versions
             os: ubuntu-latest
-            python: 3.9
-            toxenv: py39-test-all-latest-cov
+            python: '3.10'
+            toxenv: py310-test-all-latest-cov
             toxposargs: --cov-report=xml:${GITHUB_WORKSPACE}/coverage.xml
 
           - name: oldest supported versions
             os: ubuntu-latest
-            python: 3.6
-            toxenv: py36-test-oldest
+            python: 3.7
+            toxenv: py37-test-oldest
 
           - name: macOS latest supported
             os: macos-latest
-            python: 3.9
-            toxenv: py39-test-latest
+            python: '3.10'
+            toxenv: py310-test-latest
 
           - name: Windows latest supported
             os: windows-latest
-            python: 3.9
-            toxenv: py39-test-latest
+            python: '3.10'
+            toxenv: py310-test-latest
 
     steps:
     - name: Checkout Repository
diff --git a/skypy/galaxies/tests/test_spectrum.py b/skypy/galaxies/tests/test_spectrum.py
index 5aea0f4..f529530 100644
--- a/skypy/galaxies/tests/test_spectrum.py
+++ b/skypy/galaxies/tests/test_spectrum.py
@@ -158,6 +158,39 @@ def test_kcorrect_stellar_mass():
     np.testing.assert_allclose(stellar_mass, truth)
 
 
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_kcorrect_stellar_mass_remain():
+
+    from astropy import units
+    from skypy.galaxies.spectrum import kcorrect
+    from speclite.filters import FilterResponse
+
+    # Gaussian bandpass
+    filt_lam = np.logspace(3, 4, 1000) * units.AA
+    filt_mean = 5000 * units.AA
+    filt_width = 100 * units.AA
+    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # Using the identity matrix for the coefficients yields trivial test cases
+    coeff = np.eye(5)
+    Mt = kcorrect.absolute_magnitudes(coeff, 'test-filt')
+
+    # Using the absolute magnitudes of the templates as reference magnitudes
+    # should return one solar mass for each template.
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mt, 'test-filt')
+    truth = kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+    # Solution for given magnitudes without template mixing
+    Mb = np.array([10, 20, 30, 40, 50])
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mb, 'test-filt')
+    truth = np.power(10, -0.4*(Mb-Mt)) * kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+
 def test_kcorrect_metallicity():
 
     from skypy.galaxies.spectrum import kcorrect
diff --git a/skypy/galaxies/tests/test_stellar_mass.py b/skypy/galaxies/tests/test_stellar_mass.py
index cb66798..b880923 100644
--- a/skypy/galaxies/tests/test_stellar_mass.py
+++ b/skypy/galaxies/tests/test_stellar_mass.py
@@ -3,6 +3,7 @@ import scipy.stats
 import scipy.integrate
 from scipy.special import gammaln
 import pytest
+from astropy.modeling.models import Exponential1D
 
 from skypy.galaxies import stellar_mass
 from skypy.utils import special
@@ -25,12 +26,30 @@ def test_stellar_masses():
     with pytest.raises(ValueError):
         stellar_mass.schechter_smf_mass(0., -1.4, np.array([1e10, 2e10]), 1.e7, 1.e13, size=3)
 
-    # Test that an array with the sme shape as m_star is returned if m_star is
+    # Test that an array with the same shape as m_star is returned if m_star is
     # an array and size = None
     m_star = np.array([1e10, 2e10])
     sample = stellar_mass.schechter_smf_mass(0., -1.4, m_star, 1.e7, 1.e13, size=None)
     assert m_star.shape == sample.shape
 
+    # Test m_star can be a function that returns an array of values for each redshift
+    redshift = np.linspace(0, 2, 100)
+    alpha = 0.357
+    m_min = 10 ** 7
+    m_max = 10 ** 14
+    m_star_function = Exponential1D(10**10.626, np.log(10)/0.095)
+    sample = stellar_mass.schechter_smf_mass(redshift, alpha, m_star_function, m_min, m_max)
+    assert sample.shape == redshift.shape
+
+    # Test alpha can be a function returning a scalar value.
+    sample = stellar_mass.schechter_smf_mass(redshift, lambda z: alpha, 10**10.626, m_min, m_max)
+    assert sample.shape == redshift.shape
+
+    # Sampling with an array for alpha is not implemented
+    alpha_array = np.full_like(redshift, alpha)
+    with pytest.raises(NotImplementedError, match='only scalar alpha is supported'):
+        sample = stellar_mass.schechter_smf_mass(redshift, alpha_array, m_star, m_min, m_max)
+
     # Test that sampling corresponds to sampling from the right pdf.
     # For this, we sample an array of luminosities for redshift z = 1.0 and we
     # compare it to the corresponding cdf.
diff --git a/skypy/pipeline/tests/data/bad_object.yml b/skypy/pipeline/tests/data/bad_object.yml
new file mode 100644
index 0000000..a9ff708
--- /dev/null
+++ b/skypy/pipeline/tests/data/bad_object.yml
@@ -0,0 +1,1 @@
+bad_object: !astropy.cosmology.Planck15 "bad value"
diff --git a/skypy/pipeline/tests/data/test_config.yml b/skypy/pipeline/tests/data/test_config.yml
index db7bef8..e9f4ce9 100644
--- a/skypy/pipeline/tests/data/test_config.yml
+++ b/skypy/pipeline/tests/data/test_config.yml
@@ -1,9 +1,10 @@
 test_int: 1
 test_float: 1.0
 test_str: hello world
-test_func: !numpy.random.uniform
-test_func_with_arg: !len 'hello world'
-test_cosmology: !astropy.cosmology.FlatLambdaCDM
+test_func: !numpy.random.uniform []
+test_func_with_arg: !len ['hello world']
+test_object: !astropy.cosmology.Planck15
+cosmology: !astropy.cosmology.FlatLambdaCDM
   H0: 67.74
   Om0: 0.3075
 tables:
diff --git a/skypy/pipeline/tests/test_config.py b/skypy/pipeline/tests/test_config.py
index d25236c..16b0b03 100644
--- a/skypy/pipeline/tests/test_config.py
+++ b/skypy/pipeline/tests/test_config.py
@@ -3,6 +3,7 @@ import pytest
 from skypy.pipeline import load_skypy_yaml
 from skypy.pipeline._items import Call
 from astropy import units
+from astropy.cosmology.core import Cosmology
 
 
 def test_load_skypy_yaml():
@@ -18,7 +19,9 @@ def test_load_skypy_yaml():
     assert isinstance(config['test_float'], float)
     assert isinstance(config['test_str'], str)
     assert isinstance(config['test_func'], Call)
-    assert isinstance(config['test_cosmology'], Call)
+    assert isinstance(config['test_func_with_arg'], Call)
+    assert isinstance(config['test_object'], Cosmology)
+    assert isinstance(config['cosmology'], Call)
     assert isinstance(config['tables']['test_table_1']['test_column_3'], Call)
 
     # Bad function
@@ -31,6 +34,11 @@ def test_load_skypy_yaml():
     with pytest.raises(ImportError):
         load_skypy_yaml(filename)
 
+    # Bad object
+    filename = get_pkg_data_filename('data/bad_object.yml')
+    with pytest.raises(ValueError):
+        load_skypy_yaml(filename)
+
 
 def test_empty_ref():
     filename = get_pkg_data_filename('data/test_empty_ref.yml')
diff --git a/skypy/pipeline/tests/test_pipeline.py b/skypy/pipeline/tests/test_pipeline.py
index 7c9674e..3b6151f 100644
--- a/skypy/pipeline/tests/test_pipeline.py
+++ b/skypy/pipeline/tests/test_pipeline.py
@@ -1,19 +1,26 @@
 from astropy.cosmology import FlatLambdaCDM, default_cosmology
 from astropy.cosmology.core import Cosmology
 from astropy.io import fits
-from astropy.table import Table
+from astropy.io.misc.hdf5 import read_table_hdf5
+from astropy.table import Table, vstack
 from astropy.table.column import Column
 from astropy.units import Quantity
 from astropy.utils.data import get_pkg_data_filename
 import networkx
 import numpy as np
-import os
 import pytest
 from skypy.pipeline import Pipeline
 from skypy.pipeline._items import Call, Ref
 
+try:
+    import h5py # noqa
+except ImportError:
+    HAS_H5PY = False
+else:
+    HAS_H5PY = True
 
-def test_pipeline():
+
+def test_pipeline(tmp_path):
 
     # Evaluate and store the default astropy cosmology.
     config = {'test_cosmology': Call(default_cosmology.get)}
@@ -36,17 +43,22 @@ def test_pipeline():
 
     pipeline = Pipeline(config)
     pipeline.execute()
-    pipeline.write(file_format='fits')
+    output_filename = str(tmp_path / 'output.fits')
+    pipeline.write(output_filename)
     assert len(pipeline['test_table']) == size
     assert np.all(pipeline['test_table.column1'] < pipeline['test_table.column2'])
-    with fits.open('test_table.fits') as hdu:
-        assert np.all(Table(hdu[1].data) == pipeline['test_table'])
+    with fits.open(output_filename) as hdu:
+        assert np.all(Table(hdu['test_table'].data) == pipeline['test_table'])
+
+    # Test invalid file extension
+    with pytest.raises(ValueError):
+        pipeline.write('output.invalid')
 
     # Check for failure if output files already exist and overwrite is False
     pipeline = Pipeline(config)
     pipeline.execute()
     with pytest.raises(OSError):
-        pipeline.write(file_format='fits', overwrite=False)
+        pipeline.write(output_filename, overwrite=False)
 
     # Check that the existing output files are modified if overwrite is True
     new_size = 2 * size
@@ -55,8 +67,8 @@ def test_pipeline():
     config['tables']['test_table']['column3'].args = [new_string]
     pipeline = Pipeline(config)
     pipeline.execute()
-    pipeline.write(file_format='fits', overwrite=True)
-    with fits.open('test_table.fits') as hdu:
+    pipeline.write(output_filename, overwrite=True)
+    with fits.open(output_filename) as hdu:
         assert len(hdu[1].data) == new_size
 
     # Check for failure if 'column1' requires itself creating a cyclic
@@ -202,7 +214,7 @@ def test_pipeline_read():
     assert isinstance(pipeline['test_int'], int)
     assert isinstance(pipeline['test_float'], float)
     assert isinstance(pipeline['test_str'], str)
-    assert isinstance(pipeline['test_cosmology'], Cosmology)
+    assert isinstance(pipeline['cosmology'], Cosmology)
     assert isinstance(pipeline['test_table_1'], Table)
     assert isinstance(pipeline['test_table_1']['test_column_3'], Column)
 
@@ -237,7 +249,43 @@ def test_column_quantity():
     np.testing.assert_array_less(pipeline['test_table.lengths_in_cm'], 100)
 
 
-def teardown_module(module):
+@pytest.mark.skipif(not HAS_H5PY, reason='Requires h5py')
+def test_hdf5(tmp_path):
+    size = 100
+    string = size*'a'
+    config = {'tables': {
+              'test_table': {
+                'column1': Call(np.random.uniform, [], {
+                  'size': size}),
+                'column2': Call(np.random.uniform, [], {
+                  'low': Ref('test_table.column1')}),
+                'column3': Call(list, [string], {})}}}
+
+    pipeline = Pipeline(config)
+    pipeline.execute()
+    output_filename = str(tmp_path / 'output.hdf5')
+    pipeline.write(output_filename)
+    hdf_table = read_table_hdf5(output_filename, 'tables/test_table', character_as_bytes=False)
+    assert np.all(hdf_table == pipeline['test_table'])
+
+
+def test_depends():
+
+    # Regression test for GitHub Issue #464
+    # Previously the .depends keyword was also being passed to functions as a
+    # keyword argument. This was because Pipeline was executing Item.infer to
+    # handle additional function arguments from context before handling
+    # additional dependencies specified using the .depends keyword. The
+    # .depends keyword is now handled first.
 
-    # Remove fits file generated in test_pipeline
-    os.remove('test_table.fits')
+    config = {'tables': {
+                'table_1': {
+                  'column1': Call(np.random.uniform, [0, 1, 10])},
+                'table_2': {
+                    '.init': Call(vstack, [], {
+                      'tables': [Ref('table_1')],
+                      '.depends': ['table_1.complete']})}}}
+
+    pipeline = Pipeline(config)
+    pipeline.execute()
+    assert np.all(pipeline['table_1'] == pipeline['table_2'])
diff --git a/skypy/pipeline/tests/test_skypy.py b/skypy/pipeline/tests/test_skypy.py
index 31d2969..14a66db 100644
--- a/skypy/pipeline/tests/test_skypy.py
+++ b/skypy/pipeline/tests/test_skypy.py
@@ -1,12 +1,15 @@
+from astropy.table import Table
 from astropy.utils.data import get_pkg_data_filename
 from contextlib import redirect_stdout
 from io import StringIO
 import pytest
 from skypy import __version__ as skypy_version
+from skypy.pipeline import load_skypy_yaml
+from skypy.pipeline._items import Call
 from skypy.pipeline.scripts import skypy
 
 
-def test_skypy():
+def test_skypy(tmp_path):
 
     # No arguments
     with pytest.raises(SystemExit) as e:
@@ -26,20 +29,81 @@ def test_skypy():
     assert version.getvalue().strip() == skypy_version
     assert e.value.code == 0
 
-    # Missing positional argument 'config'
+    # Missing positional argument 'output'
     with pytest.raises(SystemExit) as e:
-        skypy.main(['--format', 'fits'])
+        skypy.main(['config.filename'])
     assert e.value.code == 2
 
+    # Process empty config file
+    config_filename = get_pkg_data_filename('data/empty_config.yml')
+    output_filename = str(tmp_path / 'empty.fits')
+    assert skypy.main([config_filename, output_filename]) == 0
+
+    # Process test config file
+    config_filename = get_pkg_data_filename('data/test_config.yml')
+    output_filename = str(tmp_path / 'test.fits')
+    assert skypy.main([config_filename, output_filename]) == 0
+
     # Invalid file format
+    output_filename = str(tmp_path / 'test.invalid')
     with pytest.raises(SystemExit) as e:
-        skypy.main(['--format', 'invalid', 'config.filename'])
+        skypy.main([config_filename, output_filename])
     assert e.value.code == 2
 
-    # Process empty config file
-    filename = get_pkg_data_filename('data/empty_config.yml')
-    assert skypy.main([filename]) == 0
 
-    # Process test config file
-    filename = get_pkg_data_filename('data/test_config.yml')
-    assert skypy.main([filename]) == 0
+def test_logging(capsys, tmp_path):
+
+    # Run skypy with default verbosity and check log is empty
+    config_filename = get_pkg_data_filename('data/test_config.yml')
+    output_filename = str(tmp_path / 'logging.fits')
+    skypy.main([config_filename, output_filename])
+    out, err = capsys.readouterr()
+    assert not err
+
+    # Run again with increased verbosity and capture log. Force an exception by
+    # not using the "--overwrite" flag when the output file already exists.
+    with pytest.raises(SystemExit):
+        skypy.main([config_filename, output_filename, '--verbose'])
+    out, err = capsys.readouterr()
+
+    # Determine all DAG jobs and function calls from config
+    config = load_skypy_yaml(config_filename)
+    cosmology = config.pop('cosmology', None)
+    tables = config.pop('tables', {})
+    config.update({k: v.pop('.init', Call(Table)) for k, v in tables.items()})
+    columns = [f'{t}.{c}' for t, cols in tables.items() for c in cols]
+    functions = [f for f in config.values() if isinstance(f, Call)]
+    functions += [f for t, cols in tables.items() for f in cols.values() if isinstance(f, Call)]
+
+    # Check all jobs appear in the log
+    for job in list(config) + list(tables) + columns:
+        log_string = f"[INFO] skypy.pipeline: Generating {job}"
+        assert log_string in err
+
+    # Check all functions appear in the log
+    for f in functions:
+        log_string = f"[INFO] skypy.pipeline: Calling {f.function.__name__}"
+        assert log_string in err
+
+    # Check cosmology appears in the log
+    if cosmology:
+        assert "[INFO] skypy.pipeline: Setting cosmology" in err
+
+    # Check writing output file is in the log
+    assert f"[INFO] skypy: Writing {output_filename}" in err
+
+    # Check error for existing output file is in the log
+    try:
+        # New error message introduced in astropy PR #12179
+        from astropy.utils.misc import NOT_OVERWRITING_MSG
+        error_string = NOT_OVERWRITING_MSG.format(output_filename)
+    except ImportError:
+        # Fallback on old error message from astropy v4.x
+        error_string = f"[ERROR] skypy: File {output_filename!r} already exists."
+    assert error_string in err
+
+    # Run again with decreased verbosity and check the log is empty
+    with pytest.raises(SystemExit):
+        skypy.main([config_filename, output_filename, '-qq'])
+    out, err = capsys.readouterr()
+    assert not err
diff --git a/skypy/utils/tests/test_photometry.py b/skypy/utils/tests/test_photometry.py
index b617b5f..c6ec04d 100644
--- a/skypy/utils/tests/test_photometry.py
+++ b/skypy/utils/tests/test_photometry.py
@@ -161,3 +161,72 @@ def test_template_spectra():
                                                   Planck15, resolution=2)
     np.testing.assert_allclose(m_true, m_interp, rtol=1e-5)
     assert not np.all(m_true == m_interp)
+
+
+@pytest.mark.skipif(HAS_SPECLITE, reason='test requires that speclite is not installed')
+def test_speclite_not_installed():
+    """
+    Regression test for #436
+    Test that mag_ab raises the correct exception if speclite is not insalled.
+    """
+    from skypy.utils.photometry import mag_ab
+    wavelength = np.linspace(1, 10, 100)
+    spectrum = np.ones(10)
+    filter = 'bessell-B'
+    with pytest.raises(ImportError):
+        mag_ab(wavelength, spectrum, filter)
+
+
+def test_magnitude_error_rykoff():
+    from skypy.utils.photometry import magnitude_error_rykoff
+
+    # Test broadcasting to same shape given array for each parameter and
+    # test for correct result.
+    magnitude = np.full((2, 1, 1, 1, 1), 21)
+    magnitude_limit = np.full((3, 1, 1, 1), 21)
+    magnitude_zp = np.full((5, 1, 1), 21)
+    a = np.full((7, 1), np.log(200))
+    b = np.zeros(11)
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    # test result
+    assert np.allclose(error, 0.25 / np.log(10))
+    # test shape
+    assert error.shape == (2, 3, 5, 7, 11)
+
+    # second test for result
+    magnitude = 20
+    magnitude_limit = 22.5
+    magnitude_zp = 25
+    b = 2
+    a = np.log(10) - 1.5 * b
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b)
+    assert np.isclose(error, 0.25 / np.log(10) / np.sqrt(10))
+
+    # test that error limit is returned if error is larger than error_limit
+    # The following set-up would give a value larger than 10
+    magnitude = 30
+    magnitude_limit = 25
+    magnitude_zp = 30
+    a = 0.5
+    b = 1.0
+    error_limit = 1
+    error = magnitude_error_rykoff(magnitude, magnitude_limit, magnitude_zp, a, b, error_limit)
+    assert error == error_limit
+
+
+def test_logistic_completeness_function():
+    from skypy.utils.photometry import logistic_completeness_function
+
+    # Test that arguments broadcast correctly
+    m = np.full((2, 1, 1), 21)
+    m95 = np.full((3, 1), 22)
+    m50 = np.full(5, 23)
+    p = logistic_completeness_function(m, m95, m50)
+    assert p.shape == np.broadcast(m, m95, m50).shape
+
+    # Test result of completeness function for different given magnitudes
+    m95 = 24
+    m50 = 25
+    m = [np.finfo(np.float64).min, m95, m50, 2*m50-m95, np.finfo(np.float64).max]
+    p = logistic_completeness_function(m, m95, m50)
+    assert np.allclose(p, [1, 0.95, 0.5, 0.05, 0])
diff --git a/skypy/utils/tests/test_special.py b/skypy/utils/tests/test_special.py
index 565e70d..92b7f10 100644
--- a/skypy/utils/tests/test_special.py
+++ b/skypy/utils/tests/test_special.py
@@ -1,7 +1,7 @@
 import numpy as np
 import numpy.testing as npt
 from astropy.utils.data import get_pkg_data_filename
-from astropy.tests.helper import raises
+import pytest
 
 from skypy.utils.special import gammaincc
 
@@ -80,13 +80,13 @@ def test_gammaincc_precision():
     npt.assert_allclose(g, v, rtol=1e-10, atol=0, err_msg='\n'.join(lines))
 
 
-@raises(ValueError)
 def test_gammaincc_neg_x_scalar():
     # negative x raises an exception
-    gammaincc(0.5, -1.0)
+    with pytest.raises(ValueError):
+        gammaincc(0.5, -1.0)
 
 
-@raises(ValueError)
 def test_gammaincc_neg_x_array():
     # negative x in array raises an exception
-    gammaincc(0.5, [3.0, 2.0, 1.0, 0.0, -1.0])
+    with pytest.raises(ValueError):
+        gammaincc(0.5, [3.0, 2.0, 1.0, 0.0, -1.0])