Repository: skypyproject/skypy

Issue:
SHAM redshift slices
## Description
Implement SHAM on different redshift slices (snapshot) of a simulation.
## Inputs
 - Halo catalogues at different redshifts containing the halo proxies (halo mass, maximum circular velocity, xyzzy positions).
 - Galaxy catalogue containing the galaxy proxies (stellar mass, luminosity, RA/DEC positions).

## Outputs
 - Galaxy mock catalogues including both halo and galactic properties.

## References
https://arxiv.org/pdf/1509.06404

Discussion:
skypybot: Hello @rrjbca! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`skypy/pipeline/scripts/skypy.py`](https://github.com/skypyproject/skypy/blob/2b60f826290b04b753be235b3d85c5f35d12a588/skypy/pipeline/scripts/skypy.py):

> [Line 87:80](https://github.com/skypyproject/skypy/blob/2b60f826290b04b753be235b3d85c5f35d12a588/skypy/pipeline/scripts/skypy.py#L87): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (81 > 79 characters)



Patch:
diff --git a/examples/abundance_matching.yml b/examples/abundance_matching.yml
index dc7addf..2b17330 100644
--- a/examples/abundance_matching.yml
+++ b/examples/abundance_matching.yml
@@ -1,10 +1,13 @@
+lightcone:
+  z_min: 1
+  z_max: 2
+  n_slice: 10
 cosmology: !astropy.cosmology.default_cosmology.get
 A_s: 2.1982e-09
 n_s: 0.969453
-redshift: 1
 M_star: !astropy.modeling.models.Linear1D [-0.70798041, -20.37196157]
 growth_function: !skypy.power_spectrum.growth_function
-  redshift: $redshift
+  redshift: $slice_z_mid
   cosmology: $cosmology
 wavenumber: !numpy.logspace [-4, 2, 100]
 power_spectrum: !skypy.power_spectrum.eisenstein_hu
@@ -32,7 +35,7 @@ tables:
         m_min: 1.0E+9
         noise: True
       galaxy_kwargs:
-        redshift: $redshift
+        redshift: $slice_z_mid
         M_star: $M_star
         alpha: -0.5
         m_lim: 35
diff --git a/skypy/pipeline/__init__.py b/skypy/pipeline/__init__.py
index 48f02f8..3651ef9 100644
--- a/skypy/pipeline/__init__.py
+++ b/skypy/pipeline/__init__.py
@@ -5,3 +5,4 @@ dependencies and handle their outputs.
 
 from ._config import *  # noqa
 from ._pipeline import *  # noqa
+from ._lightcone import *  # noqa
diff --git a/skypy/pipeline/_lightcone.py b/skypy/pipeline/_lightcone.py
new file mode 100644
index 0000000..468e314
--- /dev/null
+++ b/skypy/pipeline/_lightcone.py
@@ -0,0 +1,112 @@
+from astropy.cosmology import default_cosmology, z_at_value
+from astropy.table import Table, vstack
+from copy import deepcopy
+import numpy as np
+from skypy.pipeline import Pipeline
+
+__all__ = [
+    'Lightcone',
+]
+
+
+class Lightcone:
+    r'''Class for running a lightcone simulation.
+
+    This class runs a simulaion pipeline for a set of independent redshift
+    slices and returns their combined results.
+    '''
+
+    def __init__(self, configuration):
+        '''Construct the pipeline.
+
+        Parameters
+        ----------
+        configuration : dict-like
+            Configuration for the lightcone simulation.
+
+        Notes
+        -----
+        'configuration' should contain an entry 'lightcone' which is a
+        dictionary defining 'z_min', 'z_max' and 'n_slice'. These are the
+        minimum and maximum redshift of the simulation and the number of
+        redshift slices respectively. The redshift range is subdivided into
+        slices of equal comoving distance and each slice is simulated as a
+        SkyPy Pipeline taking the other entries in 'configuraion'.
+        '''
+
+        self.config = deepcopy(configuration)
+        self.lightcone_config = self.config.pop('lightcone')
+
+    def execute(self, parameters={}):
+
+        # Lightcone parameters
+        z_min = self.lightcone_config['z_min']
+        z_max = self.lightcone_config['z_max']
+        n_slice = self.lightcone_config['n_slice']
+
+        params = {'z_min': z_min,
+                  'z_max': z_min,
+                  'slice_z_min': None,
+                  'slice_z_max': None,
+                  'slice_z_mid': None, }
+
+        # Additional user-defined parameters
+        params.update(parameters)
+
+        # Update config with ligthcone parameters and user parameters
+        if 'parameters' in self.config:
+            self.config['parameters'].update(params)
+        else:
+            self.config['parameters'] = params
+
+        # SkyPy Pipeline object
+        pipeline = Pipeline(self.config)
+
+        # Initialise empty tables
+        self.tables = {k: Table() for k in pipeline.table_config.keys()}
+
+        # Cosmology from pipeline
+        if pipeline.cosmology:
+            self.cosmology = pipeline.get_value(pipeline.cosmology)
+        else:
+            self.cosmology = default_cosmology.get()
+
+        # Calculate equispaced comoving distance slices in redshift space
+        chi_min = self.cosmology.comoving_distance(z_min)
+        chi_max = self.cosmology.comoving_distance(z_max)
+        chi = np.linspace(chi_min, chi_max, n_slice + 1)
+        chi_mid = (chi[:-1] + chi[1:]) / 2
+        z = [z_at_value(self.cosmology.comoving_distance, c, z_min, z_max) for c in chi[1:-1]]
+        z_mid = [z_at_value(self.cosmology.comoving_distance, c, z_min, z_max) for c in chi_mid]
+        redshift_slices = zip([z_min] + z, z + [z_max], z_mid)
+
+        # Simulate redshift slices and append results to tables
+        for slice_z_min, slice_z_max, slice_z_mid in redshift_slices:
+            slice_params = {'slice_z_min': slice_z_min,
+                            'slice_z_max': slice_z_max,
+                            'slice_z_mid': slice_z_mid}
+            pipeline.execute(parameters=slice_params)
+            for k, v in self.tables.items():
+                self.tables[k] = vstack((v, pipeline[k]))
+
+    def write(self, file_format=None, 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.
+        overwrite : bool
+            Whether to overwrite any existing files without warning.
+
+        References
+        ----------
+        .. [1] https://docs.astropy.org/en/stable/io/unified.html
+
+        '''
+        if file_format:
+            for name, data in self.tables.items():
+                filename = '.'.join((name, file_format))
+                data.write(filename, overwrite=overwrite)
diff --git a/skypy/pipeline/_pipeline.py b/skypy/pipeline/_pipeline.py
index 81e917b..5ff37b9 100644
--- a/skypy/pipeline/_pipeline.py
+++ b/skypy/pipeline/_pipeline.py
@@ -82,6 +82,9 @@ class Pipeline:
         self.config.update({k: v.pop('.init', default_table)
                             for k, v in self.table_config.items()})
 
+        # Initalise state with parameters
+        self.state = copy(self.parameters)
+
         # Create a Directed Acyclic Graph of all jobs and dependencies
         self.dag = networkx.DiGraph()
 
diff --git a/skypy/pipeline/scripts/skypy.py b/skypy/pipeline/scripts/skypy.py
index e157bfb..177d1d8 100644
--- a/skypy/pipeline/scripts/skypy.py
+++ b/skypy/pipeline/scripts/skypy.py
@@ -36,7 +36,7 @@ values as keyword arguments.
 
 import argparse
 from skypy import __version__ as skypy_version
-from skypy.pipeline import Pipeline
+from skypy.pipeline import Lightcone, Pipeline, load_skypy_yaml
 import sys
 
 
@@ -55,8 +55,12 @@ def main(args=None):
         args = sys.argv[1:]
 
     args = parser.parse_args(args or ['--help'])
+    config = load_skypy_yaml(args.config)
 
-    pipeline = Pipeline.read(args.config)
+    if 'lightcone' in config:
+        pipeline = Lightcone(config)
+    else:
+        pipeline = Pipeline(config)
     pipeline.execute()
     pipeline.write(file_format=args.format, overwrite=args.overwrite)
     return(0)


Test patch:
diff --git a/skypy/pipeline/tests/data/lightcone_config.yml b/skypy/pipeline/tests/data/lightcone_config.yml
new file mode 100644
index 0000000..8edef76
--- /dev/null
+++ b/skypy/pipeline/tests/data/lightcone_config.yml
@@ -0,0 +1,21 @@
+lightcone:
+  z_min: 1
+  z_max: 2
+  n_slice: 5
+cosmology: !astropy.cosmology.FlatLambdaCDM
+  H0: 67.74
+  Om0: 0.3075
+tables:
+  test_table_1:
+    test_column_1: !numpy.random.uniform
+      low: $z_min
+      high: $z_max
+      size: 10
+    test_column_2: !numpy.random.uniform
+      low: $slice_z_min
+      high: $slice_z_max
+      size: 10
+    test_column_3: !numpy.random.uniform
+      low: $slice_z_mid
+      high: $slice_z_max
+      size: 10
diff --git a/skypy/pipeline/tests/test_lightcone.py b/skypy/pipeline/tests/test_lightcone.py
new file mode 100644
index 0000000..393cace
--- /dev/null
+++ b/skypy/pipeline/tests/test_lightcone.py
@@ -0,0 +1,75 @@
+from astropy.cosmology import FlatLambdaCDM
+from astropy.io import fits
+from astropy.table import Table
+import numpy as np
+import os
+import pytest
+from skypy.pipeline import Lightcone
+
+
+def test_lightcone():
+
+    # Test that lightcone parameters are parsed and handled correctly
+    z_min, z_max, n_slice = 1, 2, 5
+    nz = 100
+    config = {'lightcone': {'z_min': z_min, 'z_max': z_max, 'n_slice': n_slice},
+              'tables':
+              {'test_table':
+               {'z1': (np.random.uniform, ['$slice_z_min', '$slice_z_max', nz]),
+                'z2': (np.random.uniform, ['$slice_z_mid', '$slice_z_max', nz])
+                }
+               }
+              }
+
+    # Sample values depending on the boundaries and midpoints of each slice
+    lightcone = Lightcone(config)
+    lightcone.execute()
+    chi1 = lightcone.cosmology.comoving_distance(lightcone.tables['test_table']['z1'])
+    chi2 = lightcone.cosmology.comoving_distance(lightcone.tables['test_table']['z2'])
+
+    # Calculate the comoving distance to the slice boundaries and midpoints
+    chi_min = lightcone.cosmology.comoving_distance(z_min)
+    chi_max = lightcone.cosmology.comoving_distance(z_max)
+    chi_edge = np.linspace(chi_min, chi_max, n_slice + 1)
+    chi_mid = (chi_edge[:-1] + chi_edge[1:]) / 2
+
+    # Check sampled values are consistent with slice boundaries and midpoints
+    for i, (low, hi, mid) in enumerate(zip(chi_edge[:-1], chi_edge[1:], chi_mid)):
+        np.testing.assert_array_less(low, chi1[i*nz:(i+1)*nz])
+        np.testing.assert_array_less(chi1[i*nz:(i+1)*nz], hi)
+        np.testing.assert_array_less(mid, chi2[i*nz:(i+1)*nz])
+        np.testing.assert_array_less(chi2[i*nz:(i+1)*nz], hi)
+
+    # Repeat test with non-default cosmology
+    config['parameters'] = {'H0': 70, 'Om0': 0.3}
+    config['cosmology'] = (FlatLambdaCDM, {'H0': '$H0', 'Om0': '$Om0'})
+    lightcone = Lightcone(config)
+    lightcone.execute()
+    chi1 = lightcone.cosmology.comoving_distance(lightcone.tables['test_table']['z1'])
+    chi2 = lightcone.cosmology.comoving_distance(lightcone.tables['test_table']['z2'])
+
+    # Calculate the comoving distance to the slice boundaries and midpoints
+    chi_min = lightcone.cosmology.comoving_distance(z_min)
+    chi_max = lightcone.cosmology.comoving_distance(z_max)
+    chi_edge = np.linspace(chi_min, chi_max, n_slice + 1)
+    chi_mid = (chi_edge[:-1] + chi_edge[1:]) / 2
+
+    # Check sampled values are consistent with slice boundaries and midpoints
+    for i, (low, hi, mid) in enumerate(zip(chi_edge[:-1], chi_edge[1:], chi_mid)):
+        np.testing.assert_array_less(low, chi1[i*nz:(i+1)*nz])
+        np.testing.assert_array_less(chi1[i*nz:(i+1)*nz], hi)
+        np.testing.assert_array_less(mid, chi2[i*nz:(i+1)*nz])
+        np.testing.assert_array_less(chi2[i*nz:(i+1)*nz], hi)
+
+    # Test writing to file and overwrite failure
+    lightcone.write(file_format='fits')
+    with fits.open('test_table.fits') as hdu:
+        assert np.all(Table(hdu[1].data) == lightcone.tables['test_table'])
+    with pytest.raises(OSError):
+        lightcone.write(file_format='fits', overwrite=False)
+
+
+def teardown_module(module):
+
+    # Remove fits file generated in test_lightcone
+    os.remove('test_table.fits')
diff --git a/skypy/pipeline/tests/test_pipeline.py b/skypy/pipeline/tests/test_pipeline.py
index ef4b58c..a8a034a 100644
--- a/skypy/pipeline/tests/test_pipeline.py
+++ b/skypy/pipeline/tests/test_pipeline.py
@@ -1,6 +1,9 @@
 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.table.column import Column
+from astropy.utils.data import get_pkg_data_filename
 import networkx
 import numpy as np
 import os
@@ -179,6 +182,20 @@ def test_pipeline_cosmology():
     assert default_cosmology.get() == initial_default
 
 
+def test_pipeline_read():
+
+    # Test reading config from a file
+    filename = get_pkg_data_filename('data/test_config.yml')
+    pipeline = Pipeline.read(filename)
+    pipeline.execute()
+    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['test_table_1'], Table)
+    assert isinstance(pipeline['test_table_1']['test_column_3'], Column)
+
+
 def teardown_module(module):
 
     # Remove fits file generated in test_pipeline
diff --git a/skypy/pipeline/tests/test_skypy.py b/skypy/pipeline/tests/test_skypy.py
index 31d2969..8bc191c 100644
--- a/skypy/pipeline/tests/test_skypy.py
+++ b/skypy/pipeline/tests/test_skypy.py
@@ -43,3 +43,7 @@ def test_skypy():
     # Process test config file
     filename = get_pkg_data_filename('data/test_config.yml')
     assert skypy.main([filename]) == 0
+
+    # Process lightcone config file
+    filename = get_pkg_data_filename('data/lightcone_config.yml')
+    assert skypy.main([filename]) == 0