Complex Library Generation

In addition to the basic library generation demonstrated previously, Synference also supports more complex scenarios, such as libraries with more supplementary parameters, custom observation transformations, and more complicated model setups. Below is an example of how to generate a more complex model library using Synference.

Firstly we’ll just import the necessary modules and set up the synthesizer model.

[1]:
import numpy as np
from astropy.cosmology import Planck18 as cosmo
from synthesizer.emission_models import Greybody
from synthesizer.emission_models.attenuation import (
    Calzetti2000,
)
from synthesizer.emission_models.stellar.pacman_model import (
    PacmanEmission,
)
from synthesizer.grid import Grid
from synthesizer.instruments import FilterCollection, Instrument
from synthesizer.parametric import ZDist
from tqdm import tqdm
from unyt import K

And our Synference components - the GalaxyBasis class, and some utility functions we will use later.

[2]:
from synference import (
    GalaxyBasis,
    calculate_balmer_decrement,
    calculate_beta,
    calculate_colour,
    calculate_d4000,
    calculate_mass_weighted_age,
    calculate_muv,
    calculate_sfh_quantile,
    draw_from_hypercube,
    generate_random_DB_sfh,
)
/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

We’ll set up Synthesizer similarly to before - see the basic library generation example for more details.

For this example we’ll use a set of filter used in wide area surveys, including VISTA, Subaru Hyper Suprime-Cam, Euclid, and Spitzer IRAC.

Note that if you’re running this step on a cluster node without internet access, you’ll need to create the instrument file beforehand and pass in the HDF5 file path insead.

[3]:
filter_codes = [
    "Paranal/VISTA.Z",
    "Paranal/VISTA.Y",
    "Paranal/VISTA.J",
    "Paranal/VISTA.H",
    "Paranal/VISTA.Ks",
    "Subaru/HSC.g",
    "Subaru/HSC.r",
    "Subaru/HSC.i",
    "Subaru/HSC.z",
    "Subaru/HSC.Y",
    "CFHT/MegaCam.u",
    "CFHT/MegaCam.g",
    "CFHT/MegaCam.r",
    "CFHT/MegaCam.i",
    "CFHT/MegaCam.z",
    "Euclid/VIS.vis",
    "Euclid/NISP.Y",
    "Euclid/NISP.J",
    "Euclid/NISP.H",
    "Spitzer/IRAC.I1",
    "Spitzer/IRAC.I2",
]

filterset = FilterCollection(filter_codes=filter_codes)

instrument = Instrument("EuclidDeep", filters=filterset)

grid = Grid("test_grid")
---------------------------------------------------------------------------
TimeoutError                              Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:1348, in AbstractHTTPHandler.do_open(self, http_class, req, **http_conn_args)
   1347 try:
-> 1348     h.request(req.get_method(), req.selector, req.data, headers,
   1349               encode_chunked=req.has_header('Transfer-encoding'))
   1350 except OSError as err: # timeout error

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:1303, in HTTPConnection.request(self, method, url, body, headers, encode_chunked)
   1302 """Send a complete request to the server."""
-> 1303 self._send_request(method, url, body, headers, encode_chunked)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:1349, in HTTPConnection._send_request(self, method, url, body, headers, encode_chunked)
   1348     body = _encode(body, 'body')
-> 1349 self.endheaders(body, encode_chunked=encode_chunked)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:1298, in HTTPConnection.endheaders(self, message_body, encode_chunked)
   1297     raise CannotSendHeader()
-> 1298 self._send_output(message_body, encode_chunked=encode_chunked)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:1058, in HTTPConnection._send_output(self, message_body, encode_chunked)
   1057 del self._buffer[:]
-> 1058 self.send(msg)
   1060 if message_body is not None:
   1061
   1062     # create a consistent interface to message_body

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:996, in HTTPConnection.send(self, data)
    995 if self.auto_open:
--> 996     self.connect()
    997 else:

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/http/client.py:962, in HTTPConnection.connect(self)
    961 sys.audit("http.client.connect", self, self.host, self.port)
--> 962 self.sock = self._create_connection(
    963     (self.host,self.port), self.timeout, self.source_address)
    964 # Might fail in OSs that don't implement TCP_NODELAY

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/socket.py:857, in create_connection(address, timeout, source_address)
    856 try:
