Tutorial: convolution on the fly
On-the-fly convolution is a method to run population-synthesis simulations during the convolution. SSPC provides the total mass formed in stars at a particular time, and optionally a metallicity distribution.
This method of convolution is useful if the population-synthesis simulations are in some way affected by the environment.
This notebook will provide a boilerplate setup, but the user should provide their own actual implementation in the ‘on-the-fly’ function by e.g. using binary_c-python to evolve a population.
This method does not require any input data to be present, as it will generate that as the convolution runs.
[1]:
import os
import copy
import numpy as np
import astropy.units as u
from syntheticstellarpopconvolve import convolve, default_convolution_config, default_convolution_instruction
from syntheticstellarpopconvolve.general_functions import temp_dir, generate_boilerplate_outputfile
#
TMP_DIR = temp_dir("notebook", "tutorial_on_the_fly_convolution", clean_path=True)
To run an on-the-fly convolution, the user should configure the convolution_instruction with
convolution_config["convolution_instructions"] = [
{
...
"convolution_type": "on-the-fly",
...
"on_the_fly_function": <callable function with appropriate arguments>,
...
"convolution_direction": "forward"
}
]
where the ‘on_the_fly_function’ calleable should look like
def dummy_on_the_fly_function(total_star_formation_in_bin):
"""
on-the-fly function that should work
"""
return {}
when no metallicity information is required, or
def dummy_on_the_fly_function(total_star_formation_in_bin, metallicity_distribution):
"""
on-the-fly function that should work
"""
return {}
if it is.
The entire set of arguments that the on-the-fly function has access to is:
available_parameters = {
# Standard info
"config": config,
"sfr_dict": sfr_dict,
"time_bin_info_dict": time_bin_info_dict,
"convolution_instruction": convolution_instruction,
# specific required info
"total_star_formation_in_bin": total_star_formation_in_bin,
"metallicity_distribution": metallicity_distribution,
**convolution_instruction.get("on_the_fly_function_extra_parameters", {}),
}
[2]:
def dummy_on_the_fly_function(total_star_formation_in_bin):
"""
on-the-fly function
"""
print(total_star_formation_in_bin)
return {}
Lets put together a small working example to run on-the-fly convolution.
[7]:
# setup output file
output_hdf5_filename = os.path.join(TMP_DIR, "output_hdf5.h5")
generate_boilerplate_outputfile(output_hdf5_filename)
# general config, SFR, convolution bin
convolution_config = copy.copy(default_convolution_config)
convolution_config["output_filename"] = output_hdf5_filename
convolution_config["SFR_info"] = {
"lookback_time_bin_edges": np.array([0, 1, 2, 3, 4, 5]) * u.yr,
"starformation_rate_array": np.array([1, 1, 1, 1, 1]) * u.Msun/u.yr,
}
# convolution_config["convolution_lookback_time_bin_edges"] = (
# np.array([0, 1, 2, 3, 4]) * u.yr
# )
convolution_config["convolution_instructions"] = [
{
**default_convolution_instruction,
"output_data_name": "BHBH",
"convolution_type": "on-the-fly",
"convolution_direction": "forward",
"on_the_fly_function": dummy_on_the_fly_function,
},
]
convolution_config["tmp_dir"] = os.path.join(TMP_DIR, "tmp")
#
convolve(convolution_config)
/home/david/projects/binary_c_root/syntheticstellarpopconvolve/syntheticstellarpopconvolve/convolve_populations.py:360: UserWarning: On-the-fly convolution is currently not fully tested
warnings.warn("On-the-fly convolution is currently not fully tested")
1.0 solMass
1.0 solMass
1.0 solMass
1.0 solMass
1.0 solMass
example implementation: binary_c-python
In the following section we will flesh out an actual setup that uses binary_c-python to evolve some systems.
Note: the monte-carlo module in binary_c-python is not yet fully finished, so we’ll make use of a custom generator setup.
[ ]:
import binarycpython
import numpy as np
from scipy.interpolate import interp1d
from scipy.integrate import cumulative_trapezoid
class KroupaIMFSampler:
def __init__(self, m_min=0.08, m_max=150.0, resolution=10000):
self.m_min = m_min
self.m_max = m_max
self.masses = np.linspace(m_min, m_max, resolution)
self.pdf_values = self._kroupa_pdf(self.masses)
self.pdf_values /= np.trapz(self.pdf_values, self.masses) # normalize
# Compute the CDF
cdf = cumulative_trapezoid(self.pdf_values, self.masses, initial=0)
cdf /= cdf[-1] # normalize to [0, 1]
# Inverse CDF (interpolator)
self.inv_cdf_interp = interp1d(cdf, self.masses, bounds_error=False, fill_value=(m_min, m_max))
def _kroupa_pdf(self, m):
"""
Unnormalized Kroupa IMF.
"""
alpha1, alpha2 = 1.3, 2.3
m_break = 0.5
pdf = np.where(m < m_break,
m**(-alpha1),
(m_break**(alpha2 - alpha1)) * m**(-alpha2))
return pdf
def sample(self, size=1):
u = np.random.uniform(0, 1, size)
return self.inv_cdf_interp(u)
def generate_population(self, total_mass_limit=1000.0):
"""
Generator that yields (star_id, mass) pairs until the total mass limit is exceeded.
"""
current_mass = 0.0
i = 1
while current_mass < total_mass_limit:
mass = self.sample(1)[0]
if current_mass + mass > total_mass_limit:
break
system = {'M1': mass}
yield system
current_mass += mass
i += 1
def binarycpython_on_the_fly_function(total_star_formation_in_bin):
"""
Example on-the-fly function that makes use
"""
# Load system generator
sampler = KroupaIMFSampler()
system_gen_obj = sampler.generate_population(total_mass_limit=10)
# Set up population
test_pop = Population(verbosity=1)
# Load settings
test_pop.set(**settings)
test_pop.set(num_cores=16)
# Add situational settings
test_pop.set(
C_logging_code=custom_logging_string,
parse_function=parse_function,
evolution_type="custom_generator",
custom_generator=system_gen_obj,
)
test_pop.set(data_dir=output_dir, base_filename="output_sn.dat")
# # Get version info of the binary_c build
# version_info = test_pop.return_binary_c_version_info(parsed=True)
# create local tmp_dir
test_pop.set(tmp_dir=os.path.join(test_pop.custom_options["data_dir"], "local_tmp_dir"))
test_pop.set(log_args_dir=test_pop.population_options["tmp_dir"])
os.makedirs(test_pop.population_options["tmp_dir"], exist_ok=True)
test_pop.evolve()
# read out data and return
return {}
###################
###################
#
# setup output file
output_hdf5_filename = os.path.join(TMP_DIR, "bcp_output_hdf5.h5")
generate_boilerplate_outputfile(output_hdf5_filename)
# general config, SFR, convolution bin
convolution_config = copy.copy(default_convolution_config)
convolution_config["output_filename"] = output_hdf5_filename
convolution_config["SFR_info"] = {
"lookback_time_bin_edges": np.array([0, 1, 2, 3, 4, 5]) * u.yr,
"starformation_rate_array": np.array([1, 1, 1, 1, 1]) * u.Msun/u.yr,
}
# convolution_config["convolution_lookback_time_bin_edges"] = (
# np.array([0, 1, 2, 3, 4]) * u.yr
# )
convolution_config["convolution_instructions"] = [
{
**default_convolution_instruction,
"output_data_name": "bcp_on_the_fly",
"convolution_type": "on-the-fly",
"convolution_direction": "forward",
"on_the_fly_function": binarycpython_on_the_fly_function,
},
]
convolution_config["tmp_dir"] = os.path.join(TMP_DIR, "tmp")
#
convolve(convolution_config)