Convolution of binned data
The convolution of binned data is somewhat different than that of event-based, as the data has a somewhat different meaning. With binned data, the value in the columns indicates the center-value of the bin (usually). As such, we need to inform the convolution about the bin edges, which internally allows us to select the matching indices for a particular bin. Moreover, because each bin spans a range of delay times, as opposed to the singular value for event-based data, these bins can fall in more than one starformation bin. We use the bin edge info to calculate, given a particular target convolution time, the overlap fraction with the star formation bins (see below).
Apart from some different configuration in the convolution_instructions, binned data is convolved in the same way as event-based data.
One example of binned data is the Ensemble data type generated by binary_c.
Ensemble-based data is stored as nested dictionaries. An example of ensemble-based data looks like this:
"Xyield": {
"time": {
"-0.1": {
"source": {
"Wind": {
"isotope": {
"Al27": 1.3202421292393783e-08,
"Ar36": 1.7624018781546946e-08,
"Ar38": 3.502033439864038e-09,
"Ar40": 5.758546201573555e-12,
"B10": 2.4295643555965993e-13,
"B11": 1.0109986571758494e-12,
"Be9": 3.7843822119497306e-14,
[...]
}
}
[...]
}
}
[...]
}
}
With binary_c-python this type of data can be generated through the options explained in the ensemble-data logging notebook.
To use this type of data, however, one must first transform it to a different shape. In particular one must inflate the ensemble, turning it from a nested dictionairy to a rectangular data format. How to do so is covered in XXX (TODO: refer to notebook).
Because the binned data spans a range of delay times, we can not just use the SFR of one of the SFR bins, but rather we should calculate the average SFR by calculating how much the delay time bin overlaps with SFR bins, what the rates and bin-widths are, and then calculate a weighted average. In its most general form this is written as
Here \(R_{j}^{d}\) is the averaged starformation rate over for delay-time bin \(j\), \(w_i\) indicates the fraction that delay-time data bin \(j\) overlaps with SFR bin \(i\), \(R_{i}^{s}\) indicates the SFR rate in SFR bin \(i\), \(\Delta t_{i}^{s}\) denotes the bin width of SFR bin \(i\), and \(\Delta t^{d}\) indicates the cropped width of delay-time bin \(j\),
This cropped width is relevant when the delay-time bin exceeds beyond the last star formation rate bin edge.
based on the example below, can be expressed as follows
More generally, this can be written as,
TODO: describe how this is done internally with fraction of bin overlap
[1]:
import copy
import os
import h5py
import pandas as pd
import numpy as np
import json
import logging
import astropy.units as u
from syntheticstellarpopconvolve import convolve, default_convolution_config, default_convolution_instruction
from syntheticstellarpopconvolve.general_functions import calculate_bin_edges, temp_dir, generate_boilerplate_outputfile
#
TMP_DIR = temp_dir("notebook", "tutorial_binned_data_convolution", clean_path=True)
Loading up binned data
Lets start with setting up some dummy data. Here we use a bin-width of 1, and a bin-center of 0.25+N.
Based on the bin-centers we can automatically figure out what the edges are with the calculate_bin_edges function.
[2]:
# set up some data
records = [
{"time": 0.25, "value": 10, "probability": 1},
{"time": 1.25, "probability": 2, "value": 20},
{"time": 2.25, "probability": 3, "value": 30},
{"time": 0.25, "value": 11, "probability": 1.1},
{"time": 3.25, "probability": 4, "value": 40},
]
# load into dataframe
example_dataframe = pd.DataFrame.from_records(records)
sorted_unique_time_centers = np.sort(example_dataframe["time"].unique())
time_bin_edges = calculate_bin_edges(sorted_unique_time_centers)
# create file
output_hdf5_filename = os.path.join(TMP_DIR, "output_hdf5.h5")
generate_boilerplate_outputfile(output_hdf5_filename)
# store the data frame in the hdf5file
example_dataframe.to_hdf(output_hdf5_filename, key="input_data/binned_example")
We now need to set up the input file and set up the configuration, the starformation rate dictionary and the convolution steps. These steps are all unchanged with respect to non-binned data.
[3]:
# general configuration
convolution_config = copy.copy(default_convolution_config)
convolution_config["output_filename"] = output_hdf5_filename
convolution_config["tmp_dir"] = TMP_DIR
# set convolution time bins
convolution_config["convolution_lookback_time_bin_edges"] = np.arange(0, 3, 1) * u.yr
# set up SFR dict
convolution_config["SFR_info"] = {
'lookback_time_bin_edges': np.arange(0, 10, 1) * u.yr,
'starformation_rate_array': np.arange(0, 9) ** 2 * u.Msun / u.yr
}
To properly configure the code to convolve binned data, we need to indicate so in the convolution_instructions with the entry
convolution_config["convolution_instructions"] = [
{
...
"contains_binned_data": True,
...
}
]
and provide information about the delay time bins through
convolution_config["convolution_instructions"] = [
{
...
"delay_time_data_bin_info_dict": {
"delay_time_data_bin_edges": time_bin_edges * u.yr
},
...
}
]
[4]:
###
# convolution instructions
convolution_config["convolution_instructions"] = [
{
**default_convolution_instruction,
"convolution_type": "integrate",
"input_data_name": "binned_example",
"output_data_name": "binned_example",
"contains_binned_data": True,
"delay_time_data_bin_info_dict": {
"delay_time_data_bin_edges": time_bin_edges * u.yr
},
"data_column_dict": {
"normalized_yield": {"column_name": "probability", "unit": 1 / u.Msun},
"delay_time": {"column_name": "time", "unit": u.yr},
},
},
]
After this, the convolution handled in the same way as non-binned/event-based data.
[5]:
# convolve
convolve(config=convolution_config)
print("finished convolution")
# read out content and integrate until today
with h5py.File(convolution_config["output_filename"], "r") as output_hdf5_file:
#
main_group = "output_data/binned_example/binned_example/convolution_results/"
# loop over the formation-time bins
formation_time_bin_keys = list(output_hdf5_file[main_group].keys())
formation_time_bin_keys = sorted(
formation_time_bin_keys, key=lambda x: float(x.split(" ")[0])
)
for formation_time_bin_key in formation_time_bin_keys:
print("=================================")
print(f"formation_time_bin_key: {formation_time_bin_key}")
print("===============
###########
# Read out data
# convert units
unit_dict = json.loads(
output_hdf5_file[f"{main_group}/{formation_time_bin_key}"].attrs["units"]
)
unit_dict = {key: u.Unit(val) for key, val in unit_dict.items()}
print(unit_dict)
#
yield_result = output_hdf5_file[f"{main_group}/{formation_time_bin_key}/yield"][
()
]
print(yield_result)
Cell In[5], line 20
print("===============
^
SyntaxError: EOL while scanning string literal
Post-processing binned data
Post-processing of binned data is handled similar to the non-binned (event-based) data.
Convolution-by-sampling of binned data
Convolution-by-sampling of binned data is in essence similar to convolution-by-integration, but there are some subleties to consider:
filtering of data TODO: CHANGE
Assigning new values to initially binned quantities
Because the data in each bin actually represents a population of stars within that bin, one should keep in mind what the purpose of this convolution method is. The point of convolution-by-sampling is that one can generate instances of actual systems, which can then be assigned additional properties like positions, velocities, etc.
This method will generate N systems for each entry in the data (which consititutes a bin), but each of them will still have the same value. One idea is to assign new values to the systems, determined by some value randomly (or differently informed) chosen between the edges of the bin of that quantity.
for example, if \(10\) systems with \(\rm{log10Teff} = 3.0\) are sampled, and \(\rm{log10Teff}\) bins have size \(1.0\) and the relevant bin-edges are \(2.5\) and \(3.5\), then one can instead assign \(10\) random values between \(2.5\) and \(3.5\) to these \(10\) systems.
This is possible with the function sample_around_bin_center available in syntheticstellarpopconvolve.general_functions
def sample_around_bin_center(bin_edges, values):
"""
Basic function to handle sampling around bincenter given bin edges and values.
Note: this does not handle values that fall outside of the bins well
"""
bin_widths = np.diff(bin_edges)
indices = (
np.digitize(values, bin_edges) - 1
)
# get random values and scale
random_arr = np.random.random(indices.shape) - 0.5
random_arr = random_arr * bin_widths[indices]
# Add to values
sampled_values = values + random_arr
return sampled_values
Related to the previous section, if we want to assign formation lookback times, event lookback times, and filter out anything in the future, we can do that manually in the post-convolution function.
Below we present an example post-convolution function that handles these steps, as well as assigning values to the \(\rm{log10Teff}\) sampled accross their bins.
[19]:
from syntheticstellarpopconvolve.general_functions import sample_around_bin_center
def post_convolution_function(
config,
sfr_dict,
data_dict,
time_bin_info_dict,
convolution_results,
convolution_instruction
):
""" """
log10Teff_bin_edges = np.array([1.5, 2.5, 3.5])
convolution_results['sampled_log10Teff'] = sample_around_bin_center(
log10Teff_bin_edges,
data_dict['log10Teff'][convolution_results['sampled_indices']]
)
print('Sampled indices ', convolution_results['sampled_indices'])
print('Sampled log10Teff at bincenter ', data_dict['log10Teff'][convolution_results['sampled_indices']])
print('Sampled log10Teff values around bincenter: ', convolution_results['sampled_log10Teff'])
return convolution_results
Lets set up some data to convolve and generate samples from
[26]:
# set up data
records = [
{"time": 0.25, "log10Teff": 3., "probability": 10},
{"time": 1.25, "log10Teff": 2, "probability": 1},
]
#
example_dataframe = pd.DataFrame.from_records(records)
#
sorted_unique_time_centers = np.sort(example_dataframe["time"].unique())
time_bin_edges = calculate_bin_edges(sorted_unique_time_centers)
#
output_hdf5_filename = os.path.join(TMP_DIR, "sample_output_hdf5.h5")
generate_boilerplate_outputfile(output_hdf5_filename)
# store the data frame in the hdf5file
example_dataframe.to_hdf(output_hdf5_filename, key="input_data/binned_example_sample")
[27]:
# set up general config
convolution_config = copy.copy(default_convolution_config)
convolution_config["output_filename"] = output_hdf5_filename
convolution_config["tmp_dir"] = TMP_DIR
# set up convolution bins
convolution_config["convolution_lookback_time_bin_edges"] = np.arange(0, 3, 1) * u.yr
# set up SFR dict
convolution_config["SFR_info"] = {
'lookback_time_bin_edges': np.arange(0, 10, 1) * u.yr,
'starformation_rate_array': np.arange(0, 9) ** 2 * u.Msun / u.yr
}
# set up convolution instructions
convolution_config["convolution_instructions"] = [
{
**default_convolution_instruction,
"convolution_type": "sample",
"input_data_name": "binned_example_sample",
"output_data_name": "binned_example_sample",
"contains_binned_data": True,
"delay_time_data_bin_info_dict": {
"delay_time_data_bin_edges": time_bin_edges * u.yr
},
"data_column_dict": {
"normalized_yield": {"column_name": "probability", "unit": 1 / u.Msun},
"delay_time": {"column_name": "time", "unit": u.yr},
'log10Teff': "log10Teff",
},
'post_convolution_function': post_convolution_function,
'multiply_by_sfr_time_binsize': True
},
]
# convolve
convolve(config=convolution_config)
print("finished convolution")
{'column_name': 'probability', 'unit': <Quantity 1. 1 / solMass>}
{'column_name': 'time', 'unit': Unit("yr")}
Sampled indices [0 0 1]
Sampled log10Teff at bincenter [3. 3. 2.]
Sampled log10Teff values around bincenter: [3.05248434 2.87736634 1.60710813]
Sampled indices [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
Sampled log10Teff at bincenter [3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 2. 2. 2. 2. 2.]
Sampled log10Teff values around bincenter: [3.41162005 3.15500885 3.41386025 2.76235703 2.93290253 2.947159
3.45823447 2.73657305 2.97467374 3.12150858 2.55090045 2.50465049
2.52890603 2.81312424 2.50703613 2.70744225 2.80814476 1.95864139
2.02432386 2.31008159 2.44907214 2.40623331]
finished convolution