--> 857     raise err
    858 finally:
    859     # Break explicitly a reference cycle

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/socket.py:845, in create_connection(address, timeout, source_address)
    844     sock.bind(source_address)
--> 845 sock.connect(sa)
    846 # Break explicitly a reference cycle

TimeoutError: [Errno 110] Connection timed out

During handling of the above exception, another exception occurred:

URLError                                  Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/instruments/filters.py:2041, in Filter._make_svo_filter(self)
   2040 try:
-> 2041     with urllib.request.urlopen(self.svo_url) as f:
   2042         # Get the root of the XML tree
   2043         root = ElementTree.parse(f).getroot()

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:216, in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    215     opener = _opener
--> 216 return opener.open(url, data, timeout)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:519, in OpenerDirector.open(self, fullurl, data, timeout)
    518 sys.audit('urllib.Request', req.full_url, req.data, req.headers, req.get_method())
--> 519 response = self._open(req, data)
    521 # post-process response

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:536, in OpenerDirector._open(self, req, data)
    535 protocol = req.type
--> 536 result = self._call_chain(self.handle_open, protocol, protocol +
    537                           '_open', req)
    538 if result:

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:496, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    495 func = getattr(handler, meth_name)
--> 496 result = func(*args)
    497 if result is not None:

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:1377, in HTTPHandler.http_open(self, req)
   1376 def http_open(self, req):
-> 1377     return self.do_open(http.client.HTTPConnection, req)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/urllib/request.py:1351, in AbstractHTTPHandler.do_open(self, http_class, req, **http_conn_args)
   1350 except OSError as err: # timeout error
-> 1351     raise URLError(err)
   1352 r = h.getresponse()

URLError: <urlopen error [Errno 110] Connection timed out>

During handling of the above exception, another exception occurred:

SVOInaccessible                           Traceback (most recent call last)
Cell In[3], line 25
      1 filter_codes = [
      2     "Paranal/VISTA.Z",
      3     "Paranal/VISTA.Y",
   (...)
     22     "Spitzer/IRAC.I2",
     23 ]
---> 25 filterset = FilterCollection(filter_codes=filter_codes)
     27 instrument = Instrument("EuclidDeep", filters=filterset)
     29 grid = Grid("test_grid")

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/utils/operation_timers.py:98, in timed.<locals>.decorator.<locals>.wrapped(*args, **kwargs)
     94 tic(timer_name)
     95 try:
     96     # Return the wrapped function result unchanged so the decorator
     97     # is transparent aside from its timing side effect.
---> 98     return func(*args, **kwargs)
     99 finally:
    100     # Always stop the timer, even if the wrapped function raises,
    101     # so the timing stack remains balanced.
    102     toc(timer_name)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/instruments/filters.py:297, in FilterCollection.__init__(self, filter_codes, tophat_dict, generic_dict, filters, path, new_lam, fill_gaps, verbose)
    295 # Let's make the filters
    296 if filter_codes is not None:
--> 297     self._include_svo_filters(filter_codes)
    298 if tophat_dict is not None:
    299     self._include_top_hat_filters(tophat_dict)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/instruments/filters.py:474, in FilterCollection._include_svo_filters(self, filter_codes)
    471 # Loop over the given filter codes
    472 for f in filter_codes:
    473     # Get filter from SVO
--> 474     _filter = Filter(f, new_lam=self.lam)
    476     # Store the filter and its code
    477     self.filters[_filter.filter_code] = _filter

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/units.py:909, in accepts.<locals>.check_accepts.<locals>.wrapped(*args, **kwargs)
    906 finally:
    907     toc(f"accepts({func.__qualname__})")
--> 909 return func(*bound.args, **bound.kwargs)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/utils/operation_timers.py:98, in timed.<locals>.decorator.<locals>.wrapped(*args, **kwargs)
     94 tic(timer_name)
     95 try:
     96     # Return the wrapped function result unchanged so the decorator
     97     # is transparent aside from its timing side effect.
---> 98     return func(*args, **kwargs)
     99 finally:
    100     # Always stop the timer, even if the wrapped function raises,
    101     # so the timing stack remains balanced.
    102     toc(timer_name)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/instruments/filters.py:1770, in Filter.__init__(self, filter_code, transmission, lam_min, lam_max, lam_eff, lam_fwhm, new_lam, hdf)
   1768 # Is this an SVO filter?
   1769 elif "/" in filter_code and "." in filter_code:
