Accessing evolutionary Coeval data

When you run run_coeval with non-saturated spin temperature fluctuations (or various other options, such as inhomogeneous recombinations), or you run run_lightcone, many coeval simulation cubes at higher redshifts are computed, to generate the evolution up to the redshift you requested (or in the case of lightcones, to be interpolated over). By default, with write=True, all these boxes are saved to the cache (or a folder of your choice). What is the best way to access that data?

Looking at the output datafiles themselves (a bunch of .h5 files) is fairly confusing – the files are saved with names based on inscrutable hashes. You could just run run_coeval at the redshift you care about, with the same input parameters – this will just return the cached object. But this is a bit clunky, and requires you to know the redshifts that were calculated during the evolution.

Never fear though – in v4 we provide a RunCache object that assists with understanding what’s in your cache.

import matplotlib.pyplot as plt

import py21cmfast as p21c
from tempfile import mkdtemp
import numpy as np

# For plotting the cubes, we use the plotting submodule:
# For interacting with the cache
from py21cmfast import plotting

We’re going to make a temporary cache, so that this notebook produces the same output every time.

If you want to hold on to each computed OutputStruct, you can specify your own folder in OutputCache and read/copy/write/delete the files as needed:

cache = p21c.OutputCache(mkdtemp())

To control which outputs are cached, one can pass a CacheConfig object to the write arguments of the lightcone and coeval functions.

cacheconfig = p21c.CacheConfig.on()
cacheconfig

CacheConfig(initial_conditions=True, perturbed_field=True, spin_temp=True, ionized_box=True, brightness_temp=True, halobox=True, perturbed_halo_catalog=True, halo_catalog=True, xray_source_box=True)

Passing True or False to the write argument turns caching for all objects on and off respectively

Create an example box

We first run a coeval box to have something to read later:

inputs = p21c.InputParameters.from_template(
    ['simple', 'small'],
    random_seed=1234,
    node_redshifts=np.arange(6, 21),
)

coevals = p21c.run_coeval(
    out_redshifts=(6,16),
    inputs=inputs,
    write=cacheconfig,
    cache=cache,
    progressbar=True
)

Accessing cached data from Coeval

If you have the exact object you want to load, you can simply pass it in to the find_existing method of the OutputCache to get the path, then read it using our hdf5 submodule

# Create an *empty* brightness temperature object
brightness_temp_z16 = p21c.BrightnessTemp.new(redshift=inputs.node_redshifts[10],inputs=inputs)

# Get the path to such an object within the cache.
path = cache.find_existing(brightness_temp_z16)

# Read the output struct (this is lazy -- it doesn't actually read any simulation data until you ask for it)
brightness_temp_z16 = p21c.read_output_struct(path)

The output is an OutputStruct of whatever kind you requested:

type(brightness_temp_z16)

py21cmfast.wrapper.outputs.BrightnessTemp

You can also Construct a RunCache instance, which contains the filenames of every box needed for the run, and find it there

runcache = p21c.RunCache.from_inputs(cache=cache, inputs=inputs)
brightness_temp_z16 = runcache.get_output_struct_at_z(
        kind="BrightnessTemp",
        z=16.0,
        match_z_within=0.01,
)

To list all objects currently stored in your cache, use list_datasets() with filters for output type, input parameters, redshift and seed. Pass None for each filter to output everything

cache.list_datasets(
    kind="PerturbedField",
    inputs=None,
    redshift=None,
    all_seeds=True,
)

[PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/12.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/19.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/6.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/16.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/10.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/9.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/15.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/17.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/7.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/13.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/8.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/18.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/20.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/14.0000/PerturbedField.h5'),
 PosixPath('/tmp/tmp39gegpyh/76bf7c537b8551236989120698767cf8/1234/01a01b1fe697e22e1682177d07c45412/11.0000/PerturbedField.h5')]

Note that the directory structure in the cache is split up into several hashes which represent various parts of the input structure and redshift. This ensures we don’t re-calculate the same boxes twice unless we explicitly pass regenerate=True to our functions

We can proceed to plot the brightness temperature that we read from cache

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
plotting.coeval_sliceplot(coevals[0], ax=ax[0], fig=fig)
plotting.coeval_sliceplot(brightness_temp_z16, ax=ax[1], fig=fig)
plt.tight_layout()

Of course, we could have pulled the cached data from a different redshift:

xH_z10 = runcache.get_output_struct_at_z(z=10., kind="IonizedBox")
xH_z6 = runcache.get_output_struct_at_z(z=7., kind="IonizedBox")

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
plotting.coeval_sliceplot(xH_z10, ax=ax[0], kind="neutral_fraction")
ax[0].set_title("z=10")
plotting.coeval_sliceplot(xH_z6, ax=ax[1], kind="neutral_fraction")
ax[1].set_title("z=6")
plt.tight_layout();

In fact, we could go further than this, and plot all the evaluated brightness temperatures. We can get the list of redshifts evaluated from the inputs directly as node_redshifts:

fig, ax = plt.subplots(
    len(inputs.node_redshifts) // 3,
    3,
    figsize=(14, 4 * len(inputs.node_redshifts) // 3),
)

for i, z in enumerate(inputs.node_redshifts):
    Tbz = runcache.get_output_struct_at_z(z=z, kind="BrightnessTemp")
    plotting.coeval_sliceplot(Tbz, ax=ax.flatten()[i])
    ax.flatten()[i].set_title(f"z={z:.2f}")
plt.tight_layout()

We can also directly construct Coeval objects using the cache:

coeval_z16 = runcache.get_coeval_at_z(
    z=8.,
)

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
plotting.coeval_sliceplot(coeval_z16, ax=ax[0],kind='brightness_temp')
plotting.coeval_sliceplot(coeval_z16, ax=ax[1],kind='neutral_fraction')
plt.tight_layout()

Coevals read in this way do not necessarily have to be in the out_redshifts of the coeval run, they can also be in the node_redshifts, or even computed in some other run with the same inputs and random seed.