diff --git a/.gitignore b/.gitignore
index 25de861a1db6009fcd0c62652d0e745b4547d292..fc4b5eaa0e905f92701264354a082eb56ebfe503 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,6 +10,9 @@ lib/
bin/
obj/
binarycpython_snippets
+-1.log
+
+
# Standard template
*.org~
*.so
diff --git a/HOW_TO_CONTRIBUTE b/HOW_TO_CONTRIBUTE
index d7597eef602513e45d1c3ed6e7691acf613b3143..d9f7329c6704825e4c7ae154b9a89d628a84a7bb 100644
--- a/HOW_TO_CONTRIBUTE
+++ b/HOW_TO_CONTRIBUTE
@@ -1 +1,6 @@
-You can contribute by contacting me (David) or rob and asking what needs to be done
+You can contribute by contacting me (David) or Rob and asking what needs to be done,
+
+Or resolve one of the many TODOs that are scattered throughout the code.
+
+Probably best to discuss with us first before making pull requests though.
+
diff --git a/binarycpython/utils/custom_logging_functions.py b/binarycpython/utils/custom_logging_functions.py
index e3be2e16ba91c5bfbcaae024c8fe48fdf0aa5156..a17fcbb7adf42efcb1d0e34f233fc5eb5e1c3a40 100644
--- a/binarycpython/utils/custom_logging_functions.py
+++ b/binarycpython/utils/custom_logging_functions.py
@@ -203,7 +203,7 @@ def from_binary_c_config(config_file: str, flag: str) -> str:
# convert and chop off newline
res = res.decode("utf-8").rstrip()
- print(res)
+ # print(res)
return res
@@ -360,7 +360,7 @@ def compile_shared_lib(
)
)
res = subprocess.check_output("{command}".format(command=command), shell=True)
- print(res)
+ # print(res)
if verbose > 0:
if res:
print("Output of compilation command:\n{}".format(res))
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 3d0ea4ca3df13df4cb93a118f2591254ffc81713..65880d59e658c57742f69d8ceeafc77c2462ed2d 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -945,15 +945,11 @@ def cosmic_SFH_madau_dickinson2014(z):
# known value if we try to interpolate outside of the tables.
# There are still some open tasks and improvements that can be made:
#
-# TODO: check all the raise ValueErrors to make them more appropriate
-# TODO: Put the json checking stuff in a different function
-# TODO: Clean up the logging
# TODO: Solve the memory issues that are present.
# Are the interpolators not cleaned?
# TODO: Parallellize the setting up of the interpolators
# TODO: Generalize the code such that we can input other/newer tables
-
########################################################################
import py_rinterpolate
@@ -1128,7 +1124,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
if options["multiplicity_model"] == "Poisson":
multiplicity = Moecache["rinterpolator_multiplicity"].interpolate(
- [np.log10(options["M1"])]
+ [np.log10(options["M_1"])]
)[0]
# Fill the multiplicity array
@@ -1154,7 +1150,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
for n in range(3):
full_fractions_array[n] = (
Moecache["rinterpolator_multiplicity"].interpolate(
- [np.log10(options["M1"])]
+ [np.log10(options["M_1"])]
)[n + 1])
# Set last value
@@ -1264,9 +1260,6 @@ def build_q_table(options, m, p, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
- Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)].destroy()
- gc.collect()
-
#
if not incache:
@@ -1278,7 +1271,7 @@ def build_q_table(options, m, p, verbosity=0):
# 0
# ] # TODO: this lower range must not be lower than Mmin.
# print("build_q_table qmin: {}".format(qmin))
- qmin = options['Mmin']/options['M1']
+ qmin = options['Mmin']/options['M_1']
print("build_q_table qmin: {}".format(qmin))
# qmax = maximum_mass_ratio_for_RLOF(options[m], options[p])
@@ -1823,7 +1816,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
prob_dict['multiplicity'] = multiplicity_probability
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob dict ({})".format(
- options["M1"], prob_dict['multiplicity'], prob_dict
+ options["M_1"], prob_dict['multiplicity'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -1837,11 +1830,11 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# TODO: Create an n-part-powerlaw method that can have breakpoints and slopes. I'm using a three-part powerlaw now.
# TODO: is this actually the correct way? putting the M1 in there? Do we sample in logspace?
- M1_probability = Kroupa2001(options["M1"]) * options["M1"]
- prob_dict['M1'] = M1_probability
+ M1_probability = Kroupa2001(options["M_1"]) * options["M_1"]
+ prob_dict['M_1'] = M1_probability
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob dict ({})".format(
- options["M1"], prob_dict['M1'], prob_dict
+ options["M_1"], prob_dict['M_1'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -1876,7 +1869,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
# Make a table storing Moe's data for q distributions
- if options.get("M2", None) or options.get("M3", None) or options.get("M4", None):
+ if options.get("M_2", None) or options.get("M_3", None) or options.get("M_4", None):
if not Moecache.get("rinterpolator_q", None):
Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
table=Moecache["q_distributions"], # Contains the table of data
@@ -1913,30 +1906,30 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
###############
# Calculation for period of the binary
- if options.get("M2", None):
+ if options.get("M_2", None):
# Separation of the inner binary
options["sep"] = calc_sep_from_period(
- options["M1"], options["M2"], options["P"]
+ options["M_1"], options["M_2"], options["P"]
)
# TODO: add check for min_P with instant RLOF?
# TODO: Actually use the value above.
# Total mass inner binary:
- options["M1+M2"] = options["M1"] + options["M2"]
+ options["M_1+M_2"] = options["M_1"] + options["M_2"]
# Calculate P integral or use cached value
calc_P_integral(
- options, 1, "P_integrals", "rinterpolator_log10P", "M1", verbosity
+ options, 1, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
)
# Set probabilty for P1
p_val = Moecache["rinterpolator_log10P"].interpolate(
- [np.log10(options["M1"]), np.log10(options["P"])]
+ [np.log10(options["M_1"]), np.log10(options["P"])]
)[0]
- p_val = p_val / Moecache["P_integrals"][options["M1"]]
+ p_val = p_val / Moecache["P_integrals"][options["M_1"]]
prob_dict['P'] = p_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended period (m={}, P={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["P"], prob_dict['P'], prob_dict
+ options["M_1"], options["P"], prob_dict['P'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -1951,11 +1944,11 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# we need to construct the q table for the given M1
# subject to qmin = Mmin/M1
- if options.get("M2", None):
+ if options.get("M_2", None):
# Build the table for q
- primary_mass = options["M1"]
- secondary_mass = options["M2"]
- m_label = "M1"
+ primary_mass = options["M_1"]
+ secondary_mass = options["M_2"]
+ m_label = "M_1"
p_label = "P"
# Construct the q table
@@ -1975,9 +1968,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
prob_dict['q'] = q_prob
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
- options["M1"],
+ options["M_1"],
options["P"],
- options["M2"] / options["M1"],
+ options["M_2"] / options["M_1"],
prob_dict['q'],
prob_dict,
),
@@ -1997,19 +1990,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
if not options['ecc']==None:
# Calculate ecc integral or use cached value
calc_e_integral(
- options, "ecc_integrals", "rinterpolator_e", "M1", "P", verbosity
+ options, "ecc_integrals", "rinterpolator_e", "M_1", "P", verbosity
)
- mass_period_string = "{}_{}".format(options["M1"], options["P"])
+ mass_period_string = "{}_{}".format(options["M_1"], options["P"])
# Set probabilty for ecc
ecc_val = Moecache["rinterpolator_e"].interpolate(
- [np.log10(options["M1"]), np.log10(options["P"]), options["ecc"]]
+ [np.log10(options["M_1"]), np.log10(options["P"]), options["ecc"]]
)[0]
ecc_val = ecc_val / Moecache["ecc_integrals"][mass_period_string]
prob_dict['ecc'] = ecc_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended eccentricity (m={}, P={}, ecc={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["P"], options["ecc"], prob_dict['ecc'], prob_dict
+ options["M_1"], options["P"], options["ecc"], prob_dict['ecc'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -2035,7 +2028,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# TODO: Is this a correct assumption?
max_sep = 10.0 * options["sep"] * (1.0 + options["ecc"])
- min_P2 = calc_period_from_sep(options["M1+M2"], options["mmin"], max_sep)
+ min_P2 = calc_period_from_sep(options["M_1+M_2"], options["mmin"], max_sep)
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
@@ -2056,7 +2049,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# hence the separation between the outer star
# and inner binary
options["sep2"] = calc_sep_from_period(
- options["M3"], options["M1+M2"], options["P2"]
+ options["M_3"], options["M_1+M_2"], options["P2"]
)
# Check for cached value of P integral or calculate
@@ -2065,19 +2058,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
min_P2,
"P2_integrals",
"rinterpolator_log10P",
- "M1+M2",
+ "M_1+M_2",
verbosity,
)
# Add the probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
- [np.log10(options["M1+M2"]), np.log10(options["P2"])]
+ [np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
)[0]
- p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
+ p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
prob_dict['P2'] = p_val
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["M2"], options["P2"], prob_dict['P2'], prob_dict
+ options["M_1"], options["M_2"], options["P2"], prob_dict['P2'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -2093,9 +2086,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# subject to qmin = Mmin/(M1+M2)
# Set the variables for the masses and their names
- primary_mass = options["M1+M2"]
- secondary_mass = options["M3"]
- m_label = "M1+M2"
+ primary_mass = options["M_1+M_2"]
+ secondary_mass = options["M_3"]
+ m_label = "M_1+M_2"
p_label = "P2"
# Build q table
@@ -2112,9 +2105,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
].interpolate([secondary_mass / primary_mass])[0]
prob_dict['q2'] = q2_val
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
- options["M1+M2"],
- options["M3"],
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M_1+M_2"],
+ options["M_3"],
options["P"],
secondary_mass / primary_mass,
prob_dict['q2'],
@@ -2142,7 +2135,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# TODO: fix this here
max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
max_per3 = calc_period_from_sep(
- options["M1+M2"], options["mmin"], max_sep3
+ options["M_1+M_2"], options["mmin"], max_sep3
)
# Calculate P integral or use the cached value
@@ -2152,19 +2145,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
min_P2,
"P2_integrals",
"rinterpolator_log10P",
- "M1+M2",
+ "M_1+M_2",
verbosity,
)
# Set probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
- [np.log10(options["M1+M2"]), np.log10(options["P2"])]
+ [np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
)[0]
- p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
+ p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
prob_dict['P3'] = p_val
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
- options["M1"], options["M2"], options["P2"], prob_dict['P3'], prob_dict
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (M_1={} M_2={}, P2={}) probability ({}) to the prob list ({})".format(
+ options["M_1"], options["M_2"], options["P2"], prob_dict['P3'], prob_dict
),
verbosity,
_MS_VERBOSITY_LEVEL,
@@ -2181,9 +2174,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# Make a table storing Moe's data for q distributions
# Build the table for q2
- primary_mass = options["M1+M2"]
- secondary_mass = options["M3"]
- m_label = "M1+M2"
+ primary_mass = options["M_1+M_2"]
+ secondary_mass = options["M_3"]
+ m_label = "M_1+M_2"
p_label = "P2"
# Calculate new q table
@@ -2200,9 +2193,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
].interpolate([secondary_mass / primary_mass])[0]
prob_dict['q3'] = q3_prob
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
- options["M1+M2"],
- options["M3"],
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M_1+M_2"],
+ options["M_3"],
options["P"],
secondary_mass / primary_mass,
prob_dict['q3'],
@@ -2235,8 +2228,8 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# Some info
if multiplicity == 1:
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
- options["M1"],
+ "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
+ options["M_1"],
len(prob_dict),
str(prob_dict),
prob_dict['total_probdens']
@@ -2246,9 +2239,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
elif multiplicity == 2:
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
- options["M1"],
- options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+ "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
+ options["M_1"],
+ options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
options['ecc'] if options.get("ecc", None) else "N/A",
len(prob_dict),
@@ -2260,12 +2253,12 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
elif multiplicity == 3:
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} ({}): {} -> {}".format(
- options["M1"],
- options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+ "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} ({}): {} -> {}".format(
+ options["M_1"],
+ options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
options["ecc"] if options.get("ecc", None) else "N/A",
- options["M3"],
+ options["M_3"],
np.log10(options["P2"]),
options["ecc2"] if options.get("ecc2", None) else "N/A",
len(prob_dict),
@@ -2277,15 +2270,15 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
elif multiplicity == 4:
verbose_print(
- "M&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} M4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
- options["M1"],
- options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+ "M&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} M_4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
+ options["M_1"],
+ options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
np.log10(options["P"]),
options['ecc'] if options.get("ecc", None) else "N/A",
- options["M3"],
+ options["M_3"],
np.log10(options["P2"]),
options["ecc2"] if options.get("ecc2", None) else "N/A",
- options["M4"],
+ options["M_4"],
np.log10(options["P3"]),
options["ecc3"] if options.get("ecc3", None) else "N/A",
len(prob_dict), str(prob_dict), prob_dict['total_probdens']
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 8b958a4b9d7dccd1f5e37f2d1d4acf03b35800bf..4eceb3b62774c2839215f9a9f266b71e4b7718c1 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -20,7 +20,10 @@ import types
from io import StringIO
from typing import Union, Any
-from collections import defaultdict
+from collections import (
+ defaultdict,
+ OrderedDict,
+)
import h5py
import numpy as np
@@ -169,17 +172,17 @@ def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
del new_dict[key]
else:
- print(
- "Error key: {} value: {} type: {} and key: {} value: {} type: {} are not of the same type and cannot be merged".format(
- key,
- dict_1[key],
- type(dict_1[key]),
- key,
- dict_2[key],
- type(dict_2[key]),
- )
+ msg = "Error key: {} value: {} type: {} and key: {} value: {} type: {} are not of the same type and cannot be merged".format(
+ key,
+ dict_1[key],
+ type(dict_1[key]),
+ key,
+ dict_2[key],
+ type(dict_2[key]),
)
- raise ValueError
+
+ print(msg)
+ raise ValueError(msg)
# This is where the keys are the same
else:
@@ -210,26 +213,37 @@ def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
return new_dict
-def get_moe_distefano_dataset(options):
+def get_moe_distefano_dataset(options, verbosity=0):
"""
Function to get the default moe and Distefano dataset or accept a userinput.
"""
if not options.get("file", None):
- print("Using the default Moe and de Stefano 2017 datafile")
+ verbose_print(
+ "Using the default Moe and de Stefano 2017 datafile",
+ verbosity,
+ 1,
+ )
json_data = copy.deepcopy(moe_distefano_data.moe_distefano_2017_data)
else:
if not os.path.isfile(options["file"]):
- print(
+ 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"):
- print("Provided filename does not end with .json")
+ verbose_print(
+ "Provided filename is not a json file",
+ verbosity,
+ 1,
+ )
else:
# Read input data and Clean up the data if there are whitespaces around the keys
@@ -1447,7 +1461,7 @@ def inspect_dict(
type(input_dict[key]) (except if the value is a dict)
"""
- structure_dict = {}
+ structure_dict = OrderedDict() # TODO: check if this still works
#
for key, value in input_dict.items():
@@ -1489,7 +1503,7 @@ def merge_dicts(dict_1: dict, dict_2: dict) -> dict:
"""
# Set up new dict
- new_dict = {}
+ new_dict = OrderedDict() # TODO: check if this still works
#
keys_1 = dict_1.keys()
@@ -1531,6 +1545,11 @@ def merge_dicts(dict_1: dict, dict_2: dict) -> dict:
type(dict_2[key]) in [int, float, np.float64]
):
new_dict[key] = dict_1[key] + dict_2[key]
+ # Exceptions:
+ elif (type(dict_1[key]) in [dict, OrderedDict]) and (
+ type(dict_2[key]) in [dict, OrderedDict]
+ ):
+ new_dict[key] = merge_dicts(dict_1[key], dict_2[key])
else:
print(
@@ -1584,7 +1603,6 @@ def merge_dicts(dict_1: dict, dict_2: dict) -> dict:
#
return new_dict
-
def update_dicts(dict_1: dict, dict_2: dict) -> dict:
"""
Function to update dict_1 with values of dict_2 in a recursive way.
@@ -1607,7 +1625,7 @@ def update_dicts(dict_1: dict, dict_2: dict) -> dict:
"""
# Set up new dict
- new_dict = {}
+ new_dict = OrderedDict() # TODO: check if this still works
#
keys_1 = dict_1.keys()
@@ -1677,6 +1695,141 @@ def update_dicts(dict_1: dict, dict_2: dict) -> dict:
return new_dict
+def count_keys_recursive(input_dict):
+ """
+ Function to count the total amount of keys in a dictionary
+ """
+
+ local_count = 0
+ for key in input_dict.keys():
+ local_count += 1
+ if isinstance(input_dict[key], (dict, OrderedDict)):
+ local_count += count_keys_recursive(input_dict[key])
+ return local_count
+
+
+def recursive_change_key_to_float(input_dict):
+ """
+ Function to recursively change the key to float
+
+ This only works if the dict contains just sub-dicts or numbers/strings.
+ Does not work with lists as values
+ """
+
+ new_dict = OrderedDict() # TODO: check if this still works
+
+ for key in input_dict:
+ if isinstance(input_dict[key], (dict, OrderedDict)):
+ try:
+ num_key = float(key)
+ new_dict[num_key] = recursive_change_key_to_float(copy.deepcopy(input_dict[key]))
+ except ValueError:
+ new_dict[key] = recursive_change_key_to_float(copy.deepcopy(input_dict[key]))
+ else:
+ try:
+ num_key = float(key)
+ new_dict[num_key] = input_dict[key]
+ except ValueError:
+ new_dict[key] = input_dict[key]
+
+ return new_dict
+
+def recursive_change_key_to_string(input_dict):
+ """
+ Function to recursively change the key back to a string but this time in a format that we decide
+ """
+
+ new_dict = OrderedDict() # TODO: check if this still works
+
+ for key in input_dict:
+ if isinstance(input_dict[key], (dict, OrderedDict)):
+ if isinstance(key, (int, float)):
+ string_key = "{:g}".format(key)
+ new_dict[string_key] = recursive_change_key_to_string(copy.deepcopy(input_dict[key]))
+ else:
+ new_dict[key] = recursive_change_key_to_string(copy.deepcopy(input_dict[key]))
+ else:
+ if isinstance(key, (int, float)):
+ string_key = "{:g}".format(key)
+ new_dict[string_key] = input_dict[key]
+ else:
+ new_dict[key] = input_dict[key]
+
+ return new_dict
+
+
+# New method to create a ordered dictionary
+def custom_sort_dict(input_dict):
+ """
+ Returns a dictionary that is ordered, but can handle numbers better than normal OrderedDict
+
+ When the keys of the current dictionary are of mixed type, we first find all the unique types.
+ Sort that list of type names. Then find the values that fit that type.
+ Sort those and append them to the sorted keys list.
+ This is done until all the keys are sorted.
+
+ All objects other than dictionary types are directly return as they are
+ """
+
+ # If the new input is a dictionary, then try to sort it
+ if isinstance(input_dict, (dict, OrderedDict)):
+ new_dict = OrderedDict()
+
+ keys = input_dict.keys()
+
+ # Check if types are the same
+ all_types_keys = []
+ for key in keys:
+ if not type(key) in all_types_keys:
+ all_types_keys.append(type(key))
+
+ # If there are multiple types, then we loop over them and do a piecewise sort
+ if len(all_types_keys) > 1:
+ msg = "Different types in the same dictionary key set"
+
+ # Create a string repr of the type name to sort them afterwards
+ str_types = {repr(el):el for el in all_types_keys}
+
+ # Set up sorted keys list
+ sorted_keys = []
+
+ for key_str_type in sorted(str_types.keys()):
+ cur_type = str_types[key_str_type]
+
+ cur_list = [key for key in keys if isinstance(key, cur_type)]
+ cur_sorted_list = sorted(cur_list)
+
+ sorted_keys = sorted_keys + cur_sorted_list
+ else:
+ sorted_keys = sorted(keys)
+
+ for key in sorted_keys:
+ new_dict[key] = custom_sort_dict(copy.deepcopy(input_dict[key]))
+
+ return new_dict
+ return input_dict
+
+
+def multiply_values_dict(input_dict, factor):
+ """
+ Function that goes over dictionary recursively and multiplies the value if possible by a factor
+
+ If the key equals "general_info", the multiplication gets skipped
+ """
+
+ for key in input_dict:
+ if not key == 'general_info':
+ if isinstance(input_dict[key], (dict, OrderedDict)):
+ input_dict[key] = multiply_values_dict(input_dict[key], factor)
+ else:
+ if isinstance(input_dict[key], (int, float)):
+ input_dict[key] = input_dict[key] * factor
+
+ return input_dict
+
+
+#####
+
def extract_ensemble_json_from_string(binary_c_output: str) -> dict:
"""
Function to extract the ensemble_json information from a raw binary_c output string
@@ -1788,7 +1941,7 @@ def binaryc_json_serializer(obj: Any) -> Any:
Function objects will be turned into str representations of themselves
Args:
- obj: obj being process
+ obj: The object that might not be serializable
Returns:
Either string representation of object if the object is a function, or the object itself
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 5e5c9498304cf6172c05a74aa2fe1c078596df46..27d32cc9fcba345b64548731a2ea69247c0ae56e 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -37,10 +37,12 @@ import argparse
import importlib.util
import multiprocessing
from typing import Union, Any
+from collections import (
+ OrderedDict,
+)
import setproctitle
import py_rinterpolate
-
from binarycpython.utils.grid_options_defaults import (
grid_options_defaults_dict,
moe_distefano_default_options,
@@ -64,8 +66,13 @@ from binarycpython.utils.functions import (
update_dicts,
extract_ensemble_json_from_string,
get_moe_distefano_dataset,
+ recursive_change_key_to_float,
+ custom_sort_dict,
+ recursive_change_key_to_string,
+ multiply_values_dict,
)
+
# from binarycpython.utils.hpc_functions import (
# get_condor_version,
# get_slurm_version,
@@ -79,6 +86,9 @@ from binarycpython.utils.distribution_functions import (
LOG_LN_CONVERTER,
fill_data,
get_max_multiplicity,
+ Arenou2010_binary_fraction,
+ raghavan2010_binary_fraction,
+ Moe_de_Stefano_2017_multiplicity_fractions,
)
from binarycpython import _binary_c_bindings
@@ -115,6 +125,14 @@ class Population:
# Custom options
self.custom_options = {}
+ # Load M&s options
+ self.grid_options['m&s_options'] = copy.deepcopy(moe_distefano_default_options)
+
+ # Write M&S 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") as f:
+ f.write(json.dumps(self.grid_options['m&s_options'], indent=4))
+
# Argline dict
self.argline_dict = {}
@@ -930,6 +948,33 @@ class Population:
for _ in range(amt_cores):
job_queue.put("STOP")
+ def format_ensemble_results(self, ensemble_dictionary):
+ """
+ Function to handle all the steps of formatting the ensemble output again.
+
+ Input:
+ ensemble_dictionary: dictionary containing all the ensemble results
+ """
+
+ original_ensemble_results = ensemble_dictionary
+
+ float_format_ensemble_results = recursive_change_key_to_float(original_ensemble_results)
+ del original_ensemble_results
+ gc.collect()
+
+ # Then sort the dictionary
+ sorted_ensemble_results = custom_sort_dict(float_format_ensemble_results)
+ del float_format_ensemble_results
+ gc.collect()
+
+ # Then Change the keys back to a string but with a %g format.
+ reformatted_ensemble_results = recursive_change_key_to_string(sorted_ensemble_results)
+ del sorted_ensemble_results
+ gc.collect()
+
+ # Put back in the dictionary
+ return reformatted_ensemble_results
+
def _evolve_population_grid(self):
"""
Function that handles running the population using multiprocessing.
@@ -982,7 +1027,7 @@ class Population:
# Handle the results by merging all the dictionaries. How that merging happens exactly is
# described in the merge_dicts description.
- combined_output_dict = {}
+ combined_output_dict = OrderedDict()
sentinel = object()
for output_dict in iter(result_queue.get, sentinel):
@@ -990,17 +1035,46 @@ class Population:
if result_queue.empty():
break
+ # Extra ensemble result manipulation:
+ combined_output_dict['ensemble_results']['ensemble'] = self.format_ensemble_results(combined_output_dict['ensemble_results'].get('ensemble', {}))
+ gc.collect()
+
+ # Take into account that we run this on multiple cores
+ combined_output_dict["_total_probability_weighted_mass_run"] = combined_output_dict["_total_probability_weighted_mass_run"]
+
# Put the values back as object properties
self.grid_results = combined_output_dict["results"]
+
+ #################################
+ # Put Ensemble results
self.grid_ensemble_results = combined_output_dict[
"ensemble_results"
] # Ensemble results are also passed as output from that dictionary
- # Add some metadata
- self.grid_ensemble_results["population_id"] = self.grid_options[
+ # Add metadata
+ self.grid_ensemble_results["metadata"] = {}
+ self.grid_ensemble_results["metadata"]["population_id"] = self.grid_options[
"_population_id"
]
+ self.grid_ensemble_results["metadata"][
+ "total_probability_weighted_mass"
+ ] = combined_output_dict["_total_probability_weighted_mass_run"]
+ self.grid_ensemble_results["metadata"]["factored_in_probability_weighted_mass"] = False
+ if self.grid_options['ensemble_factor_in_probability_weighted_mass']:
+ multiply_values_dict(self.grid_ensemble_results['ensemble'], 1/self.grid_ensemble_results["metadata"]['total_probability_weighted_mass'])
+ self.grid_ensemble_results["metadata"]["factored_in_probability_weighted_mass"] = True
+
+ # Add settings of the populations
+ all_info = self.return_all_info(
+ include_population_settings=True,
+ include_binary_c_defaults=True,
+ include_binary_c_version_info=True,
+ include_binary_c_help_all=True,
+ )
+ self.grid_ensemble_results["metadata"]["settings"] = json.loads(json.dumps(all_info, default=binaryc_json_serializer))
+ ##############################
+ # Update grid options
self.grid_options["_failed_count"] = combined_output_dict["_failed_count"]
self.grid_options["_failed_prob"] = combined_output_dict["_failed_prob"]
self.grid_options["_failed_systems_error_codes"] = list(
@@ -1319,6 +1393,10 @@ class Population:
ensemble_raw_output
) # Load this into a dict so that we can combine it later
+ # Extra ensemble result manipulation:
+ # We need to reformat and multiply by a factor
+ ensemble_json["ensemble"] = self.format_ensemble_results(ensemble_json["ensemble"])
+ gc.collect()
else:
# If we do not allow this, automatically we will export this to the data_dir, in
# some formatted way
@@ -1421,14 +1499,9 @@ class Population:
stream_logger.debug(f"Process-{self.process_ID} is finished.")
# Clean up the interpolators if they exist
+
# TODO: make a cleanup function for the individual threads
# TODO: make sure this is necessary. Actually its probably not, because we have a centralized queue
- verbose_print(
- "Process {}: Cleaning up interpolators".format(ID),
- self.grid_options["verbosity"],
- 1,
- )
- self._clean_interpolators()
return
@@ -1665,9 +1738,6 @@ class Population:
# Unload/free custom_logging_code
# TODO: cleanup custom logging code.
- # Clean up the interpolator functions
- self._clean_interpolators()
-
# Also remove the rest of the contents
keys_moecache = list(Moecache.keys())
for key in keys_moecache:
@@ -1714,6 +1784,7 @@ class Population:
# Import packages
code_string += "import math\n"
code_string += "import numpy as np\n"
+ code_string += "from collections import OrderedDict\n"
code_string += "from binarycpython.utils.distribution_functions import *\n"
code_string += "from binarycpython.utils.spacing_functions import *\n"
code_string += "from binarycpython.utils.useful_funcs import *\n"
@@ -2031,7 +2102,7 @@ class Population:
code_string += indent * (depth + 1) + "#" * 40 + "\n"
code_string += (
indent * (depth + 1)
- + "# Code below is for finalising the handling of this iteration of the parameter\n"
+ + "# Code below is for finalising the handling of this iteration of the parameter {}\n".format(grid_variable["name"])
)
# Set phasevol
@@ -2146,7 +2217,15 @@ class Population:
+ "# Code below will get evaluated for every generated system\n"
)
- # Calculate value
+ # Factor in the custom weight input
+ code_string += (
+ indent * (depth + 1)
+ + "\n"
+ )
+ code_string += (
+ indent * (depth + 1)
+ + "# Weigh the probability by a custom weighting factor\n"
+ )
code_string += (
indent * (depth + 1)
+ 'probability = self.grid_options["weight"] * probabilities_list[{}]'.format(
@@ -2154,22 +2233,54 @@ class Population:
)
+ "\n"
)
- # TODO: ask rob if just replacing this with probability is enough
+
+ # Take into account the multiplicity fraction:
+ code_string += (
+ indent * (depth + 1)
+ + "\n"
+ )
+ code_string += (
+ indent * (depth + 1)
+ + "# Factor the multiplicity fraction into the probability\n"
+ )
+ code_string += (
+ indent * (depth + 1)
+ + 'probability = probability * self._calculate_multiplicity_fraction(parameter_dict)'
+ + "\n"
+ )
+
+ # Add division by amount of repeats
+ code_string += (
+ indent * (depth + 1)
+ + "\n"
+ )
+ code_string += (
+ indent * (depth + 1)
+ + "# Divide the probability by the amount of repeats\n"
+ )
+
code_string += (
indent * (depth + 1)
- # + 'repeat_probability = probability / self.grid_options["repeat"]'
+ 'probability = probability / self.grid_options["repeat"]'
+ "\n"
)
- # For each repeat of the system this has to be done yes.
+ # Now we yield the system self.grid_options["repeat"] times.
+ code_string += (
+ indent * (depth + 1)
+ + "\n"
+ )
+ code_string += (
+ indent * (depth + 1)
+ + "# Loop over the repeats\n"
+ )
code_string += (
indent * (depth + 1) + 'for _ in range(self.grid_options["repeat"]):' + "\n"
)
code_string += indent * (depth + 2) + "_total_starcount += 1\n"
- # set probability and phasevol values
+ # set probability and phasevol values into the system dict
code_string += (
indent * (depth + 2)
+ 'parameter_dict["{}"] = {}'.format("probability", "probability")
@@ -2196,16 +2307,11 @@ class Population:
# Handling of what is returned, or what is not.
# TODO: think of whether this is a good method
# code_string += indent * (depth + 2) + "if (self.grid_options['multiplicity'] >= 2): print('phasevol_q: ',phasevol_q); print('phasevol_log10per: ',phasevol_log10per);\n"
- # code_string += indent * (depth + 2) + "print('phasevol_lnm1: ',phasevol_lnm1); print('phasevol_multiplicity: ',phasevol_multiplicity);\n"
- # code_string += indent * (depth + 2) + "print(probabilities_list)\n"
- # code_string += indent * (depth + 2) + "print(parameter_dict)\n"
- # code_string += indent * (depth + 2) + "print('YOO IK GA LEKKER NOG EEN RONDJE')\n"
code_string += indent * (depth + 2) + "yield(parameter_dict)\n"
# If its a dry run, dont do anything with it
else:
# code_string += indent * (depth + 2) + "if (self.grid_options['multiplicity'] >= 2): print(phasevol_q)\n"
- # code_string += indent * (depth + 2) + "print(parameter_dict)\n"
code_string += indent * (depth + 2) + "pass\n"
code_string += indent * (depth + 1) + "#" * 40 + "\n"
@@ -3066,11 +3172,6 @@ class Population:
self.logger.addHandler(handler_file)
self.logger.addHandler(handler_stdout)
- # def join_result_dicts(self):
- # """
- # Function to join the result dictionaries
- # """
-
def _check_binary_c_error(self, binary_c_output, system_dict):
"""
Function to check whether binary_c throws an error and handle accordingly.
@@ -3145,71 +3246,45 @@ class Population:
f.write(argstring + "\n")
else:
verbose_print(
- "binary_c had 0 output. Which is strange",
+ "binary_c had 0 output. Which is strange. Make sure you mean to do this. If there is actually ensemble output being generated then its fine.",
self.grid_options["verbosity"],
2,
)
- ################################################################################################
- def Moe_de_Stefano_2017(self, options=None):
+
+ def set_moe_distefano_settings(self, options=None):
"""
+ Function to set user input configurations for the moe & distefano methods
+
+ If nothing is passed then we just use the default options
"""
- # Require input options
- if not options:
- msg = "Please provide a options dictionary."
- verbose_print(
- "\tMoe_de_Stefano_2017: {}".format(msg),
- self.grid_options["verbosity"],
- 0,
- )
- raise ValueError(msg)
-
- # Create the tmp dir
- ms_tmp_dir = os.path.join(self.grid_options["tmp_dir"], "moe_distefano")
- os.makedirs(ms_tmp_dir, exist_ok=True)
-
- # TODO: put this in defaults
- default_options = copy.deepcopy(moe_distefano_default_options)
-
- # Take the option dictionary that was given and override.
- options = update_dicts(default_options, options)
+ if options:
+ # Take the option dictionary that was given and override.
+ options = update_dicts(self.grid_options['m&s_options'], options)
+ self.grid_options['m&s_options'] = copy.deepcopy(dict(options))
- # Write options to a file
- with open(os.path.join(ms_tmp_dir, "moeopts.dat"), "w") as f:
- f.write(json.dumps(options, indent=4))
-
- json_data = get_moe_distefano_dataset(options)
-
- # entry of log10M1 is a list containing 1 dict. We can take the dict out of the list
-
- json_data["log10M1"] = json_data["log10M1"][0]
-
- # Get all the masses
- logmasses = sorted(json_data["log10M1"].keys())
- if not logmasses:
- msg = "The table does not contain masses."
- verbose_print(
- "\tMoe_de_Stefano_2017: {}".format(msg),
- self.grid_options["verbosity"],
- 0,
- )
- raise ValueError(msg)
+ # 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") as f:
+ f.write(json.dumps(self.grid_options['m&s_options'], indent=4))
- # Write to file
- with open(os.path.join(ms_tmp_dir, "moe.log"), "w") as logfile:
- logfile.write("logâ‚â‚€Masses(M☉) {}\n".format(logmasses))
+ def _load_moe_distefano_data(self):
+ """
+ Function to load the moe & distefano data
+ """
- # 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"])
- dlog10P = float(logperiods[1]) - float(logperiods[0])
+ # Only if the grid is loaded and moecache contains information
+ if (not self.grid_options['_loaded_ms_data']) and not Moecache:
+ # Load the json data
+ json_data = get_moe_distefano_dataset(self.grid_options['m&s_options'], verbosity=self.grid_options['verbosity'])
- current_logperiods = sorted(json_data["log10M1"][logmass]["logP"])
+ # entry of log10M1 is a list containing 1 dict. We can take the dict out of the list
+ json_data["log10M1"] = json_data["log10M1"][0]
- if not (logperiods == current_logperiods):
- msg = "Period values are not consistent throughout the dataset"
+ # Get all the masses
+ logmasses = sorted(json_data["log10M1"].keys())
+ if not logmasses:
+ msg = "The table does not contain masses."
verbose_print(
"\tMoe_de_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
@@ -3217,15 +3292,22 @@ class Population:
)
raise ValueError(msg)
- ############################################################
- # log10period binwidth : of course this assumes a fixed
- # binwidth, so we check for this too.
+ # 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") as logfile:
+ logfile.write("logâ‚â‚€Masses(M☉) {}\n".format(logmasses))
- 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"
+ # 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"])
+ dlog10P = float(logperiods[1]) - float(logperiods[0])
+
+ current_logperiods = sorted(json_data["log10M1"][logmass]["logP"])
+
+ if not (logperiods == current_logperiods):
+ msg = "Period values are not consistent throughout the dataset"
verbose_print(
"\tMoe_de_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
@@ -3233,69 +3315,100 @@ class Population:
)
raise ValueError(msg)
- # Write to file
- with open(
- os.path.join(self.grid_options["tmp_dir"], "moe_distefano", "moe.log"), "a"
- ) 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"],
- ]
- )
+ ############################################################
+ # log10period binwidth : of course this assumes a fixed
+ # binwidth, so we check for this too.
- ############################################################
- # a small log10period which we can shift just outside the
- # table to force integration out there to zero
- epslog10P = 1e-8 * dlog10P
+ 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_de_Stefano_2017: {}".format(msg),
+ self.grid_options["verbosity"],
+ 0,
+ )
+ raise ValueError(msg)
- ############################################################
- # loop over either binary or triple-outer periods
- first = 1
+ # 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"
+ ) as logfile:
+ logfile.write("logâ‚â‚€Periods(days) {}\n".format(logperiods))
- # Go over the periods
- for logperiod in logperiods:
- # print("logperiod: {}".format(logperiod))
- ############################################################
- # 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
+ # Fill the global dict
+ for logmass in logmasses:
+ # Create the multiplicity table
+ if not Moecache.get("multiplicity_table", None):
+ Moecache["multiplicity_table"] = []
- if first:
- first = 0
+ # 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"],
+ ]
+ )
- # Create the multiplicity table
- if not Moecache.get("period_distributions", None):
- Moecache["period_distributions"] = []
+ ############################################################
+ # 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:
+ # print("logperiod: {}".format(logperiod))
############################################################
- # lower bound the period distributions to zero probability
- Moecache["period_distributions"].append(
- [
- float(logmass),
- float(logperiod) - 0.5 * dlog10P - epslog10P,
- 0.0,
- 0.0,
- ]
- )
+ # 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) - 0.5 * dlog10P,
+ float(logperiod),
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_bin_frac_p_dist"
]
@@ -3306,132 +3419,118 @@ class Population:
/ 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 'normalize'
+ qdata = 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 'normalize'
+ ecc_data = 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(logperiod),
- json_data["log10M1"][logmass]["logP"][logperiod][
+ float(logperiods[-1]) + 0.5 * dlog10P, # TODO: why this shift?
+ json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_bin_frac_p_dist"
]
/ dlog10P,
- json_data["log10M1"][logmass]["logP"][logperiod][
+ json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"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 'normalize'
- qdata = 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 'normalize'
- ecc_data = fill_data(
- eccs, json_data["log10M1"][logmass]["logP"][logperiod]["e"]
+ Moecache["period_distributions"].append(
+ [
+ float(logmass),
+ float(logperiods[-1]) + 0.5 * dlog10P + epslog10P,
+ 0.0,
+ 0.0,
+ ]
)
- # 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?
- 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,
- ]
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Length period_distributions table: {}".format(
+ len(Moecache["period_distributions"])
+ ),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
)
- Moecache["period_distributions"].append(
- [
- float(logmass),
- float(logperiods[-1]) + 0.5 * dlog10P + epslog10P,
- 0.0,
- 0.0,
- ]
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Length multiplicity table: {}".format(
+ len(Moecache["multiplicity_table"])
+ ),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
+ )
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Length q table: {}".format(
+ len(Moecache["q_distributions"])
+ ),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
+ )
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Length ecc table: {}".format(
+ len(Moecache["ecc_distributions"])
+ ),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
)
- verbose_print(
- "\tMoe_de_Stefano_2017: Length period_distributions table: {}".format(
- len(Moecache["period_distributions"])
- ),
- self.grid_options["verbosity"],
- _MS_VERBOSITY_LEVEL,
- )
- verbose_print(
- "\tMoe_de_Stefano_2017: Length multiplicity table: {}".format(
- len(Moecache["multiplicity_table"])
- ),
- self.grid_options["verbosity"],
- _MS_VERBOSITY_LEVEL,
- )
- verbose_print(
- "\tMoe_de_Stefano_2017: Length q table: {}".format(
- len(Moecache["q_distributions"])
- ),
- self.grid_options["verbosity"],
- _MS_VERBOSITY_LEVEL,
- )
- verbose_print(
- "\tMoe_de_Stefano_2017: Length ecc table: {}".format(
- len(Moecache["ecc_distributions"])
- ),
- self.grid_options["verbosity"],
- _MS_VERBOSITY_LEVEL,
- )
-
- # Write to logfile
- with open(os.path.join(ms_tmp_dir, "moecache.json"), "w") as cache_filehandle:
- cache_filehandle.write(json.dumps(Moecache, indent=4))
-
+ # Write to logfile
+ 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") as cache_filehandle:
+ cache_filehandle.write(json.dumps(Moecache, indent=4))
- ############################################################
- # construct the grid here
- #
+ # Signal that the data has been loaded
+ self.grid_options['_loaded_ms_data'] = True
- ############################################################
- # make a version of the options which can be passed to
- # the probability density function in the gridcode
+ def _set_moe_distefano_distributions(self):
+ """
+ Function to set the moe & distefano distribution
+ """
############################################################
# first, the multiplicity, this is 1,2,3,4, ...
# for singles, binaries, triples, quadruples, ...
- max_multiplicity = get_max_multiplicity(options["multiplicity_modulator"])
+ max_multiplicity = get_max_multiplicity(self.grid_options['m&s_options']["multiplicity_modulator"])
verbose_print(
"\tMoe_de_Stefano_2017: Max multiplicity = {}".format(max_multiplicity),
self.grid_options["verbosity"],
@@ -3448,11 +3547,11 @@ class Population:
valuerange=[1, max_multiplicity],
resolution=1,
spacingfunc="range(1, 5)",
- precode='self.grid_options["multiplicity"] = multiplicity; options={}'.format(
- options
+ precode='self.grid_options["multiplicity"] = multiplicity; self.bse_options["multiplicity"] = multiplicity; options={}'.format(
+ self.grid_options['m&s_options']
),
condition="({}[multiplicity-1] > 0)".format(
- str(options["multiplicity_modulator"])
+ str(self.grid_options['m&s_options']["multiplicity_modulator"])
),
gridtype="edge",
dphasevol=-1,
@@ -3472,19 +3571,19 @@ class Population:
longname="Primary mass",
resolution="options['resolutions']['M'][0]",
spacingfunc="const(np.log({}), np.log({}), {})".format(
- options["ranges"]["M"][0],
- options["ranges"]["M"][1],
- options["resolutions"]["M"][0],
+ self.grid_options['m&s_options']["ranges"]["M"][0],
+ self.grid_options['m&s_options']["ranges"]["M"][1],
+ self.grid_options['m&s_options']["resolutions"]["M"][0],
),
valuerange=[
- "np.log({})".format(options["ranges"]["M"][0]),
- "np.log({})".format(options["ranges"]["M"][1]),
+ "np.log({})".format(self.grid_options['m&s_options']["ranges"]["M"][0]),
+ "np.log({})".format(self.grid_options['m&s_options']["ranges"]["M"][1]),
],
gridtype="centred",
dphasevol="dlnm1",
- precode='M_1 = np.exp(lnm1); options["M1"]=M_1',
+ precode='M_1 = np.exp(lnm1); options["M_1"]=M_1',
probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
- str(options)[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
+ str(dict(self.grid_options['m&s_options']))[1:-1], "'multiplicity': multiplicity", "'M_1': M_1"
),
)
@@ -3495,23 +3594,23 @@ class Population:
name="log10per",
parameter_name="orbital_period",
longname="log10(Orbital_Period)",
- resolution=options["resolutions"]["logP"][0],
+ resolution=self.grid_options['m&s_options']["resolutions"]["logP"][0],
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=[options["ranges"]["logP"][0], options["ranges"]["logP"][1]],
+ valuerange=[self.grid_options['m&s_options']["ranges"]["logP"][0], self.grid_options['m&s_options']["ranges"]["logP"][1]],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
- options["resolutions"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][1],
+ self.grid_options['m&s_options']["resolutions"]["logP"][0],
),
precode="""orbital_period = 10.0**log10per
qmin={}/M_1
qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
""".format(
- options.get("Mmin", 0.07)
+ self.grid_options['m&s_options'].get("Mmin", 0.07)
),
) # TODO: change the maximum_mass_ratio_for_RLOF
@@ -3521,14 +3620,14 @@ qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
parameter_name="M_2",
longname="Mass ratio",
valuerange=[
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "options.get('Mmin', 0.07)/M_1",
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "qmax",
],
- resolution=options["resolutions"]["M"][1],
+ resolution=self.grid_options['m&s_options']["resolutions"]["M"][1],
probdist=1,
gridtype="centred",
dphasevol="dq",
@@ -3537,35 +3636,35 @@ M_2 = q * M_1
sep = calc_sep_from_period(M_1, M_2, orbital_period)
""",
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", [None, None])[0]
- else "{}/M_1".format(options.get("Mmin", 0.07)),
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", [None, None])[1]
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", [None, None])[0]
+ else "{}/M_1".format(self.grid_options['m&s_options'].get("Mmin", 0.07)),
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", [None, None])[1]
else "qmax",
- options["resolutions"]["M"][1],
+ self.grid_options['m&s_options']["resolutions"]["M"][1],
),
)
# (optional) binaries: eccentricity
- if options["resolutions"]["ecc"][0] > 0:
+ if self.grid_options['m&s_options']["resolutions"]["ecc"][0] > 0:
self.add_grid_variable(
name="ecc",
parameter_name="eccentricity",
longname="Eccentricity",
- resolution=options["resolutions"]["ecc"][0],
+ resolution=self.grid_options['m&s_options']["resolutions"]["ecc"][0],
probdist=1,
gridtype="centred",
dphasevol="decc",
precode="eccentricity=ecc",
valuerange=[
- options["ranges"]["ecc"][0], # Just fail if not defined.
- options["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["ranges"]["ecc"][0], # Just fail if not defined.
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["ecc"][0], # Just fail if not defined.
- options["ranges"]["ecc"][1],
- options["resolutions"]["ecc"][0],
+ self.grid_options['m&s_options']["ranges"]["ecc"][0], # Just fail if not defined.
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["resolutions"]["ecc"][0],
),
)
@@ -3576,45 +3675,44 @@ sep = calc_sep_from_period(M_1, M_2, orbital_period)
name="log10per2",
parameter_name="orbital_period_triple",
longname="log10(Orbital_Period2)",
- resolution=options["resolutions"]["logP"][1],
+ resolution=self.grid_options['m&s_options']["resolutions"]["logP"][1],
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=[
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
+ self.grid_options['m&s_options']["ranges"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][1],
],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
- options["resolutions"]["logP"][1],
+ self.grid_options['m&s_options']["ranges"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][1],
+ self.grid_options['m&s_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(
- options.get("Mmin", 0.07)
+ self.grid_options['m&s_options'].get("Mmin", 0.07)
),
)
# Triples: mass ratio
# Note, the mass ratio is M_outer/M_inner
- print("M2 resolution {}".format(options["resolutions"]["M"][2]))
self.add_grid_variable(
name="q2",
parameter_name="M_3",
longname="Mass ratio outer/inner",
valuerange=[
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "options.get('Mmin', 0.07)/(M_1+M_2)",
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "q2max",
],
- resolution=options["resolutions"]["M"][2],
+ resolution=self.grid_options['m&s_options']["resolutions"]["M"][2],
probdist=1,
gridtype="centred",
dphasevol="dq2",
@@ -3624,35 +3722,35 @@ sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
eccentricity2=0
""",
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "options.get('Mmin', 0.07)/(M_1+M_2)",
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "q2max",
- options["resolutions"]["M"][2],
+ self.grid_options['m&s_options']["resolutions"]["M"][2],
),
)
# (optional) triples: eccentricity
- if options["resolutions"]["ecc"][1] > 0:
+ if self.grid_options['m&s_options']["resolutions"]["ecc"][1] > 0:
self.add_grid_variable(
name="ecc2",
parameter_name="eccentricity2",
longname="Eccentricity of the triple",
- resolution=options["resolutions"]["ecc"][1],
+ resolution=self.grid_options['m&s_options']["resolutions"]["ecc"][1],
probdist=1,
gridtype="centred",
dphasevol="decc2",
precode="eccentricity2=ecc2",
valuerange=[
- options["ranges"]["ecc"][0], # Just fail if not defined.
- options["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["ranges"]["ecc"][0], # Just fail if not defined.
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["ecc"][0], # Just fail if not defined.
- options["ranges"]["ecc"][1],
- options["resolutions"]["ecc"][1],
+ self.grid_options['m&s_options']["ranges"]["ecc"][0], # Just fail if not defined.
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["resolutions"]["ecc"][1],
),
)
@@ -3662,44 +3760,43 @@ eccentricity2=0
name="log10per3",
parameter_name="orbital_period_quadruple",
longname="log10(Orbital_Period3)",
- resolution=options["resolutions"]["logP"][2],
+ resolution=self.grid_options['m&s_options']["resolutions"]["logP"][2],
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=[
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
+ self.grid_options['m&s_options']["ranges"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][1],
],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["logP"][0],
- options["ranges"]["logP"][1],
- options["resolutions"]["logP"][2],
+ self.grid_options['m&s_options']["ranges"]["logP"][0],
+ self.grid_options['m&s_options']["ranges"]["logP"][1],
+ self.grid_options['m&s_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(
- options.get("Mmin", 0.07)
+ self.grid_options['m&s_options'].get("Mmin", 0.07)
),
)
# Quadruple: mass ratio : M_outer / M_inner
- print("M3 resolution {}".format(options["resolutions"]["M"][3]))
self.add_grid_variable(
name="q3",
parameter_name="M_4",
longname="Mass ratio outer low/outer high",
valuerange=[
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", None)
- else "options.get('Mmin', 0.07)/(m3)",
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
+ else "options.get('Mmin', 0.07)/(M_3)",
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "q3max",
],
- resolution=options["resolutions"]["M"][3],
+ resolution=self.grid_options['m&s_options']["resolutions"]["M"][3],
probdist=1,
gridtype="centred",
dphasevol="dq3",
@@ -3707,41 +3804,41 @@ q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
M_4 = q3 * M_3
sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
eccentricity3=0
- """,
+""",
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["q"][0]
- if options.get("ranges", {}).get("q", None)
- else "options.get('Mmin', 0.07)/(m3)",
- options["ranges"]["q"][1]
- if options.get("ranges", {}).get("q", None)
+ self.grid_options['m&s_options']["ranges"]["q"][0]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
+ else "options.get('Mmin', 0.07)/(M_3)",
+ self.grid_options['m&s_options']["ranges"]["q"][1]
+ if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
else "q3max",
- options["resolutions"]["M"][2],
+ self.grid_options['m&s_options']["resolutions"]["M"][2],
),
)
# (optional) triples: eccentricity
- if options["resolutions"]["ecc"][2] > 0:
+ if self.grid_options['m&s_options']["resolutions"]["ecc"][2] > 0:
self.add_grid_variable(
name="ecc3",
parameter_name="eccentricity3",
longname="Eccentricity of the triple+quadruple/outer binary",
- resolution=options["resolutions"]["ecc"][2],
+ resolution=self.grid_options['m&s_options']["resolutions"]["ecc"][2],
probdist=1,
gridtype="centred",
dphasevol="decc3",
precode="eccentricity3=ecc3",
valuerange=[
- options["ranges"]["ecc"][
+ self.grid_options['m&s_options']["ranges"]["ecc"][
0
], # Just fail if not defined.
- options["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
],
spacingfunc="const({}, {}, {})".format(
- options["ranges"]["ecc"][
+ self.grid_options['m&s_options']["ranges"]["ecc"][
0
], # Just fail if not defined.
- options["ranges"]["ecc"][1],
- options["resolutions"]["ecc"][2],
+ self.grid_options['m&s_options']["ranges"]["ecc"][1],
+ self.grid_options['m&s_options']["resolutions"]["ecc"][2],
),
)
@@ -3751,13 +3848,13 @@ eccentricity3=0
# the real probability. The trick we use is to strip the options_dict as a string
# and add some keys to it:
- probdist_addition = "Moe_de_Stefano_2017_pdf({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens']".format(
- str(options)[1:-1],
+ updated_options = "{{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}".format(
+ str(dict(self.grid_options['m&s_options']))[1:-1],
'"multiplicity": multiplicity',
- '"M1": M_1',
- '"M2": M_2',
- '"M3": M_3',
- '"M4": M_4',
+ '"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',
@@ -3766,11 +3863,20 @@ eccentricity3=0
'"ecc3": eccentricity3',
)
+ probdist_addition = "Moe_de_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
+ # Add the options to the grid
+ # self.grid_options["_grid_variables"][self.last_grid_variable()][
+ # "precode"
+ # ] += "; self.grid_options['m&s_ensemble_options'] = {}".format(o)
+
verbose_print(
"\tMoe_de_Stefano_2017: Added final call to the pdf function".format(
max_multiplicity
@@ -3782,6 +3888,25 @@ eccentricity3=0
# Signal that the M&S grid has been set
self.grid_options['_set_ms_grid'] = True
+ ################################################################################################
+ def Moe_de_Stefano_2017(self, options=None):
+ """
+ Function to handle setting the user input settings,
+ set up the data and load that into interpolators and
+ and then set the distribution functions
+
+ Takes a dictionary as its only argument
+ """
+
+ # Set the userinput
+ self.set_moe_distefano_settings(options=options)
+
+ # Load the data
+ self._load_moe_distefano_data()
+
+ # construct the grid here
+ self._set_moe_distefano_distributions()
+
def _clean_interpolators(self):
"""
Function to clean up the interpolators after a run
@@ -3799,3 +3924,93 @@ eccentricity3=0
Moecache[key].destroy()
del Moecache[key]
gc.collect()
+
+ ##### Unsorted functions
+ 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
+ """
+
+ # Just return 1 if no option has been chosen
+ if self.grid_options['multiplicity_fraction_function'] == 0:
+ print("Chosen not to use any multiplicity fraction.")
+ 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'] == 1:
+ # Arenou 2010 will be used
+ verbose_print(
+ "_calculate_multiplicity_fraction: Using Arenou 2010 to calculate multiplicity fractions",
+ self.grid_options["verbosity"],
+ 1,
+ )
+
+ binary_fraction = 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'] == 2:
+ # Raghavan 2010 will be used
+ verbose_print(
+ "_calculate_multiplicity_fraction: Using Rhagavan 2010 to calculate multiplicity fractions",
+ self.grid_options["verbosity"],
+ 1,
+ )
+
+ binary_fraction = 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'] == 3:
+ # We need to check several things now here:
+
+ # First, are the options for the M&S grid set? On start it is filled with the default settings
+ if not self.grid_options['m&s_options']:
+ msg = "The M&S options do not seem to be set properly. The value is {}".format(self.grid_options['m&s_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"],
+ 1,
+ )
+
+ # Load the data
+ self._load_moe_distefano_data()
+
+ # record the prev value
+ prev_M1_value_ms = self.grid_options['m&s_options'].get('M_1', None)
+
+ # Set value of M1 of the current system
+ self.grid_options['m&s_options']['M_1'] = system_dict['M_1']
+
+ # Calculate the multiplicity fraction
+ multiplicity_fraction_list = Moe_de_Stefano_2017_multiplicity_fractions(self.grid_options['m&s_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['m&s_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)
+
+ verbose_print(
+ "Multiplicity: {} multiplicity_fraction: {}".format(system_dict['multiplicity'], multiplicity_fraction_dict[system_dict['multiplicity']]),
+ self.grid_options["verbosity"],
+ 1,
+ )
+
+ return multiplicity_fraction_dict[system_dict['multiplicity']]
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index 667b552b89e36e7f8544a2e89ead0bf9d7331ee1..5d9240ecb6823bd2e5b7031faf7c9fefd1b4f143 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -11,8 +11,6 @@ There are several other functions in this module, mostly to generate help texts
With this its also possible to automatically generate a document containing all the setting names + descriptions.
All the options starting with _ should not be changed by the user except when you really know what you're doing (which is probably hacking the code :P)
-
-# TODO: add moe distefano
"""
import os
@@ -24,26 +22,24 @@ _MS_VERBOSITY_LEVEL = 5
_MS_VERBOSITY_INTERPOLATOR_LEVEL = 6
_MS_VERBOSITY_INTERPOLATOR_EXTRA_LEVEL = 7
-
# Options dict
grid_options_defaults_dict = {
##########################
# general (or unordered..)
##########################
"amt_cores": 1, # total amount of cores used to evolve the population
- "binary": 0, # FLag on whether the systems are binary systems or single systems.
"parse_function": None, # FUnction to parse the output with.
+ "multiplicity_fraction_function": 0, # Which multiplicity fraction function to use. 0: None, 1: Arenou 2010, 2: Rhagavan 2010, 3: M&S 2017
"tmp_dir": temp_dir(), # Setting the temp dir of the program
"_main_pid": -1, # Placeholder for the main process id of the run.
"combine_ensemble_with_thread_joining": True, # Flag on whether to combine everything and return it to the user or if false: write it to data_dir/ensemble_output_{popuation_id}_{thread_id}.json
- # "output_dir":
"_commandline_input": "",
"log_runtime_systems": 0, # whether to log the runtime of the systems (1 file per thread. stored in the tmp_dir)
"_actually_evolve_system": True, # Whether to actually evolve the systems of just act as if. for testing. used in _process_run_population_grid
- "max_queue_size": 1000, # Maximum size of the system call queue. Can't be too big!
- "_set_ms_grid": False, # Whether the M&S grid has been loaded
+ "max_queue_size": 1000, # Maximum size of the system call queue.
"run_zero_probability_system": True, # Whether to run the zero probability systems
"_zero_prob_stars_skipped": 0,
+ "ensemble_factor_in_probability_weighted_mass": False, # Whether to multiply the ensemble results by 1/probability_weighted_mass
##########################
# Execution log:
##########################
@@ -65,6 +61,12 @@ grid_options_defaults_dict = {
), # TODO: make this more robust
"_binary_c_dir": os.environ["BINARY_C"],
##########################
+ # M&S internal settings
+ ##########################
+ "_loaded_ms_data": False, # Holds flag whether the M&S data is loaded into memory
+ "_set_ms_grid": False, # Whether the M&S grid has been loaded
+ "m&s_options": None, # Holds the M&S options.
+ ##########################
# Custom logging
##########################
"C_auto_logging": None, # Should contain a dictionary where the kes are they headers
@@ -484,7 +486,11 @@ grid_options_descriptions = {
"max_queue_size": "Maximum size of the queue that is used to feed the processes. Don't make this too big! Default: 1000. Input: int",
"_set_ms_grid": "Internal flag whether the M&S grid has been loaded",
"run_zero_probability_system": "Whether to run the zero probability systems. Default: True. Input: boolean",
- "_zero_prob_stars_skipped": "Internal counter to track how many systems are skipped because they have 0 probability"
+ "_zero_prob_stars_skipped": "Internal counter to track how many systems are skipped because they have 0 probability",
+ "ensemble_factor_in_probability_weighted_mass": "Flag to multiply all the ensemble results with 1/probability_weighted_mass",
+ "multiplicity_fraction_function": "Which multiplicity fraction function to use. 0: None, 1: Arenou 2010, 2: Rhagavan 2010, 3: M&S 2017",
+ "m&s_options": "Internal variable that holds the M&S options. Don't write to this your self",
+ "_loaded_ms_data": "Internal variable storing whether the M&S data has been loaded into memory",
}
###
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 8424ee5026658de08bfdbdc2de9017a52d9c8006..cb25228599ceb5b3605588e6989f2041502d7907 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -1,15 +1,18 @@
-# Configuration file for the Sphinx documentation builder.
-#
-# This file only contains a selection of the most common options. For a full
-# list see the documentation:
-# https://www.sphinx-doc.org/en/master/usage/configuration.html
-# https://brendanhasz.github.io/2019/01/05/sphinx.html
+"""
+Configuration file for the Sphinx documentation builder.
+
+This file only contains a selection of the most common options. For a full
+list see the documentation:
+https://www.sphinx-doc.org/en/master/usage/configuration.html
+https://brendanhasz.github.io/2019/01/05/sphinx.html
+"""
+
# -- Path setup --------------------------------------------------------------
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
-#
+
import os
import sys
@@ -127,4 +130,4 @@ print("Done")
from functions import write_custom_footer
print("Generating custom footer")
write_custom_footer()
-print("Done")
\ No newline at end of file
+print("Done")
diff --git a/examples/notebooks/basic_example.ipynb b/examples/notebooks/basic_example.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..ab0bd54d76b898f29e6b0bbb07651c415bd192c9
--- /dev/null
+++ b/examples/notebooks/basic_example.ipynb
@@ -0,0 +1,128 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "66ca8e31-7cbf-4484-bb6c-c9235e14bda2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pandas as pd\n",
+ "import numpy as np\n",
+ "\n",
+ "import binarycpython\n",
+ "from binarycpython.utils.custom_logging_functions import binary_c_log_code\n",
+ "from binarycpython.utils.run_system_wrapper import run_system\n",
+ "from binarycpython.utils.functions import example_parse_output\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "31851c7f-4ccf-4b5a-9f74-8a3e06c6b2d8",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "<AxesSubplot:xlabel='time'>"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbuElEQVR4nO3dfXRc9X3n8fdXj2PLtvwkG2PZSCBjYwzGRhAoCe0GU0zAOOlJTqBNqYNjH5qFpqG7PbBkN+np6YaUnD0pB1rqEwiwB0xch6Y2D6ENhJrNssQ2DQ9+wsYPIAN+BFlI1vN3/5greSwkW9LcmXs19/M6R0d3fvfOzFc/zXz00+/eudfcHRERKXxFURcgIiL5ocAXEUkIBb6ISEIo8EVEEkKBLyKSECVRF3AqkydP9pqamqjLEBEZMTZv3nzY3av6WxfrwK+pqWHTpk1RlyEiMmKY2b6B1mlKR0QkIRT4IiIJocAXEUmIWM/hi4j01dHRQUNDA62trVGXEqlUKkV1dTWlpaWDvk8sA9/MlgBL6urqoi5FRGKmoaGBsWPHUlNTg5lFXU4k3J0jR47Q0NBAbW3toO8Xyykdd1/v7isrKyujLkVEYqa1tZVJkyYlNuwBzIxJkyYN+b+cWAa+iMipJDnsewynDxT4Iqex70gz/2fn4ajLEMmaAl/kNH733pf42kOvRl2GjHAvvfQS119/PQDr1q3jnnvuyXsNsdxpKyIyUrg77k5R0eDHzzfccAM33HBDDqvqn0b4IiJDtHfvXmbPns3NN9/MvHnzWL58OfX19Zx//vl897vf7d3uF7/4BXPmzGHhwoU89dRTve2PPPIIt912GwDLli1j7dq1vevGjBkDwAcffMCVV17JRRddxLx583j55ZezrlsjfBEZsf5q/Ra2vn8s1Mece+Y4vrvk/NNut3PnTh599FEuu+wyjh49ysSJE+nq6uKqq67ijTfe4Nxzz2XFihW8+OKL1NXV8dWvfnVIdTzxxBNcc8013H333XR1ddHS0jLcH6lXLEf4ZrbEzFY1NjZGXYqISL/OOussLrvsMgDWrFnDwoULWbBgAVu2bGHr1q1s376d2tpaZs2ahZnxta99bUiPf8kll/CTn/yE733ve7z55puMHTs265pjOcJ39/XA+vr6+hVR1yIi8TWYkXiuVFRUALBnzx5++MMfsnHjRiZMmMCyZcuGdHx8SUkJ3d3dAHR3d9Pe3g7AlVdeyYYNG3jmmWdYtmwZd9xxBzfffHNWNcdyhC8iMlIcO3aMiooKKisrOXDgAM899xwAc+bMYe/evbzzzjsArF69ut/719TUsHnzZiB99E5HRwcA+/btY+rUqaxYsYJvfOMbvPbaa1nXGssRvojISDF//nwWLFjAnDlzmDFjBldccQWQPtfNqlWruO666xg9ejSf+9znaGpq+tT9V6xYwdKlS5k/fz6LFy/u/c/hpZde4t5776W0tJQxY8bw2GOPZV2ruXvWD5Ir9fX1rgugSNRq7nwGgL33XBdxJQKwbds2zjvvvKjLiIX++sLMNrt7fX/ba0pHRCQhFPgiIgmhwBeRESfOU9H5Mpw+UOCLyIiSSqU4cuRIokO/53z4qVRqSPfTUToiMqJUV1fT0NDAoUOHoi4lUj1XvBoKBb7IAA42tVI5avCXj5P8KC0tHdJVnuSEvAW+mX0RuA4YBzzk7v+ar+cWGY5L/+YFrjl/atRliIQmqzl8M3vYzA6a2Vt92heb2Q4z22VmdwK4+8/dfQVwKzC0swiJROT5LQeiLkEkNNnutH0EWJzZYGbFwAPAtcBc4CYzm5uxyXeC9SKxleQdglK4sgp8d98AHO3TfCmwy913u3s78CSw1NJ+ADzn7gOeFMLMVprZJjPblPSdMhId5b0UolwcljkdeC/jdkPQdjuwCPiymd060J3dfZW717t7fVVVVQ7KEzm97iDxda1sKSR522nr7vcB9w1mWzNbAiypq6vLbVEiA+gORvgGaLAvhSIXI/z9wIyM29VB26C5+3p3X1lZWRlqYSKD5UHMF2mILwUkF4G/EZhlZrVmVgbcCKwbygPoilcStZ45fAW+FJJsD8tcDbwCzDazBjNb7u6dwG3A88A2YI27bxnK42qEL1HTHL4Uoqzm8N39pgHanwWezeaxRaLUrRG+FKBYnjxNUzoStZ4RfpHyXgpILANfUzoSNU9fUxrTCF8KSCwDXyRqmsOXQhTLwNeUjkSt59h7zeFLIYll4GtKR6KmOXwpRLEMfJGonQh8Jb4UDgW+SD96PnilnbZSSGIZ+JrDl6hpSkcKUSwDX3P4ErXek6cp8KWAxDLwRaLmmsOXAqTAF+mHTp4mhSiWga85fIlaZzCnUxTLd4jI8MTy5aw5fIlaZ1f63AqlSnwpIHo1i/Sjoys9wi8p1pSOFA4Fvkg/OrvTI/wSjfClgOjVLNKPnhF+qUb4UkAU+CL96JnDLynWW0QKRyxfzTpKR6LWc5ROiT5qKwUkloGvo3Qkah09R+lohC8FRK9mkX506igdKUAKfJF+6CgdKUR6NYv0Q0fpSCFS4Iv0o3eErzl8KSB6NYv0o70zHfjlJXqLSOHQq1mkH8fbuwAYVVoccSUi4Yll4Os4fIlaS0cQ+GUKfCkcsQx8HYcvUWsNRvia0pFColezSD+OByN8HaMjhUSBL9KPnsBvCz5xK1IIFPgi/Tjeng76o5+0R1yJSHgU+CL9aA1G+B8ea424EpHwKPBF+tHc3gnAu0dbIq5EJDwKfJF+NB7vAGDfkXTgV+jwTCkACnyRfvQEfg8djy+FQIEv0o9jfQJfpBAo8EX6cPeTRvjVE0ZFWI1IePIW+GZ2tpk9ZGZr8/WcIsNxvKOLji5n8pgy7v3yhfzuuVVRlyQSiqwC38weNrODZvZWn/bFZrbDzHaZ2Z0A7r7b3Zdn83wi+fBxS3p0/xe/P5uv1M+IuBqR8GQ7wn8EWJzZYGbFwAPAtcBc4CYzm5vl84jkzcGmNgCmjC2PuBKRcGUV+O6+ATjap/lSYFcwom8HngSWDvYxzWylmW0ys02HDh3KpjyRYTkQfNhq6rhUxJWIhCsXc/jTgfcybjcA081skpk9CCwws7sGurO7r3L3enevr6rS3Knk38Eg8KeM0whfCktJvp7I3Y8Atw5mWzNbAiypq6vLbVEi/ThwrI3iImNShQJfCksuRvj7gcw9XdVB26DpfPgSpQ+PtVI1ppziIp0cWQpLLgJ/IzDLzGrNrAy4EVg3lAfQFa8kSu8eaWHmxNFRlyESumwPy1wNvALMNrMGM1vu7p3AbcDzwDZgjbtvGcrjaoQvUdpzpJmayQp8KTxZzeG7+00DtD8LPJvNY4tE4ZO2Tg41tVEzuSLqUkRCF8tTK2hKR6Ky93AzALWTFPhSeGIZ+JrSkai8faAJgLopYyKuRCR8sQx8kai8ub+RUaXFnF2lwJfCE8vA15SOROWt/Y3MPXOcDsmUghTLwNeUjkShq9vZ+v4xLpiu150UplgGvkgUtrzfSHN7Fwtmjo+6FJGcUOCLBH696wgAv3PO5IgrEcmNWAa+5vAlCv/3ncOcO3UMVTotshSoWAa+5vAl3z5p6+TVPUf5bJ3O0CqFK5aBL5JvL2w7QHtnN9decEbUpYjkjAJfBHjmjQ+YOq6ci2dOiLoUkZyJZeBrDl/y6fAnbfxqx0Guu+BMinT8vRSwWAa+5vAln9Zseo+OLucPPzMz6lJEciqWgS+SLx1d3Tz+/97l8rMn6fw5UvAU+JJo//zafvZ/fJzln62NuhSRnFPgS2J1dHVz34s7ubC6kqvOmxJ1OSI5p8CXxFq1YTcNHx3n21efi5l21krhi2Xg6ygdybWt7x/jR798m+svnMZ/mq3RvSRDLANfR+lILrV1dnHHmt9SOaqMv146L+pyRPImq2vaioxE972wk+0fNvHjm+uZUFEWdTkieRPLEb5Irmze9xH/8NI7fOXiahbNnRp1OSJ5pcCXxDja3M7tT7zGtMpR/Pclc6MuRyTvNKUjidDV7Xzryf/g8CftrP3TyxmXKo26JJG8U+BLwXN3/vrprby88zD/80sXcGH1+KhLEomEAl8KWmtHF9/5+Vus3dzALVfUctOlM6IuSSQysQx8M1sCLKmrq4u6FBnB3mxo5C9/9gbbPjjGn101i28vmqUPWEmixTLw3X09sL6+vn5F1LXIyLP3cDP3/2oXP3utgUkV5Ty8rJ7Pz9EROSKxDHyRoero6uaFbQd4/NV3eXnnYcpLilh+RS1/tmiWdtCKBBT4MmI1tXbw728f4t+2HuDF7Qdpau1kWmWKP180iz+8dCZTxqWiLlEkVhT4MmK0dnTx2rsf8Zs9R3l191E27TtKR5czsaKMxeefwbUXnMGVs6ooKdbHS0T6o8CXWGtq7eCRX+9lw85DvP5eI+1d3ZjB3Gnj+PoVtSw6byoXnzWBYl2aUOS0FPgSK40tHew40MSOA01sff8Yv9x2gENNbcyfMZ6vX1HDZ86eyMVnTaRylOblRYZKgS+R+Ki5nT1Hmtl9qJm3DzSx/cMm3v6wiQ+PtfZuMy5VwsKzJvDtRecyf8b46IoVKRAKfMkJd+doczvvfXScvYeb2XO4mb1Hmtl7pIW9h5tpPN7Ru21ZSRF1VWO4/JxJzD5jLLOnjmX2GWOZVpnScfMiIVLgy7C0d3bzYWMr+z8+zv6Pj/N+8JV5u7Wju3d7MzizchQ1k0dz/YXTqJ1cwVmTKqidXEHNpNHa0SqSBwp8OUlLeycHj7VxsKmNg02tJy0famoLbrfyUUvHp+47eUw508enmHPGWD4/ewrTJ4xi+vhR1E6uYMbE0aRKiyP4iUSkR94C38wqgL8H2oGX3P3xfD13UnV1O43HOzja3N779VFLxnJzO0f73G5u7/rU45QWG1Vjyqkal2LmpNHU10xgytgU08anmD5+FGeOH8W0ypQCXSTmsgp8M3sYuB446O7zMtoXA38HFAM/dvd7gD8A1rr7ejP7KaDAP42Orm6a2zppau3kWGsHx473fO/gWGtn8D3d3ti73JHe/ngHTW2dAz52RVkxEyrKmBh81VWNYUJFGZPGlDFlbIopY8uZMq6cKWNTjB9VSpEOexQZ8bId4T8C3A881tNgZsXAA8DVQAOw0czWAdXAm8Fmnx5GjlDd3U5rZxetHd20dnQFX91BWxctbV00t3fS3NZFS/A9fbuTlvYuPmnr7G3PXN/S1kV7V/dpn39sqoRxqVLGjSplXKqEGRNHB7dLGJsqZeLo0pOCfWJFGRNGl2k0LpJAWQW+u28ws5o+zZcCu9x9N4CZPQksJR3+1cBvOcWVtsxsJbASYObMmcOqq/F4B20dXbR1dtN2UhgH3zu7ON7eRWtnN22ZIR2sy9y+rXPg9W0d3YMK5b5SpUVUlJVQUV7C6LJiKspLGJsqYVplitFlJVSUp9sqyooZXVbCmFQJlaNKe4O8J+DHlJfoA0ciMmi5mMOfDryXcbsB+AxwH3C/mV0HrB/ozu6+ClgFUF9f78Mp4EsP/Jrdh5uHdB8zSJUUkyotIlVaTKq0mPKSnuUiJo0p+9T6nnWp0mJSJUUntZWXFpMqKc4I7xJGlxdTUaaQFpFo5G2nrbs3A18fzLbZng//9qvqaGnvorykmLKSok+FcSoI495gLi2irLhIx3yLSEHLReDvBzIvK1QdtA1atufD/9KC6uHcTUSkoOXi0y4bgVlmVmtmZcCNwLocPI+IiAxBVoFvZquBV4DZZtZgZsvdvRO4DXge2AascfctQ3zcJWa2qrGxMZvyREQkQ7ZH6dw0QPuzwLNZPK4ucSgiErJYnsBEI3wRkfDFMvDdfb27r6ysrIy6FBGRghHLwBcRkfDFMvA1pSMiEr5YBr6mdEREwhfLwBcRkfAp8EVEEiKWga85fBGR8MUy8DWHLyISvlgGvoiIhE+BLyKSELEMfM3hi4iEL5aBrzl8EZHwxTLwRUQkfAp8EZGEUOCLiCSEAl9EJCFiGfg6SkdEJHyxDHwdpSMiEr5YBr6IiIRPgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgkRy8DXcfgiIuGLZeDrOHwRkfDFMvBFRCR8CnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEHkLfDM728weMrO1+XpOERE5YVCBb2YPm9lBM3urT/tiM9thZrvM7M5TPYa773b35dkUKyIiw1cyyO0eAe4HHutpMLNi4AHgaqAB2Ghm64Bi4Pt97n+Lux/MuloRERm2QQW+u28ws5o+zZcCu9x9N4CZPQksdffvA9cPtyAzWwmsBJg5c+ZwH0ZERPrIZg5/OvBexu2GoK1fZjbJzB4EFpjZXQNt5+6r3L3e3eurqqqyKE9ERDINdkona+5+BLh1MNua2RJgSV1dXW6LEhFJkGxG+PuBGRm3q4O2rOl8+CIi4csm8DcCs8ys1szKgBuBdeGUJSIiYRvsYZmrgVeA2WbWYGbL3b0TuA14HtgGrHH3LWEUpUscioiEb7BH6dw0QPuzwLOhVpR+3PXA+vr6+hVhP7aISFLF8tQKGuGLiIQvloGvnbYiIuGLZeCLiEj4Yhn4mtIREQlfLANfUzoiIuGLZeCLiEj4FPgiIgkRy8DXHL6ISPhiGfiawxcRCV8sA19ERMKnwBcRSQgFvohIQsQy8LXTVkQkfLEMfO20FREJXywDX0REwqfAFxFJCAW+iEhCxDLwtdNWRCR8sQx87bQVEQlfLANfRETCp8AXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEiGXg6zh8EZHwxTLwdRy+iEj4Yhn4IiISPgW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhSqIuQGQkONrcztX/69+jLkMS5NlvfY7S4nDH5Ap8kdP44oLpfNzSgeNRlyKSlbwFvpl9EbgOGAc85O7/mq/nFsnGJTUTuaRmYtRliGRtUP8vmNnDZnbQzN7q077YzHaY2S4zu/NUj+HuP3f3FcCtwFeHX7KIiAzHYEf4jwD3A4/1NJhZMfAAcDXQAGw0s3VAMfD9Pve/xd0PBsvfCe4nIiJ5NKjAd/cNZlbTp/lSYJe77wYwsyeBpe7+feD6vo9hZgbcAzzn7q8N9FxmthJYCTBz5szBlCciIoOQzS7g6cB7GbcbgraB3A4sAr5sZrcOtJG7r3L3enevr6qqyqI8ERHJlLedtu5+H3DfYLY1syXAkrq6utwWJSKSINmM8PcDMzJuVwdtWdP58EVEwpdN4G8EZplZrZmVATcC68IpS0REwjbYwzJXA68As82swcyWu3sncBvwPLANWOPuW8IoSpc4FBEJn7nH99ODZnYI2DfMu08GDodYTtjiXh+oxjDEvT6If41xrw/iVeNZ7t7vES+xDvxsmNkmd6+Puo6BxL0+UI1hiHt9EP8a414fjIwaQWfLFBFJDAW+iEhCFHLgr4q6gNOIe32gGsMQ9/og/jXGvT4YGTUW7hy+iIicrJBH+CIikkGBLyKSEAUX+EM5R38OnnuGmf3KzLaa2RYz+1bQPtHM/s3MdgbfJwTtZmb3BbW+YWYLMx7rT4Ltd5rZn4RcZ7GZ/YeZPR3crjWzV4M6fhp8chozKw9u7wrW12Q8xl1B+w4zuybk+sab2Voz225m28zs8jj1oZl9O/j9vmVmq80sFXUf9nfNijD7zMwuNrM3g/vcZ2YWUo33Br/nN8zsn81sfMa6fvtnoPf4QL+DbOrLWPcXZuZmNjm4HUkfZs3dC+aL9Ln43wHOBsqA14G5eXz+acDCYHks8DYwF/hb4M6g/U7gB8HyF4DnAAMuA14N2icCu4PvE4LlCSHWeQfwBPB0cHsNcGOw/CDwp8HyN4EHg+UbgZ8Gy3ODvi0HaoM+Lw6xvkeBbwTLZcD4uPQh6TPC7gFGZfTdsqj7ELgSWAi8ldEWWp8Bvwm2teC+14ZU4+8DJcHyDzJq7Ld/OMV7fKDfQTb1Be0zSJ9RYB8wOco+zPr1m+8nzOkPA5cDz2fcvgu4K8J6/oX0BWJ2ANOCtmnAjmD5H4GbMrbfEay/CfjHjPaTtsuypmrgBeDzwNPBi+9wxpuutw+DF/nlwXJJsJ317dfM7UKor5J0oFqf9lj0ISdOCz4x6JOngWvi0IdADSeHaSh9FqzbntF+0nbZ1Nhn3ZeAx4PlfvuHAd7jp3odZ1sfsBaYD+zlROBH1ofZfBXalM5Qz9GfM8G/7guAV4Gp7v5BsOpDYGqwPFC9ufw5fgT8JdAd3J4EfOzpcyP1fa7eOoL1jcH2uayvFjgE/MTS004/NrMKYtKH7r4f+CHwLvAB6T7ZTLz6sEdYfTY9WM5lrQC3kB75DqfGU72Oh83MlgL73f31Pqvi2oenVGiBHwtmNgb4GfDn7n4sc52n/7xHciysmV0PHHT3zVE8/yCVkP63+h/cfQHQTHo6olfEfTgBWEr6D9OZQAWwOIpahiLKPhsMM7sb6AQej7qWHmY2GvhvwP+IupawFFrg5+wc/YNlZqWkw/5xd38qaD5gZtOC9dOAnuv7DlRvrn6OK4AbzGwv8CTpaZ2/A8abWc/FcDKfq7eOYH0lcCSH9UF65NPg7q8Gt9eS/gMQlz5cBOxx90Pu3gE8Rbpf49SHPcLqs/3Bck5qNbNlpC+L+kfBH6bh1HiEgX8Hw3UO6T/srwfvmWrgNTM7Yxj15bQPBy3fc0i5/CI9OtxN+pfUs0Pn/Dw+v5G+0PuP+rTfy8k7z/42WL6Ok3f8/CZon0h6HntC8LUHmBhyrb/HiZ22/8TJO7u+GSz/Z07e4bgmWD6fk3eo7SbcnbYvA7OD5e8F/ReLPgQ+A2wBRgfP+Sjpy3dG3od8eg4/tD7j0zscvxBSjYuBrUBVn+367R9O8R4f6HeQTX191u3lxBx+ZH2Y1Wsk30+Y8x8ovff8bdJ78u/O83N/lvS/zW8Avw2+vkB6fvEFYCfwy4wXgAEPBLW+CdRnPNYtwK7g6+s5qPX3OBH4Zwcvxl3Bm6Y8aE8Ft3cF68/OuP/dQd07CPloA+AiYFPQjz8P3jix6UPgr4DtwFvA/w5CKdI+BFaT3qfQQfq/pOVh9hlQH/y87wD302enehY17iI9593zfnnwdP3DAO/xgX4H2dTXZ/1eTgR+JH2Y7ZdOrSAikhCFNocvIiIDUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+SMDSZ+n8ZrB8ppmtjbomkTDpsEyRQHD+o6fdfV7UtYjkQsnpNxFJjHuAc8zst6Q/rHSeu88LPvr/RdLnzZlF+uRpZcAfA22kPzF51MzOIf1hnCqgBVjh7tvz/UOIDERTOiIn3Am84+4XAf+1z7p5wB8AlwB/A7R4+uRurwA3B9usAm5394uB/wL8fT6KFhksjfBFBudX7t4ENJlZI7A+aH8TuDA4Q+rvAP+UcSGj8vyXKTIwBb7I4LRlLHdn3O4m/T4qIn0+9ovyXJfIoGlKR+SEJtKXphwyT1/3YI+ZfQV6r3k6P8ziRLKlwBcJuPsR4NfBRazvHcZD/BGw3MxeJ30K5aVh1ieSLR2WKSKSEBrhi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQ/x8q6Ro0odz3FAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# You can decide to `write` your own logging_line, which allows you to write a more complex logging statement with conditionals.\n",
+ "logging_line = 'Printf(\"MY_STELLAR_DATA time=%g radius=%g\\\\n\",stardata->model.time,stardata->star[0].radius)'\n",
+ "\n",
+ "# Generate entire shared lib code around logging lines\n",
+ "custom_logging_code = binary_c_log_code(logging_line, verbose=-1)\n",
+ "\n",
+ "# Run system. all arguments can be given as optional arguments. the custom_logging_code is \n",
+ "# one of them and will be processed automatically.\n",
+ "output = run_system(\n",
+ " M_1=1, metallicity=0.002, M_2=0.1, separation=0, orbital_period=100000000000, \\\n",
+ " custom_logging_code=custom_logging_code,\n",
+ ")\n",
+ "\n",
+ "df = pd.DataFrame.from_dict(example_parse_output(output, 'MY_STELLAR_DATA'), dtype=np.float64)\n",
+ "\n",
+ "df.plot('time', 'radius', logy=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "db0f3990-8680-4ea7-80d3-cdb430950c7f",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "<AxesSubplot:xlabel='time'>"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbuElEQVR4nO3dfXRc9X3n8fdXj2PLtvwkG2PZSCBjYwzGRhAoCe0GU0zAOOlJTqBNqYNjH5qFpqG7PbBkN+np6YaUnD0pB1rqEwiwB0xch6Y2D6ENhJrNssQ2DQ9+wsYPIAN+BFlI1vN3/5greSwkW9LcmXs19/M6R0d3fvfOzFc/zXz00+/eudfcHRERKXxFURcgIiL5ocAXEUkIBb6ISEIo8EVEEkKBLyKSECVRF3AqkydP9pqamqjLEBEZMTZv3nzY3av6WxfrwK+pqWHTpk1RlyEiMmKY2b6B1mlKR0QkIRT4IiIJocAXEUmIWM/hi4j01dHRQUNDA62trVGXEqlUKkV1dTWlpaWDvk8sA9/MlgBL6urqoi5FRGKmoaGBsWPHUlNTg5lFXU4k3J0jR47Q0NBAbW3toO8Xyykdd1/v7isrKyujLkVEYqa1tZVJkyYlNuwBzIxJkyYN+b+cWAa+iMipJDnsewynDxT4Iqex70gz/2fn4ajLEMmaAl/kNH733pf42kOvRl2GjHAvvfQS119/PQDr1q3jnnvuyXsNsdxpKyIyUrg77k5R0eDHzzfccAM33HBDDqvqn0b4IiJDtHfvXmbPns3NN9/MvHnzWL58OfX19Zx//vl897vf7d3uF7/4BXPmzGHhwoU89dRTve2PPPIIt912GwDLli1j7dq1vevGjBkDwAcffMCVV17JRRddxLx583j55ZezrlsjfBEZsf5q/Ra2vn8s1Mece+Y4vrvk/NNut3PnTh599FEuu+wyjh49ysSJE+nq6uKqq67ijTfe4Nxzz2XFihW8+OKL1NXV8dWvfnVIdTzxxBNcc8013H333XR1ddHS0jLcH6lXLEf4ZrbEzFY1NjZGXYqISL/OOussLrvsMgDWrFnDwoULWbBgAVu2bGHr1q1s376d2tpaZs2ahZnxta99bUiPf8kll/CTn/yE733ve7z55puMHTs265pjOcJ39/XA+vr6+hVR1yIi8TWYkXiuVFRUALBnzx5++MMfsnHjRiZMmMCyZcuGdHx8SUkJ3d3dAHR3d9Pe3g7AlVdeyYYNG3jmmWdYtmwZd9xxBzfffHNWNcdyhC8iMlIcO3aMiooKKisrOXDgAM899xwAc+bMYe/evbzzzjsArF69ut/719TUsHnzZiB99E5HRwcA+/btY+rUqaxYsYJvfOMbvPbaa1nXGssRvojISDF//nwWLFjAnDlzmDFjBldccQWQPtfNqlWruO666xg9ejSf+9znaGpq+tT9V6xYwdKlS5k/fz6LFy/u/c/hpZde4t5776W0tJQxY8bw2GOPZV2ruXvWD5Ir9fX1rgugSNRq7nwGgL33XBdxJQKwbds2zjvvvKjLiIX++sLMNrt7fX/ba0pHRCQhFPgiIgmhwBeRESfOU9H5Mpw+UOCLyIiSSqU4cuRIokO/53z4qVRqSPfTUToiMqJUV1fT0NDAoUOHoi4lUj1XvBoKBb7IAA42tVI5avCXj5P8KC0tHdJVnuSEvAW+mX0RuA4YBzzk7v+ar+cWGY5L/+YFrjl/atRliIQmqzl8M3vYzA6a2Vt92heb2Q4z22VmdwK4+8/dfQVwKzC0swiJROT5LQeiLkEkNNnutH0EWJzZYGbFwAPAtcBc4CYzm5uxyXeC9SKxleQdglK4sgp8d98AHO3TfCmwy913u3s78CSw1NJ+ADzn7gOeFMLMVprZJjPblPSdMhId5b0UolwcljkdeC/jdkPQdjuwCPiymd060J3dfZW717t7fVVVVQ7KEzm97iDxda1sKSR522nr7vcB9w1mWzNbAiypq6vLbVEiA+gORvgGaLAvhSIXI/z9wIyM29VB26C5+3p3X1lZWRlqYSKD5UHMF2mILwUkF4G/EZhlZrVmVgbcCKwbygPoilcStZ45fAW+FJJsD8tcDbwCzDazBjNb7u6dwG3A88A2YI27bxnK42qEL1HTHL4Uoqzm8N39pgHanwWezeaxRaLUrRG+FKBYnjxNUzoStZ4RfpHyXgpILANfUzoSNU9fUxrTCF8KSCwDXyRqmsOXQhTLwNeUjkSt59h7zeFLIYll4GtKR6KmOXwpRLEMfJGonQh8Jb4UDgW+SD96PnilnbZSSGIZ+JrDl6hpSkcKUSwDX3P4ErXek6cp8KWAxDLwRaLmmsOXAqTAF+mHTp4mhSiWga85fIlaZzCnUxTLd4jI8MTy5aw5fIlaZ1f63AqlSnwpIHo1i/Sjoys9wi8p1pSOFA4Fvkg/OrvTI/wSjfClgOjVLNKPnhF+qUb4UkAU+CL96JnDLynWW0QKRyxfzTpKR6LWc5ROiT5qKwUkloGvo3Qkah09R+lohC8FRK9mkX506igdKUAKfJF+6CgdKUR6NYv0Q0fpSCFS4Iv0o3eErzl8KSB6NYv0o70zHfjlJXqLSOHQq1mkH8fbuwAYVVoccSUi4Yll4Os4fIlaS0cQ+GUKfCkcsQx8HYcvUWsNRvia0pFColezSD+OByN8HaMjhUSBL9KPnsBvCz5xK1IIFPgi/Tjeng76o5+0R1yJSHgU+CL9aA1G+B8ea424EpHwKPBF+tHc3gnAu0dbIq5EJDwKfJF+NB7vAGDfkXTgV+jwTCkACnyRfvQEfg8djy+FQIEv0o9jfQJfpBAo8EX6cPeTRvjVE0ZFWI1IePIW+GZ2tpk9ZGZr8/WcIsNxvKOLji5n8pgy7v3yhfzuuVVRlyQSiqwC38weNrODZvZWn/bFZrbDzHaZ2Z0A7r7b3Zdn83wi+fBxS3p0/xe/P5uv1M+IuBqR8GQ7wn8EWJzZYGbFwAPAtcBc4CYzm5vl84jkzcGmNgCmjC2PuBKRcGUV+O6+ATjap/lSYFcwom8HngSWDvYxzWylmW0ys02HDh3KpjyRYTkQfNhq6rhUxJWIhCsXc/jTgfcybjcA081skpk9CCwws7sGurO7r3L3enevr6rS3Knk38Eg8KeM0whfCktJvp7I3Y8Atw5mWzNbAiypq6vLbVEi/ThwrI3iImNShQJfCksuRvj7gcw9XdVB26DpfPgSpQ+PtVI1ppziIp0cWQpLLgJ/IzDLzGrNrAy4EVg3lAfQFa8kSu8eaWHmxNFRlyESumwPy1wNvALMNrMGM1vu7p3AbcDzwDZgjbtvGcrjaoQvUdpzpJmayQp8KTxZzeG7+00DtD8LPJvNY4tE4ZO2Tg41tVEzuSLqUkRCF8tTK2hKR6Ky93AzALWTFPhSeGIZ+JrSkai8faAJgLopYyKuRCR8sQx8kai8ub+RUaXFnF2lwJfCE8vA15SOROWt/Y3MPXOcDsmUghTLwNeUjkShq9vZ+v4xLpiu150UplgGvkgUtrzfSHN7Fwtmjo+6FJGcUOCLBH696wgAv3PO5IgrEcmNWAa+5vAlCv/3ncOcO3UMVTotshSoWAa+5vAl3z5p6+TVPUf5bJ3O0CqFK5aBL5JvL2w7QHtnN9decEbUpYjkjAJfBHjmjQ+YOq6ci2dOiLoUkZyJZeBrDl/y6fAnbfxqx0Guu+BMinT8vRSwWAa+5vAln9Zseo+OLucPPzMz6lJEciqWgS+SLx1d3Tz+/97l8rMn6fw5UvAU+JJo//zafvZ/fJzln62NuhSRnFPgS2J1dHVz34s7ubC6kqvOmxJ1OSI5p8CXxFq1YTcNHx3n21efi5l21krhi2Xg6ygdybWt7x/jR798m+svnMZ/mq3RvSRDLANfR+lILrV1dnHHmt9SOaqMv146L+pyRPImq2vaioxE972wk+0fNvHjm+uZUFEWdTkieRPLEb5Irmze9xH/8NI7fOXiahbNnRp1OSJ5pcCXxDja3M7tT7zGtMpR/Pclc6MuRyTvNKUjidDV7Xzryf/g8CftrP3TyxmXKo26JJG8U+BLwXN3/vrprby88zD/80sXcGH1+KhLEomEAl8KWmtHF9/5+Vus3dzALVfUctOlM6IuSSQysQx8M1sCLKmrq4u6FBnB3mxo5C9/9gbbPjjGn101i28vmqUPWEmixTLw3X09sL6+vn5F1LXIyLP3cDP3/2oXP3utgUkV5Ty8rJ7Pz9EROSKxDHyRoero6uaFbQd4/NV3eXnnYcpLilh+RS1/tmiWdtCKBBT4MmI1tXbw728f4t+2HuDF7Qdpau1kWmWKP180iz+8dCZTxqWiLlEkVhT4MmK0dnTx2rsf8Zs9R3l191E27TtKR5czsaKMxeefwbUXnMGVs6ooKdbHS0T6o8CXWGtq7eCRX+9lw85DvP5eI+1d3ZjB3Gnj+PoVtSw6byoXnzWBYl2aUOS0FPgSK40tHew40MSOA01sff8Yv9x2gENNbcyfMZ6vX1HDZ86eyMVnTaRylOblRYZKgS+R+Ki5nT1Hmtl9qJm3DzSx/cMm3v6wiQ+PtfZuMy5VwsKzJvDtRecyf8b46IoVKRAKfMkJd+doczvvfXScvYeb2XO4mb1Hmtl7pIW9h5tpPN7Ru21ZSRF1VWO4/JxJzD5jLLOnjmX2GWOZVpnScfMiIVLgy7C0d3bzYWMr+z8+zv6Pj/N+8JV5u7Wju3d7MzizchQ1k0dz/YXTqJ1cwVmTKqidXEHNpNHa0SqSBwp8OUlLeycHj7VxsKmNg02tJy0famoLbrfyUUvHp+47eUw508enmHPGWD4/ewrTJ4xi+vhR1E6uYMbE0aRKiyP4iUSkR94C38wqgL8H2oGX3P3xfD13UnV1O43HOzja3N779VFLxnJzO0f73G5u7/rU45QWG1Vjyqkal2LmpNHU10xgytgU08anmD5+FGeOH8W0ypQCXSTmsgp8M3sYuB446O7zMtoXA38HFAM/dvd7gD8A1rr7ejP7KaDAP42Orm6a2zppau3kWGsHx473fO/gWGtn8D3d3ti73JHe/ngHTW2dAz52RVkxEyrKmBh81VWNYUJFGZPGlDFlbIopY8uZMq6cKWNTjB9VSpEOexQZ8bId4T8C3A881tNgZsXAA8DVQAOw0czWAdXAm8Fmnx5GjlDd3U5rZxetHd20dnQFX91BWxctbV00t3fS3NZFS/A9fbuTlvYuPmnr7G3PXN/S1kV7V/dpn39sqoRxqVLGjSplXKqEGRNHB7dLGJsqZeLo0pOCfWJFGRNGl2k0LpJAWQW+u28ws5o+zZcCu9x9N4CZPQksJR3+1cBvOcWVtsxsJbASYObMmcOqq/F4B20dXbR1dtN2UhgH3zu7ON7eRWtnN22ZIR2sy9y+rXPg9W0d3YMK5b5SpUVUlJVQUV7C6LJiKspLGJsqYVplitFlJVSUp9sqyooZXVbCmFQJlaNKe4O8J+DHlJfoA0ciMmi5mMOfDryXcbsB+AxwH3C/mV0HrB/ozu6+ClgFUF9f78Mp4EsP/Jrdh5uHdB8zSJUUkyotIlVaTKq0mPKSnuUiJo0p+9T6nnWp0mJSJUUntZWXFpMqKc4I7xJGlxdTUaaQFpFo5G2nrbs3A18fzLbZng//9qvqaGnvorykmLKSok+FcSoI495gLi2irLhIx3yLSEHLReDvBzIvK1QdtA1atufD/9KC6uHcTUSkoOXi0y4bgVlmVmtmZcCNwLocPI+IiAxBVoFvZquBV4DZZtZgZsvdvRO4DXge2AascfctQ3zcJWa2qrGxMZvyREQkQ7ZH6dw0QPuzwLNZPK4ucSgiErJYnsBEI3wRkfDFMvDdfb27r6ysrIy6FBGRghHLwBcRkfDFMvA1pSMiEr5YBr6mdEREwhfLwBcRkfAp8EVEEiKWga85fBGR8MUy8DWHLyISvlgGvoiIhE+BLyKSELEMfM3hi4iEL5aBrzl8EZHwxTLwRUQkfAp8EZGEUOCLiCSEAl9EJCFiGfg6SkdEJHyxDHwdpSMiEr5YBr6IiIRPgS8ikhAKfBGRhFDgi4gkhAJfRCQhFPgiIgkRy8DXcfgiIuGLZeDrOHwRkfDFMvBFRCR8CnwRkYRQ4IuIJIQCX0QkIRT4IiIJocAXEUkIBb6ISEIo8EVEEkKBLyKSEHkLfDM728weMrO1+XpOERE5YVCBb2YPm9lBM3urT/tiM9thZrvM7M5TPYa773b35dkUKyIiw1cyyO0eAe4HHutpMLNi4AHgaqAB2Ghm64Bi4Pt97n+Lux/MuloRERm2QQW+u28ws5o+zZcCu9x9N4CZPQksdffvA9cPtyAzWwmsBJg5c+ZwH0ZERPrIZg5/OvBexu2GoK1fZjbJzB4EFpjZXQNt5+6r3L3e3eurqqqyKE9ERDINdkona+5+BLh1MNua2RJgSV1dXW6LEhFJkGxG+PuBGRm3q4O2rOl8+CIi4csm8DcCs8ys1szKgBuBdeGUJSIiYRvsYZmrgVeA2WbWYGbL3b0TuA14HtgGrHH3LWEUpUscioiEb7BH6dw0QPuzwLOhVpR+3PXA+vr6+hVhP7aISFLF8tQKGuGLiIQvloGvnbYiIuGLZeCLiEj4Yhn4mtIREQlfLANfUzoiIuGLZeCLiEj4FPgiIgkRy8DXHL6ISPhiGfiawxcRCV8sA19ERMKnwBcRSQgFvohIQsQy8LXTVkQkfLEMfO20FREJXywDX0REwqfAFxFJCAW+iEhCxDLwtdNWRCR8sQx87bQVEQlfLANfRETCp8AXEUkIBb6ISEIo8EVEEkKBLyKSEAp8EZGEiGXg6zh8EZHwxTLwdRy+iEj4Yhn4IiISPgW+iEhCKPBFRBJCgS8ikhAKfBGRhFDgi4gkhAJfRCQhSqIuQGQkONrcztX/69+jLkMS5NlvfY7S4nDH5Ap8kdP44oLpfNzSgeNRlyKSlbwFvpl9EbgOGAc85O7/mq/nFsnGJTUTuaRmYtRliGRtUP8vmNnDZnbQzN7q077YzHaY2S4zu/NUj+HuP3f3FcCtwFeHX7KIiAzHYEf4jwD3A4/1NJhZMfAAcDXQAGw0s3VAMfD9Pve/xd0PBsvfCe4nIiJ5NKjAd/cNZlbTp/lSYJe77wYwsyeBpe7+feD6vo9hZgbcAzzn7q8N9FxmthJYCTBz5szBlCciIoOQzS7g6cB7GbcbgraB3A4sAr5sZrcOtJG7r3L3enevr6qqyqI8ERHJlLedtu5+H3DfYLY1syXAkrq6utwWJSKSINmM8PcDMzJuVwdtWdP58EVEwpdN4G8EZplZrZmVATcC68IpS0REwjbYwzJXA68As82swcyWu3sncBvwPLANWOPuW8IoSpc4FBEJn7nH99ODZnYI2DfMu08GDodYTtjiXh+oxjDEvT6If41xrw/iVeNZ7t7vES+xDvxsmNkmd6+Puo6BxL0+UI1hiHt9EP8a414fjIwaQWfLFBFJDAW+iEhCFHLgr4q6gNOIe32gGsMQ9/og/jXGvT4YGTUW7hy+iIicrJBH+CIikkGBLyKSEAUX+EM5R38OnnuGmf3KzLaa2RYz+1bQPtHM/s3MdgbfJwTtZmb3BbW+YWYLMx7rT4Ltd5rZn4RcZ7GZ/YeZPR3crjWzV4M6fhp8chozKw9u7wrW12Q8xl1B+w4zuybk+sab2Voz225m28zs8jj1oZl9O/j9vmVmq80sFXUf9nfNijD7zMwuNrM3g/vcZ2YWUo33Br/nN8zsn81sfMa6fvtnoPf4QL+DbOrLWPcXZuZmNjm4HUkfZs3dC+aL9Ln43wHOBsqA14G5eXz+acDCYHks8DYwF/hb4M6g/U7gB8HyF4DnAAMuA14N2icCu4PvE4LlCSHWeQfwBPB0cHsNcGOw/CDwp8HyN4EHg+UbgZ8Gy3ODvi0HaoM+Lw6xvkeBbwTLZcD4uPQh6TPC7gFGZfTdsqj7ELgSWAi8ldEWWp8Bvwm2teC+14ZU4+8DJcHyDzJq7Ld/OMV7fKDfQTb1Be0zSJ9RYB8wOco+zPr1m+8nzOkPA5cDz2fcvgu4K8J6/oX0BWJ2ANOCtmnAjmD5H4GbMrbfEay/CfjHjPaTtsuypmrgBeDzwNPBi+9wxpuutw+DF/nlwXJJsJ317dfM7UKor5J0oFqf9lj0ISdOCz4x6JOngWvi0IdADSeHaSh9FqzbntF+0nbZ1Nhn3ZeAx4PlfvuHAd7jp3odZ1sfsBaYD+zlROBH1ofZfBXalM5Qz9GfM8G/7guAV4Gp7v5BsOpDYGqwPFC9ufw5fgT8JdAd3J4EfOzpcyP1fa7eOoL1jcH2uayvFjgE/MTS004/NrMKYtKH7r4f+CHwLvAB6T7ZTLz6sEdYfTY9WM5lrQC3kB75DqfGU72Oh83MlgL73f31Pqvi2oenVGiBHwtmNgb4GfDn7n4sc52n/7xHciysmV0PHHT3zVE8/yCVkP63+h/cfQHQTHo6olfEfTgBWEr6D9OZQAWwOIpahiLKPhsMM7sb6AQej7qWHmY2GvhvwP+IupawFFrg5+wc/YNlZqWkw/5xd38qaD5gZtOC9dOAnuv7DlRvrn6OK4AbzGwv8CTpaZ2/A8abWc/FcDKfq7eOYH0lcCSH9UF65NPg7q8Gt9eS/gMQlz5cBOxx90Pu3gE8Rbpf49SHPcLqs/3Bck5qNbNlpC+L+kfBH6bh1HiEgX8Hw3UO6T/srwfvmWrgNTM7Yxj15bQPBy3fc0i5/CI9OtxN+pfUs0Pn/Dw+v5G+0PuP+rTfy8k7z/42WL6Ok3f8/CZon0h6HntC8LUHmBhyrb/HiZ22/8TJO7u+GSz/Z07e4bgmWD6fk3eo7SbcnbYvA7OD5e8F/ReLPgQ+A2wBRgfP+Sjpy3dG3od8eg4/tD7j0zscvxBSjYuBrUBVn+367R9O8R4f6HeQTX191u3lxBx+ZH2Y1Wsk30+Y8x8ovff8bdJ78u/O83N/lvS/zW8Avw2+vkB6fvEFYCfwy4wXgAEPBLW+CdRnPNYtwK7g6+s5qPX3OBH4Zwcvxl3Bm6Y8aE8Ft3cF68/OuP/dQd07CPloA+AiYFPQjz8P3jix6UPgr4DtwFvA/w5CKdI+BFaT3qfQQfq/pOVh9hlQH/y87wD302enehY17iI9593zfnnwdP3DAO/xgX4H2dTXZ/1eTgR+JH2Y7ZdOrSAikhCFNocvIiIDUOCLiCSEAl9EJCEU+CIiCaHAFxFJCAW+SMDSZ+n8ZrB8ppmtjbomkTDpsEyRQHD+o6fdfV7UtYjkQsnpNxFJjHuAc8zst6Q/rHSeu88LPvr/RdLnzZlF+uRpZcAfA22kPzF51MzOIf1hnCqgBVjh7tvz/UOIDERTOiIn3Am84+4XAf+1z7p5wB8AlwB/A7R4+uRurwA3B9usAm5394uB/wL8fT6KFhksjfBFBudX7t4ENJlZI7A+aH8TuDA4Q+rvAP+UcSGj8vyXKTIwBb7I4LRlLHdn3O4m/T4qIn0+9ovyXJfIoGlKR+SEJtKXphwyT1/3YI+ZfQV6r3k6P8ziRLKlwBcJuPsR4NfBRazvHcZD/BGw3MxeJ30K5aVh1ieSLR2WKSKSEBrhi4gkhAJfRCQhFPgiIgmhwBcRSQgFvohIQijwRUQSQoEvIpIQ/x8q6Ro0odz3FAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "39509b1c-f5fb-4edd-847a-97757c2197df",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.6.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}