-> 1770     self._make_svo_filter()
   1772 # Otherwise we haven't got a valid combination of inputs.
   1773 else:
   1774     raise exceptions.InconsistentArguments(
   1775         "Invalid combination of filter inputs. \n For a generic "
   1776         "filter provide a transmission and wavelength array. "
   (...)
   1781         "wavelength or an effective wavelength and FWHM."
   1782     )

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/synthesizer/instruments/filters.py:2052, in Filter._make_svo_filter(self)
   2049         data = root.find(".//TABLEDATA")
   2051 except URLError:
-> 2052     raise exceptions.SVOInaccessible(
   2053         (
   2054             f"The SVO Database at {self.svo_url} "
   2055             "is not responding. Is it down?"
   2056         )
   2057     )
   2059 # Throw an error if we didn't find the filter.
   2060 if field is None:

SVOInaccessible: The SVO Database at http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=Paranal/VISTA.Z is not responding. Is it down?

Now we can configure our model, which we will make more complex. Firstly we define the prior ranges for our main galaxy parameters.

This model will use a more complex dust attenuation model, with variable UV slope and 2175A bump strength, as well as a variable escape fraction for ionizing photons.

[4]:
N_models = 1000

redshift = (0.01, 14)
masses = (4, 12)  # log(M/Msun)=4 to log(M/Msun)=12
logAv = (-3, 0.7)  # Av=0.001 to Av=5
log_zmet = (-4, -1.39)  # Z=0.0001 to Z=0.04
fesc = (0.0, 1.0)  # Fraction of ionizing photons that escape the galaxy
slope = (-0.4, 1.1)  # UV slope modification to Calzetti law
bump_strength = (0.0, 3.0)  # 2175A bump strength in Calzetti law


prior_ranges = {
    "redshift": redshift,
    "log_masses": masses,
    "log_Av": logAv,  # Av in magnitudes
    "log_zmet": log_zmet,
    "fesc": fesc,
    "slope": slope,
    "bump_strength": bump_strength,
}

Star Formation History

We’ll use a non-parametric ‘Dense Basis’ SFH (Iyer et al. 2019), where we model the time at which different quantiles of stellar mass formed. Synference provides a helper module for generating these Star Formation Histories.

We will use 3 quantiles, e.g. \(t_{25}, t_{50}, t_{75}\), so we will simply add three dummy parameters to our LHC sampling of the parameter space, and then replace them afterward. We need to set the concentration parameter \(\alpha\) for the Dirichilet prior on the SFH, for which we will use \(\alpha=3\). This controls the correlation between different SFH quantiles, where lower values will have more rapdily varying star formation histories. We also set a prior on the recent SFR, in terms of the sSFR, which normalizes by the stellar mass. This allows for a range from quiescent to highly star-forming galaxies.

[5]:
tx_alpha = 3

for i in range(3):
    j = 100 * (i + 1) / (4)
    prior_ranges[f"sfh_quantile_{j:.0f}"] = (0, 1)

prior_ranges["ssfr"] = (-14, -7)

Now we will sample these parameters from our hypercube. Note that we set a log prior on dust attenuation \(A_V\), but we want to sample linear dust attenuation, so we get the model to ‘unlog’ the parameter.

[6]:
all_param_dict = draw_from_hypercube(prior_ranges, N_models, unlog_keys=["log_Av"])

Now we will create our metallicity distributions.

[7]:
Z_dists = [
    ZDist.DeltaConstant(log10metallicity=log_z)
    for log_z in tqdm(all_param_dict["log_zmet"], desc="Creating ZDist")
]
Creating ZDist: 100%|██████████| 1000/1000 [00:00<00:00, 486860.59it/s]

Now we can create the star-formation history. We are using a specific prior in specific star-formation rate here, but we need to provide log SFR to the function, so we calculate this inside the loop.

[8]:
# Draw SFH params from prior
sfh_models = []
for i in tqdm(range(N_models), desc="Generating SFH models"):
    z = all_param_dict["redshift"][i]
    logmass = all_param_dict["log_masses"][i]
    logssfr = all_param_dict["ssfr"][i]
    logsfr = logmass + logssfr
    sfh, tx = generate_random_DB_sfh(
        Nparam=3,
        tx_alpha=tx_alpha,
        redshift=z,
        logsfr=logsfr,
        logmass=logmass,
    )
    for j in range(3):
        all_param_dict[f"sfh_quantile_{100 * (j + 1) / (3 + 1):.0f}"][i] = tx[j]
    sfh_models.append(sfh)
