"""
Module containing the functions to handle the Moe&Distrefano data
This class object is an extension to the population grid object
"""
# pylint: disable=E1101
import os
import gc
import copy
import json
import py_rinterpolate
from binarycpython.utils.functions import (
verbose_print,
)
from binarycpython.utils.dicts import update_dicts, normalize_dict
from binarycpython.utils.population_extensions.distribution_functions import (
Moecache,
LOG_LN_CONVERTER,
)
from binarycpython.utils.population_extensions.grid_options_defaults import (
_MOE2017_VERBOSITY_LEVEL,
)
from binarycpython.utils import moe_di_stefano_2017_data
[docs]class Moe_di_Stefano_2017:
"""
Extension to the population grid object that contains functionality to handle handle the Moe & distefano distributions
"""
def __init__(self, **kwargs):
"""
Init function for the Moe_di_Stefano_2017 class
"""
return
[docs] def set_moe_di_stefano_settings(self, options=None):
"""
Function to set user input configurations for the Moe & di Stefano methods
If nothing is passed then we just use the default options
"""
if not options:
options = {}
# Take the option dictionary that was given and override.
options = update_dicts(self.grid_options["Moe2017_options"], options)
self.grid_options["Moe2017_options"] = copy.deepcopy(options)
# Write options to a file
os.makedirs(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
"moeopts.dat",
),
"w",
encoding="utf-8",
) as f:
f.write(
json.dumps(
self.grid_options["Moe2017_options"], indent=4, ensure_ascii=False
)
)
f.close()
def _load_moe_di_stefano_data(self):
"""
Function to load the moe & di stefano data
"""
# Only if the grid is loaded and Moecache contains information
if not self.grid_options["_loaded_Moe2017_data"]: # and not Moecache:
if self.grid_options["_Moe2017_JSON_data"]:
# Use the existing (perhaps modified) JSON data
json_data = self.grid_options["_Moe2017_JSON_data"]
else:
# Load the JSON data from a file
json_data = self.get_moe_di_stefano_dataset(
self.grid_options["Moe2017_options"],
verbosity=self.grid_options["verbosity"],
)
# entry of log10M1 is a list containing 1 dict.
# We can take the dict out of the list
if isinstance(json_data["log10M1"], list):
json_data["log10M1"] = json_data["log10M1"][0]
# save this data in case we want to modify it later
self.grid_options["_Moe2017_JSON_data"] = json_data
# Get all the masses
logmasses = sorted(json_data["log10M1"].keys())
if not logmasses:
msg = "The table does not contain masses."
verbose_print(
"\tMoe_di_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
0,
)
raise ValueError(msg)
# Write to file
os.makedirs(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
"moe.log",
),
"w",
encoding="utf-8",
) as logfile:
logfile.write("log₁₀Masses(M☉) {}\n".format(logmasses))
# Get all the periods and see if they are all consistently present
logperiods = []
for logmass in logmasses:
if not logperiods:
logperiods = sorted(json_data["log10M1"][logmass]["logP"].keys())
dlog10P = float(logperiods[1]) - float(logperiods[0])
current_logperiods = sorted(json_data["log10M1"][logmass]["logP"])
if logperiods != current_logperiods:
msg = (
"Period values are not consistent throughout the dataset logperiods = "
+ " ".join(str(x) for x in logperiods)
+ "\nCurrent periods = "
+ " ".join(str(x) for x in current_logperiods)
)
verbose_print(
"\tMoe_di_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
0,
)
raise ValueError(msg)
############################################################
# log10period binwidth : of course this assumes a fixed
# binwidth, so we check for this too.
for i in range(len(current_logperiods) - 1):
if not dlog10P == (
float(current_logperiods[i + 1]) - float(current_logperiods[i])
):
msg = "Period spacing is not consistent throughout the dataset"
verbose_print(
"\tMoe_di_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
0,
)
raise ValueError(msg)
# save the logperiods list in the cache:
# this is used in the renormalization integration
Moecache["logperiods"] = logperiods
# Write to file
os.makedirs(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano", "moe.log"),
"a",
encoding="utf-8",
) as logfile:
logfile.write("log₁₀Periods(days) {}\n".format(logperiods))
# Fill the global dict
for logmass in logmasses:
# Create the multiplicity table
if not Moecache.get("multiplicity_table", None):
Moecache["multiplicity_table"] = []
# multiplicity as a function of primary mass
Moecache["multiplicity_table"].append(
[
float(logmass),
json_data["log10M1"][logmass]["f_multi"],
json_data["log10M1"][logmass]["single star fraction"],
json_data["log10M1"][logmass]["binary star fraction"],
json_data["log10M1"][logmass]["triple/quad star fraction"],
]
)
############################################################
# a small log10period which we can shift just outside the
# table to force integration out there to zero
epslog10P = 1e-8 * dlog10P
############################################################
# loop over either binary or triple-outer periods
first = 1
# Go over the periods
for logperiod in logperiods:
############################################################
# distributions of binary and triple star fractions
# as a function of mass, period.
#
# Note: these should be per unit log10P, hence we
# divide by dlog10P
if first:
first = 0
# Create the multiplicity table
if not Moecache.get("period_distributions", None):
Moecache["period_distributions"] = []
############################################################
# lower bound the period distributions to zero probability
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod) - 0.5 * dlog10P - epslog10P,
0.0,
0.0,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod) - 0.5 * dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod),
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
############################################################
# distributions as a function of mass, period, q
#
# First, get a list of the qs given by Moe
#
qs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["q"])
# Fill the data and 'normalise'
qdata = self.fill_data(
qs, json_data["log10M1"][logmass]["logP"][logperiod]["q"]
)
# Create the multiplicity table
if not Moecache.get("q_distributions", None):
Moecache["q_distributions"] = []
for q in qs:
Moecache["q_distributions"].append(
[float(logmass), float(logperiod), float(q), qdata[q]]
)
############################################################
# eccentricity distributions as a function of mass, period, ecc
eccs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["e"])
# Fill the data and 'normalise'
ecc_data = self.fill_data(
eccs, json_data["log10M1"][logmass]["logP"][logperiod]["e"]
)
# Create the multiplicity table
if not Moecache.get("ecc_distributions", None):
Moecache["ecc_distributions"] = []
for ecc in eccs:
Moecache["ecc_distributions"].append(
[
float(logmass),
float(logperiod),
float(ecc),
ecc_data[ecc],
]
)
############################################################
# upper bound the period distributions to zero probability
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiods[-1])
+ 0.5 * dlog10P, # TODO: why this shift? to center it?
json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiods[-1]) + 0.5 * dlog10P + epslog10P,
0.0,
0.0,
]
)
verbose_print(
"\tMoe_di_Stefano_2017: Length period_distributions table: {}".format(
len(Moecache["period_distributions"])
),
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
verbose_print(
"\tMoe_di_Stefano_2017: Length multiplicity table: {}".format(
len(Moecache["multiplicity_table"])
),
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
verbose_print(
"\tMoe_di_Stefano_2017: Length q table: {}".format(
len(Moecache["q_distributions"])
),
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
verbose_print(
"\tMoe_di_Stefano_2017: Length ecc table: {}".format(
len(Moecache["ecc_distributions"])
),
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
# Write to log file
os.makedirs(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
"moecache.json",
),
"w",
encoding="utf-8",
) as cache_filehandle:
cache_filehandle.write(
json.dumps(Moecache, indent=4, ensure_ascii=False)
)
# Signal that the data has been loaded
self.grid_options["_loaded_Moe2017_data"] = True
def _set_moe_di_stefano_distributions(self):
"""
Function to set the Moe & di Stefano distribution
"""
############################################################
# first, the multiplicity, this is 1,2,3,4, ...
# for singles, binaries, triples, quadruples, ...
max_multiplicity = self.get_max_multiplicity(
self.grid_options["Moe2017_options"]["multiplicity_modulator"]
)
verbose_print(
"\tMoe_di_Stefano_2017: Max multiplicity = {}".format(max_multiplicity),
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
######
# Setting up the grid variables
# Multiplicity
self.add_grid_variable(
name="multiplicity",
parameter_name="multiplicity",
longname="multiplicity",
valuerange=[1, max_multiplicity],
samplerfunc="self.const_int(1, {n}, {n})".format(n=max_multiplicity),
precode='self.grid_options["multiplicity"] = multiplicity; self.bse_options["multiplicity"] = multiplicity; options={}'.format(
self.grid_options["Moe2017_options"]
),
condition="({}[int(multiplicity)-1] > 0)".format(
str(self.grid_options["Moe2017_options"]["multiplicity_modulator"])
),
gridtype="discrete",
probdist=1,
)
############################################################
# always require M1, for all systems
#
# log-spaced m1 with given resolution
self.add_grid_variable(
name="lnM_1",
parameter_name="M_1",
longname="Primary mass",
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"]["M"][0]
or "self.const_linear(np.log({}), np.log({}), {})".format(
self.grid_options["Moe2017_options"]["ranges"]["M"][0],
self.grid_options["Moe2017_options"]["ranges"]["M"][1],
self.grid_options["Moe2017_options"]["resolutions"]["M"][0],
),
valuerange=[
"np.log({})".format(
self.grid_options["Moe2017_options"]["ranges"]["M"][0]
),
"np.log({})".format(
self.grid_options["Moe2017_options"]["ranges"]["M"][1]
),
],
gridtype="centred",
dphasevol="dlnM_1",
precode='M_1 = np.exp(lnM_1); options["M_1"]=M_1',
probdist="self.Moe_di_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
str(dict(self.grid_options["Moe2017_options"]))[1:-1],
"'multiplicity': multiplicity",
"'M_1': M_1",
),
)
# Go to higher multiplicities
if max_multiplicity >= 2:
# binaries: period
self.add_grid_variable(
name="log10per",
parameter_name="orbital_period",
longname="log10(Orbital_Period)",
probdist=1.0,
condition='(self.grid_options["multiplicity"] >= 2)',
branchpoint=1
if max_multiplicity > 1
else 0, # Signal here to put a branchpoint if we have a max multiplicity higher than 1.
gridtype="centred",
dphasevol="({} * dlog10per)".format(LOG_LN_CONVERTER),
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"][
"logP"
][0]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
self.grid_options["Moe2017_options"]["resolutions"]["logP"][0],
),
precode="""orbital_period = 10.0**log10per
qmin={}/M_1
qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
""".format(
self.grid_options["Moe2017_options"]["Mmin"]
),
) # TODO: change the maximum_mass_ratio_for_RLOF
# binaries: mass ratio
self.add_grid_variable(
name="q",
parameter_name="M_2",
longname="Mass ratio",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/M_1",
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "qmax",
],
probdist=1,
gridtype="centred",
dphasevol="dq",
precode="""
M_2 = q * M_1
sep = calc_sep_from_period(M_1, M_2, orbital_period)
""",
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"]["M"][1]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", [None, None])[0]
else "{}/M_1".format(self.grid_options["Moe2017_options"]["Mmin"]),
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", [None, None])[1]
else "qmax",
self.grid_options["Moe2017_options"]["resolutions"]["M"][1],
),
)
# (optional) binaries: eccentricity
if self.grid_options["Moe2017_options"]["resolutions"]["ecc"][0] > 0:
self.add_grid_variable(
name="ecc",
parameter_name="eccentricity",
longname="Eccentricity",
probdist=1,
gridtype="centred",
dphasevol="decc",
precode="eccentricity=ecc",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][1],
],
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"][
"ecc"
][0]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][1],
self.grid_options["Moe2017_options"]["resolutions"]["ecc"][0],
),
)
# Now for triples and quadruples
if max_multiplicity >= 3:
# Triple: period
self.add_grid_variable(
name="log10per2",
parameter_name="orbital_period_triple",
longname="log10(Orbital_Period2)",
probdist=1.0,
condition='(self.grid_options["multiplicity"] >= 3)',
branchpoint=2
if max_multiplicity > 2
else 0, # Signal here to put a branchpoint if we have a max multiplicity higher than 1.
gridtype="centred",
dphasevol="({} * dlog10per2)".format(LOG_LN_CONVERTER),
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"][
"logP"
][1]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
self.grid_options["Moe2017_options"]["resolutions"]["logP"][1],
),
precode="""orbital_period_triple = 10.0**log10per2
q2min={}/(M_1+M_2)
q2max=maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
""".format(
self.grid_options["Moe2017_options"]["Mmin"]
),
)
# Triples: mass ratio
# Note, the mass ratio is M_outer/M_inner
self.add_grid_variable(
name="q2",
parameter_name="M_3",
longname="Mass ratio outer/inner",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_1+M_2)",
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q2max",
],
probdist=1,
gridtype="centred",
dphasevol="dq2",
precode="""
M_3 = q2 * (M_1 + M_2)
sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
eccentricity2=0
""",
samplerfunc=self.grid_options["Moe2017_options"]["samplerfuncs"][
"M"
][2]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_1+M_2)",
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q2max",
self.grid_options["Moe2017_options"]["resolutions"]["M"][2],
),
)
# (optional) triples: eccentricity
if self.grid_options["Moe2017_options"]["resolutions"]["ecc"][1] > 0:
self.add_grid_variable(
name="ecc2",
parameter_name="eccentricity2",
longname="Eccentricity of the triple",
probdist=1,
gridtype="centred",
dphasevol="decc2",
precode="eccentricity2=ecc2",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][1],
],
samplerfunc=self.grid_options["Moe2017_options"][
"samplerfuncs"
]["ecc"][1]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][1],
self.grid_options["Moe2017_options"]["resolutions"]["ecc"][
1
],
),
)
if max_multiplicity == 4:
# Quadruple: period
self.add_grid_variable(
name="log10per3",
parameter_name="orbital_period_quadruple",
longname="log10(Orbital_Period3)",
probdist=1.0,
condition='(self.grid_options["multiplicity"] >= 4)',
branchpoint=3
if max_multiplicity > 3
else 0, # Signal here to put a branchpoint if we have a max multiplicity higher than 1.
gridtype="centred",
dphasevol="({} * dlog10per3)".format(LOG_LN_CONVERTER),
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.grid_options["Moe2017_options"][
"samplerfuncs"
]["logP"][2]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["logP"][0],
self.grid_options["Moe2017_options"]["ranges"]["logP"][1],
self.grid_options["Moe2017_options"]["resolutions"]["logP"][
2
],
),
precode="""orbital_period_quadruple = 10.0**log10per3
q3min={}/(M_3)
q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
""".format(
self.grid_options["Moe2017_options"]["Mmin"]
),
)
# Quadruple: mass ratio : M_outer / M_inner
self.add_grid_variable(
name="q3",
parameter_name="M_4",
longname="Mass ratio outer low/outer high",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_3)",
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q3max",
],
probdist=1,
gridtype="centred",
dphasevol="dq3",
precode="""
M_4 = q3 * M_3
sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
eccentricity3=0
""",
samplerfunc=self.grid_options["Moe2017_options"][
"samplerfuncs"
]["M"][3]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["q"][0]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_3)",
self.grid_options["Moe2017_options"]["ranges"]["q"][1]
if self.grid_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q3max",
self.grid_options["Moe2017_options"]["resolutions"]["M"][2],
),
)
# (optional) triples: eccentricity
if (
self.grid_options["Moe2017_options"]["resolutions"]["ecc"][2]
> 0
):
self.add_grid_variable(
name="ecc3",
parameter_name="eccentricity3",
longname="Eccentricity of the triple+quadruple/outer binary",
probdist=1,
gridtype="centred",
dphasevol="decc3",
precode="eccentricity3=ecc3",
valuerange=[
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
1
],
],
samplerfunc=self.grid_options["Moe2017_options"][
"samplerfuncs"
]["ecc"][2]
or "self.const_linear({}, {}, {})".format(
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.grid_options["Moe2017_options"]["ranges"]["ecc"][
1
],
self.grid_options["Moe2017_options"]["resolutions"][
"ecc"
][2],
),
)
# Now we are at the last part.
# Here we should combine all the information that we calculate and update the options
# dictionary. This will then be passed to the Moe_di_Stefano_2017_pdf to calculate
# the real probability. The trick we use is to strip the options_dict as a string
# and add some keys to it:
updated_options = "{{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}".format(
str(dict(self.grid_options["Moe2017_options"]))[1:-1],
'"multiplicity": multiplicity',
'"M_1": M_1',
'"M_2": M_2',
'"M_3": M_3',
'"M_4": M_4',
'"P": orbital_period',
'"P2": orbital_period_triple',
'"P3": orbital_period_quadruple',
'"ecc": eccentricity',
'"ecc2": eccentricity2',
'"ecc3": eccentricity3',
)
probdist_addition = "self.Moe_di_Stefano_2017_pdf({}, verbosity=self.grid_options['verbosity'])['total_probdens']".format(
updated_options
)
# and finally the probability calculator
self.grid_options["_grid_variables"][self._last_grid_variable()][
"probdist"
] = probdist_addition
verbose_print(
"\tMoe_di_Stefano_2017: Added final call to the pdf function",
self.grid_options["verbosity"],
_MOE2017_VERBOSITY_LEVEL,
)
# Signal that the MOE2017 grid has been set
self.grid_options["_set_Moe2017_grid"] = True
################################################################################################
[docs] def Moe_di_Stefano_2017(self, options=None):
"""
Function to handle setting the user input settings,
set up the data and load that into interpolators and
then set the distribution functions
Takes a dictionary as its only argument
"""
default_options = {
"apply settings": True,
"setup grid": True,
"load data": True,
"clean cache": False,
"clean load flag": False,
"clean all": False,
}
if not options:
options = {}
options = update_dicts(default_options, options)
# clean cache?
if options["clean all"] or options["clean cache"]:
Moecache.clear()
if options["clean all"] or options["clean load flag"]:
self.grid_options["_loaded_Moe2017_data"] = False
# Set the user input
if options["apply settings"]:
self.set_moe_di_stefano_settings(options=options)
# Load the data
if options["load data"]:
self._load_moe_di_stefano_data()
# construct the grid here
if options["setup grid"]:
self._set_moe_di_stefano_distributions()
def _clean_interpolators(self):
"""
Function to clean up the interpolators after a run
We look in the Moecache global variable for items that are interpolators.
Should be called by the general cleanup function AND the thread cleanup function
"""
interpolator_keys = []
for key, value in Moecache.items():
if isinstance(value, py_rinterpolate.Rinterpolate):
interpolator_keys.append(key)
for key in interpolator_keys:
Moecache[key].destroy()
del Moecache[key]
gc.collect()
def _calculate_multiplicity_fraction(self, system_dict):
"""
Function to calculate multiplicity fraction
Makes use of the self.bse_options['multiplicity'] value. If its not set, it will raise an error
grid_options['multiplicity_fraction_function'] will be checked for the choice
TODO: add option to put a manual binary fraction in here (solve via negative numbers being the functions)
"""
# Just return 1 if no option has been chosen
if self.grid_options["multiplicity_fraction_function"] in [0, "None"]:
verbose_print(
"_calculate_multiplicity_fraction: Chosen not to use any multiplicity fraction.",
self.grid_options["verbosity"],
3,
)
return 1
# Raise an error if the multiplicity is not set
if not system_dict.get("multiplicity", None):
msg = "Multiplicity value has not been set. When using a specific multiplicity fraction function please set the multiplicity"
raise ValueError(msg)
# Go over the chosen options
if self.grid_options["multiplicity_fraction_function"] in [1, "Arenou2010"]:
# Arenou 2010 will be used
verbose_print(
"_calculate_multiplicity_fraction: Using Arenou 2010 to calculate multiplicity fractions",
self.grid_options["verbosity"],
3,
)
binary_fraction = self.Arenou2010_binary_fraction(system_dict["M_1"])
multiplicity_fraction_dict = {
1: 1 - binary_fraction,
2: binary_fraction,
3: 0,
4: 0,
}
elif self.grid_options["multiplicity_fraction_function"] in [2, "Raghavan2010"]:
# Raghavan 2010 will be used
verbose_print(
"_calculate_multiplicity_fraction: Using Raghavan (2010) to calculate multiplicity fractions",
self.grid_options["verbosity"],
3,
)
binary_fraction = self.raghavan2010_binary_fraction(system_dict["M_1"])
multiplicity_fraction_dict = {
1: 1 - binary_fraction,
2: binary_fraction,
3: 0,
4: 0,
}
elif self.grid_options["multiplicity_fraction_function"] in [3, "Moe2017"]:
# We need to check several things now here:
# First, are the options for the MOE2017 grid set? On start it is filled with the default settings
if not self.grid_options["Moe2017_options"]:
msg = "The MOE2017 options do not seem to be set properly. The value is {}".format(
self.grid_options["Moe2017_options"]
)
raise ValueError(msg)
# Second: is the Moecache filled.
if not Moecache:
verbose_print(
"_calculate_multiplicity_fraction: Moecache is empty. It needs to be filled with the data for the interpolators. Loading the data now",
self.grid_options["verbosity"],
3,
)
# Load the data
self._load_moe_di_stefano_data()
# record the prev value
prev_M1_value_ms = self.grid_options["Moe2017_options"].get("M_1", None)
# Set value of M1 of the current system
self.grid_options["Moe2017_options"]["M_1"] = system_dict["M_1"]
# Calculate the multiplicity fraction
multiplicity_fraction_list = (
self.Moe_di_Stefano_2017_multiplicity_fractions(
self.grid_options["Moe2017_options"], self.grid_options["verbosity"]
)
)
# Turn into dict
multiplicity_fraction_dict = {
el + 1: multiplicity_fraction_list[el]
for el in range(len(multiplicity_fraction_list))
}
# Set the prev value back
self.grid_options["Moe2017_options"]["M_1"] = prev_M1_value_ms
# we don't know what to do next
else:
msg = "Chosen value for the multiplicity fraction function is not known."
raise ValueError(msg)
# To make sure we normalize the dictionary
multiplicity_fraction_dict = normalize_dict(
multiplicity_fraction_dict, verbosity=self.grid_options["verbosity"]
)
verbose_print(
"Multiplicity: {} multiplicity_fraction: {}".format(
system_dict["multiplicity"],
multiplicity_fraction_dict[system_dict["multiplicity"]],
),
self.grid_options["verbosity"],
3,
)
return multiplicity_fraction_dict[system_dict["multiplicity"]]
[docs] def get_moe_di_stefano_dataset(self, options, verbosity=0):
"""
Function to get the default Moe and di Stefano dataset or accept a user input.
Returns a dict containing the (JSON) data.
"""
json_data = None
if "JSON" in options:
# use the JSON data passed in
json_data = options["JSON"]
elif "file" in options:
# use the file passed in, if provided
if not os.path.isfile(options["file"]):
verbose_print(
"The provided 'file' Moe and de Stefano JSON file does not seem to exist at {}".format(
options["file"]
),
verbosity,
1,
)
raise ValueError
if not options["file"].endswith(".json"):
verbose_print(
"Provided filename is not a json file",
verbosity,
1,
)
else:
# Read input data and Clean up the data if there are white spaces around the keys
with open(options["file"], "r", encoding="utf-8") as data_filehandle:
datafile_data = data_filehandle.read()
datafile_data = datafile_data.replace('" ', '"')
datafile_data = datafile_data.replace(' "', '"')
datafile_data = datafile_data.replace(' "', '"')
json_data = json.loads(datafile_data)
if not json_data:
# no JSON data or filename given, use the default 2017 dataset
verbose_print(
"Using the default Moe and de Stefano 2017 datafile",
verbosity,
1,
)
json_data = copy.deepcopy(moe_di_stefano_2017_data.moe_di_stefano_2017_data)
return json_data
[docs] def get_Moe_di_Stefano_2017_default_options(self):
"""
Function that returns the default options for the Moe & diStefano grid options
"""
return {
# place holder for the JSON data to be used if a file
# isn't specified
"JSON": None,
# resolution data
"resolutions": {
"M": [
20, # M1
20, # M2 (i.e. q)
0, # M3 currently unused
0, # M4 currently unused
],
"logP": [
20, # P2 (binary period)
0, # P3 (triple period) currently unused
0, # P4 (quadruple period) currently unused
],
"ecc": [
10, # e (binary eccentricity)
0, # e2 (triple eccentricity) currently unused
0, # e3 (quadruple eccentricity) currently unused
],
},
"samplerfuncs": {
"M": [None, None, None, None],
"logP": [None, None, None],
"ecc": [None, None, None],
},
"ranges": {
# stellar masses (Msun)
"M": [
self.minimum_stellar_mass()
* 1.05, # 0.08 is a tad bit above the minimum mass. Don't sample at 0.07, otherwise the first row of q values will have a phasevol of 0. Anything higher is fine.
80.0, # (rather arbitrary) upper mass cutoff
],
"q": [
None, # artificial qmin : set to None to use default
None, # artificial qmax : set to None to use default
],
"logP": [0.0, 8.0], # 0 = log10(1 day) # 8 = log10(10^8 days)
"ecc": [0.0, 0.99],
},
# minimum stellar mass
"Mmin": self.minimum_stellar_mass(), # We take the value that binary_c has set as the default
# multiplicity model (as a function of log10M1)
#
# You can use 'Poisson' which uses the system multiplicity
# given by Moe and maps this to single/binary/triple/quad
# fractions.
#
# Alternatively, 'data' takes the fractions directly
# from the data, but then triples and quadruples are
# combined (and there are NO quadruples).
"multiplicity_model": "Poisson",
# multiplicity modulator:
# [single, binary, triple, quadruple]
#
# e.g. [1,0,0,0] for single stars only
# [0,1,0,0] for binary stars only
#
# defaults to [1,1,0,0] i.e. all types
#
"multiplicity_modulator": [
1, # single
1, # binary
0, # triple
0, # quadruple
],
# given a mix of multiplicities, you can either (noting that
# here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
#
# 'norm' : normalise so the whole population is 1.0
# after implementing the appropriate fractions
# S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
#
# 'raw' : stick to what is predicted, i.e.
# S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
# without normalisation
# (in which case the total probability < 1.0 unless
# all you use single, binary, triple and quadruple)
#
# 'merge' : e.g. if you only have single and binary,
# add the triples and quadruples to the binaries, so
# binaries represent all multiple systems
# ...
# *** this is canonical binary population synthesis ***
#
# Note: if multiplicity_modulator == [1,1,1,1] this
# option does nothing (equivalent to 'raw').
#
#
# note: if you only set one multiplicity_modulator
# to 1, and all the others to 0, then normalising
# will mean that you effectively have the same number
# of stars as single, binary, triple or quad (whichever
# is non-zero) i.e. the multiplicity fraction is ignored.
# This is probably not useful except for
# testing purposes or comparing to old grids.
"normalize_multiplicities": "merge",
# q extrapolation (below 0.15 and above 0.9) method. We can choose from ['flat', 'linear', 'plaw2', 'nolowq']
"q_low_extrapolation_method": "linear",
"q_high_extrapolation_method": "linear",
}
[docs] def get_Moe_di_Stefano_2017_default_options_description(self):
"""
Function to return the descriptions for all the Moe & diStefano grid options
"""
return {
"resolutions": "",
"ranges": "",
"Mmin": "Minimum stellar mass",
"multiplicity_model": """
multiplicity model (as a function of log10M1)
You can use 'Poisson' which uses the system multiplicity
given by Moe and maps this to single/binary/triple/quad
fractions.
Alternatively, 'data' takes the fractions directly
from the data, but then triples and quadruples are
combined (and there are NO quadruples).
""",
"multiplicity_modulator": """
[single, binary, triple, quadruple]
e.g. [1,0,0,0] for single stars only
[0,1,0,0] for binary stars only
defaults to [1,1,0,0] i.e. singles and binaries
""",
"normalize_multiplicities": """
'norm': normalise so the whole population is 1.0
after implementing the appropriate fractions
S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
given a mix of multiplicities, you can either (noting that
here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
note: if you only set one multiplicity_modulator
to 1, and all the others to 0, then normalising
will mean that you effectively have the same number
of stars as single, binary, triple or quad (whichever
is non-zero) i.e. the multiplicity fraction is ignored.
This is probably not useful except for
testing purposes or comparing to old grids.
'raw' : stick to what is predicted, i.e.
S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
without normalisation
(in which case the total probability < 1.0 unless
all you use single, binary, triple and quadruple)
'merge' : e.g. if you only have single and binary,
add the triples and quadruples to the binaries, so
binaries represent all multiple systems
...
*** this is canonical binary population synthesis ***
It only takes the maximum multiplicity into account,
i.e. it doesn't multiply the resulting array by the multiplicity modulator again.
This prevents the resulting array to always be 1 if only 1 multiplicity modulator element is nonzero
Note: if multiplicity_modulator == [1,1,1,1]. this option does nothing (equivalent to 'raw').
""",
"q_low_extrapolation_method": """
q extrapolation (below 0.15) method
none
flat
linear2
plaw2
nolowq
""",
"q_high_extrapolation_method": "Same as q_low_extrapolation_method",
}