From 0df4d869c17d4114f9c80be106d2e85355f7e170 Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Thu, 15 Jul 2021 01:20:08 +0100
Subject: [PATCH] updated the grid by splitting the ms functinality and making
it available for other functions
---
binarycpython/utils/grid.py | 789 +++++++++++--------
binarycpython/utils/grid_options_defaults.py | 4 +
2 files changed, 465 insertions(+), 328 deletions(-)
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index bfb90e887..b4d138ff8 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -86,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
@@ -122,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 = {}
@@ -1054,7 +1065,7 @@ class Population:
include_binary_c_version_info=True,
include_binary_c_help_all=True,
)
- self.grid_ensemble_results["metadata"]["settings"] = copy.deepcopy(all_info)
+ self.grid_ensemble_results["metadata"]["settings"] = json.loads(json.dumps(all_info, default=binaryc_json_serializer))
##############################
# Update grid options
@@ -2099,7 +2110,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
@@ -2215,7 +2226,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(
@@ -2223,22 +2242,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")
@@ -2265,16 +2316,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"
@@ -3137,11 +3183,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.
@@ -3216,77 +3257,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):
"""
- Class method of the population class
-
- Takes a dictionary as its only argument
-
- The dictionary should contain a file to
+ 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))
+ # 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))
- 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 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 Moecache:
+ # Load the json data
+ json_data = get_moe_distefano_dataset(self.grid_options['m&s_options'])
- 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"],
@@ -3294,15 +3303,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))
+
+ # 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])
- 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"
+ 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"],
@@ -3310,69 +3326,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"
]
@@ -3383,132 +3430,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"],
@@ -3517,9 +3550,6 @@ class Population:
######
# Setting up the grid variables
- # Make a normal dict of the options again
- options = dict(options)
-
# Multiplicity
self.add_grid_variable(
name="multiplicity",
@@ -3528,11 +3558,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,
@@ -3552,19 +3582,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',
probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
- str(dict(options))[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
+ str(dict(self.grid_options['m&s_options']))[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
),
)
@@ -3575,23 +3605,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
@@ -3601,14 +3631,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",
@@ -3617,35 +3647,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],
),
)
@@ -3656,45 +3686,45 @@ 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]))
+ print("M2 resolution {}".format(self.grid_options['m&s_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",
@@ -3704,35 +3734,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],
),
)
@@ -3742,44 +3772,44 @@ 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]))
+ print("M3 resolution {}".format(self.grid_options['m&s_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)
+ 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)/(m3)",
- 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 "q3max",
],
- resolution=options["resolutions"]["M"][3],
+ resolution=self.grid_options['m&s_options']["resolutions"]["M"][3],
probdist=1,
gridtype="centred",
dphasevol="dq3",
@@ -3789,40 +3819,40 @@ 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)
+ 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)/(m3)",
- 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 "q3max",
- options["resolutions"]["M"][2],
+ self.grid_options['m&s_options']["resolutions"]["M"][2],
),
)
# TODO: Ask about which periods should be used for the calc_sep_from_period
# (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],
),
)
@@ -3832,8 +3862,8 @@ 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(dict(options))[1:-1],
+ updated_options = "{{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}".format(
+ str(dict(self.grid_options['m&s_options']))[1:-1],
'"multiplicity": multiplicity',
'"M1": M_1',
'"M2": M_2',
@@ -3847,11 +3877,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
@@ -3863,6 +3902,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
@@ -3880,3 +3938,78 @@ 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
+ # TODO: add verbose logging
+ print("Using Arenou 2010 to calculate multiplicity fractions")
+ 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
+ # TODO: add verbose logging
+ print("Using Rhagavan 2010 to calculate multiplicity fractions")
+ 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:
+ msg = "Moecache is empty. It needs to be filled with the data for the interpolators"
+ raise ValueError(msg)
+
+ # 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
+ ding = Moe_de_Stefano_2017_multiplicity_fractions(self.grid_options['m&s_options'], self.grid_options["verbosity"])
+ print(ding)
+
+ # Set the prev value back
+ self.grid_options['m&s_options'][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)
+
+ # TODO: add extra verbose logging
+ print("Multiplicity: {} multiplicity_fraction: {}".format(system_dict['multiplicity'], multiplicity_fraction_dict[system_dict['multiplicity']]))
+
+ return multiplicity_fraction_dict[system_dict['multiplicity']]
\ No newline at end of file
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index 707e0714d..52e8c7cf6 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -46,6 +46,8 @@ grid_options_defaults_dict = {
"_zero_prob_stars_skipped": 0,
"ensemble_factor_in_probability_weighted_mass": False, # Whether to multiply the ensemble results by 1/probability_weighted_mass
"multiplicity_fraction_function": 0, # Which multiplicity fraction function to use. 0: None, 1: Arenou 2010, 2: Rhagavan 2010, 3: M&S 2017
+ "_m&s_options": None, # Holds the M&S options.
+ "_loaded_ms_data": False, # Holds flag whether the M&S data is loaded into memory
##########################
# Execution log:
##########################
@@ -489,6 +491,8 @@ grid_options_descriptions = {
"_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",
}
###
--
GitLab