Generating SFH models:   4%|▍         | 39/1000 [00:00<00:10, 95.84it/s]
Starting dense_basis. Failed to load FSPS, only GP-SFH module will be available.
Generating SFH models: 100%|██████████| 1000/1000 [00:03<00:00, 327.95it/s]

Now we will set up two functions, which are used to convert parameters from one form to another inside the model, and allow the SFH model to be serialized and re-created. This is a flexible system which should allow complex transformations.

These are required because synference automatically looks for varying parameters which are stored on the Synthesizer galaxies and emitters, but if they are in a complex form (e.g. the Dense Basis tuple), then we need to explain how to understand them to the code.

[9]:
def make_db_tuple(params):
    """Constructs the DenseBasis tuple from the SFH."""
    nquant = 0
    for key in params:
        if key.startswith("sfh_quantile_"):
            nquant += 1

    mass_quantiles = np.linspace(0, 1, nquant + 2)[1:-1]  # Exclude the 0 and 1 quantiles

    db_tuple = [params["log_mass"], params["log_sfr"], nquant] + [
        params[f"sfh_quantile_{int(q * 100)}"] for q in mass_quantiles
    ]
    return db_tuple  # Return a tuple of (log_mass, SFR, nquant, [quantiles...])


def db_sf_convert(param, param_dict, Nparam_SFH=3):
    """Converts from a DenseBasis tuple back to parameters."""
    db_tuple = param_dict["db_tuple"]
    # dp_tuple has the folliwng
    # mass, sfr, tx_alpha, *sfh_quantiles
    if param.startswith("sfh_quantile_"):
        # Convert the SFH quantile parameters to the Dense Basis SFH format
        j = int(np.round(int(param.split("_")[-1]) / 100 * (Nparam_SFH + 1)))
        return db_tuple[j + 2]  # +3 because first three are mass, sfr, tx_alpha
    elif param == "log_sfr":
        # Convert log_sfr to the Dense Basis SFH format
        return db_tuple[1]
    elif param == "log_masses":
        # Convert log_masses to the Dense Basis SFH format
        return db_tuple[0]
    elif param == "tx_alpha":
        # Convert tx_alpha to the Dense Basis SFH format
        return db_tuple[2]
    elif param == "log_ssfr":
        # Convert log_ssfr to the Dense Basis SFH format
        return db_tuple[1] - db_tuple[0]
    else:
        raise ValueError(f"Unknown parameter {param.str} in db_tuple conversion.")

Now we can set up some parameter transformations. These are for when we want to sample a parameter which is not what is used directly in Synthesizer. The simplest example is our use of the V-band attenuation \(A_V\), whereas in reality to generate a model we must provide the V-band optical depth \(\tau_V\), to Synthesizer. So we provide a function to Synference to allow it to do this conversion, which is simply \(A_V = \tau_V * 1.086\)

These generally take the form of a dictionary, where the key is the parameter name required by the Synthesizer model, and value is a two component tuple, where the first value is a new name (or a list of new names, for multiple parameters), and the second value the conversion function.

We will also save the inverse functions in a separate dictionary. These are only used if we want to recreate this simulator later, to generate SEDs from a set of input parameters. This is useful to recover SEDs for observations when performing inference.

[10]:
alt_parametrizations = {
    "tau_v": ("Av", lambda x: x["tau_v"] * 1.086),
    "db_tuple": (
        ["log_sfr"] + [f"sfh_quantile_{100 * (j + 1) / (3 + 1):.0f}" for j in range(3)],
        db_sf_convert,
    ),
}

param_transforms_to_save = {
    "tau_v": lambda x: x["Av"] / 1.086,
    "db_tuple": make_db_tuple,
}

Emission Models

Now we will set up our complex emission model which supports our priors, with variable escape fraction and flexible attenuation law.

The basic concept is simple: Any emission model parameter set with a string, rather than an explicit value, will be inherited from the emission model, or emitter. So in this case, for ‘tau_v’, ‘fesc’, ‘bump_strength’, and ‘slope’, Synthesizer will look for these parameters to be set on the individual Galaxy or Star instances.

We will also set the emission key we will save, which is the root of the emission model, named ‘total’.

[11]:
dust_emission = Greybody(temperature=40 * K, emissivity=1.5)
dust_curve = Calzetti2000(slope="slope", ampl="bump_strength")

print("Creating emission model.")
emission_model = PacmanEmission(
    grid=grid, tau_v="tau_v", dust_curve=dust_curve, dust_emission=dust_emission, fesc="fesc"
)

emission_key = emission_model.label

print(f"Root emission model label is: {emission_key}")
Creating emission model.
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 6
      2 dust_curve = Calzetti2000(slope="slope", ampl="bump_strength")
      4 print("Creating emission model.")
      5 emission_model = PacmanEmission(
----> 6     grid=grid, tau_v="tau_v", dust_curve=dust_curve, dust_emission=dust_emission, fesc="fesc"
      7 )
      9 emission_key = emission_model.label
     11 print(f"Root emission model label is: {emission_key}")

NameError: name 'grid' is not defined

To ensure these parameters are set on each galaxy, we will create our final input, the galaxy_params dictionary. This is simply a dictionary of parameter name and value array pairs for every galaxy.

[12]:
galaxy_params = {
    "fesc": all_param_dict["fesc"],
    "tau_v": all_param_dict["Av"] / 1.086,
    "bump_strength": all_param_dict["bump_strength"],
    "slope": all_param_dict["slope"],
}

Now we can instantiate the GalaxyBasis, into which we will pass these inputs. This won’t do much until we call the correct function to build the library. Note that we set ‘build_library’ = False, because we have already generated our full parameters sample. If we wanted we could also pass in a smaller set of parameter values instead, and set build_library=True, and the code would generate all the combinations of those parameters.

[13]:
basis = GalaxyBasis(
    model_name="sps_Euclid_test",
    redshifts=all_param_dict["redshift"],
    grid=grid,
    emission_model=emission_model,
    sfhs=sfh_models,
    cosmo=cosmo,
    instrument=instrument,
    metal_dists=Z_dists,
    galaxy_params=galaxy_params,
    alt_parametrizations=alt_parametrizations,
    redshift_dependent_sfh=True,
    build_library=False,
    log_stellar_masses=all_param_dict["log_masses"],
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 4
      1 basis = GalaxyBasis(
      2     model_name="sps_Euclid_test",
      3     redshifts=all_param_dict["redshift"],
----> 4     grid=grid,
      5     emission_model=emission_model,
      6     sfhs=sfh_models,
      7     cosmo=cosmo,
      8     instrument=instrument,
      9     metal_dists=Z_dists,
     10     galaxy_params=galaxy_params,
     11     alt_parametrizations=alt_parametrizations,
     12     redshift_dependent_sfh=True,
     13     build_library=False,
     14     log_stellar_masses=all_param_dict["log_masses"],
     15 )

NameError: name 'grid' is not defined

Something else we can do here to improve the utility of our model is add more parameters to be saved and stored. Synference provides a set of these parameters, to save things like the surviving stellar mass, UV magnitude, \(\beta\) slope, D4000 break strength, UVJ colors, etc. We can see the full list here:

[14]:
from synference import SUPP_FUNCTIONS

SUPP_FUNCTIONS()
[14]:
Available supplementary functions:
  - calculate_MUV
  - calculate_Ndot_ion
  - calculate_agn_fraction
  - calculate_balmer_decrement
  - calculate_beta
  - calculate_burstiness
  - calculate_colour
  - calculate_d4000
  - calculate_flux_weighted_age
  - calculate_line_ew
  - calculate_line_flux
  - calculate_lum_weighted_age
  - calculate_mass_weighted_age
  - calculate_muv
  - calculate_sfh_quantile
  - calculate_sfr
  - calculate_surviving_mass
  - calculate_xi_ion0

We can pass in these functions to our GalaxyBasis.create_mock_library function, and they will be run for every galaxy and the output stored in the library.

The functions should take the galaxy as the first argument, and then the following arguments (if any) will be set by position as demonstrated below. The keys will be the parameter names. The return should be either a single value (float, string, unyt_quantity) or a dictionary which maps to the emission model names (e.g to record \(\beta\) slope for different components)

[15]:
supp_params = {
    "mUV": (calculate_muv, cosmo),
    "sfh_quant_25": (calculate_sfh_quantile, 0.25, True),  # Calculate SFH quantile at 25%
    "sfh_quant_50": (calculate_sfh_quantile, 0.50, True),  # Calculate SFH quantile at 50%
    "sfh_quant_75": (calculate_sfh_quantile, 0.75, True),  # Calculate SFH quantile at 75%
    "UV": (calculate_colour, "U", "V", emission_key, True),  # Calculate UV colour (rest-frame)
    "VJ": (calculate_colour, "V", "J", emission_key, True),  # Calculate VJ colour (rest-frame)
    "d4000": (calculate_d4000, emission_key),  # Calculate D4000 using the emission model
    "beta": (calculate_beta, emission_key),
    "balmer_decrement": (calculate_balmer_decrement, emission_key),
    "mass_weighted_age": calculate_mass_weighted_age,
}
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 6
      1 supp_params = {
      2     "mUV": (calculate_muv, cosmo),
      3     "sfh_quant_25": (calculate_sfh_quantile, 0.25, True),  # Calculate SFH quantile at 25%
      4     "sfh_quant_50": (calculate_sfh_quantile, 0.50, True),  # Calculate SFH quantile at 50%
      5     "sfh_quant_75": (calculate_sfh_quantile, 0.75, True),  # Calculate SFH quantile at 75%
----> 6     "UV": (calculate_colour, "U", "V", emission_key, True),  # Calculate UV colour (rest-frame)
      7     "VJ": (calculate_colour, "V", "J", emission_key, True),  # Calculate VJ colour (rest-frame)
      8     "d4000": (calculate_d4000, emission_key),  # Calculate D4000 using the emission model
      9     "beta": (calculate_beta, emission_key),
     10     "balmer_decrement": (calculate_balmer_decrement, emission_key),
     11     "mass_weighted_age": calculate_mass_weighted_age,
     12 }

NameError: name 'emission_key' is not defined

Library Generation

Now we can run the final method to create the output catalogue. There are several things to note here. We can set the number of processes to use, to make use of multiple threads. We can also set the batch size, which will split the generation into multiple HDF5 files which are later combined. This is useful to avoid running out of RAM with large libraries.

We can also set the output type - in this case it is “photometry”, but if instead we made in “spectra”, we would generate a library of spectra which we could infer from as well.

[16]:
basis.create_mock_library(
    emission_model_key=emission_key,
    out_name="library_Euclid_test",
    out_dir="./",
    overwrite=True,
    n_proc=1,
    verbose=False,
    batch_size=10_000,
    parameter_transforms_to_save=param_transforms_to_save,
    cat_type="photometry",
    **supp_params,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[16], line 1
----> 1 basis.create_mock_library(
      2     emission_model_key=emission_key,
      3     out_name="library_Euclid_test",
      4     out_dir="./",
      5     overwrite=True,
      6     n_proc=1,
      7     verbose=False,
      8     batch_size=10_000,
      9     parameter_transforms_to_save=param_transforms_to_save,
     10     cat_type="photometry",
     11     **supp_params,
     12 )

NameError: name 'basis' is not defined

Spectroscopic Grid Creation

We can also build a library of spectra, which we can train a model from using an embedding network. At the default setting, the spectra would be at the wavelength range and resolution of our SPS grid, which is likely higher than required. We can change the wavelength array on our instrument or grid to a more reasonable choice.

from unyt import Angstrom

from synference import generate_constant_R

new_lam = generate_constant_R(300, start=100 * Angstrom, stop=100_000 * Angstrom)
[17]:
basis.create_mock_library(
    emission_model_key=emission_key,
    out_name="spectral_library_Euclid_test",
    out_dir="./",
    overwrite=True,
    n_proc=4,
    verbose=False,
    batch_size=10_000,
    parameter_transforms_to_save=param_transforms_to_save,
    cat_type="spectra",
    **supp_params,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[17], line 1
----> 1 basis.create_mock_library(
      2     emission_model_key=emission_key,
      3     out_name="spectral_library_Euclid_test",
      4     out_dir="./",
      5     overwrite=True,
      6     n_proc=4,
      7     verbose=False,
      8     batch_size=10_000,
      9     parameter_transforms_to_save=param_transforms_to_save,
     10     cat_type="spectra",
     11     **supp_params,
     12 )

NameError: name 'basis' is not defined