From 65f804c8903e2c8f77552a900206c5bb0ef70c23 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Wed, 13 Oct 2021 14:09:31 +0100
Subject: [PATCH] added options to the Moe loader so we can change the JSON
data and reload it after clearing the cache
the period distribution renormalization now gets the low and high limits from the loaded data
update_dicts() now copies the type of dict_1 for the new dict, rather than using OrderedDict (maybe merge_dicts should do this too?)
---
binarycpython/utils/distribution_functions.py | 36 ++++--
binarycpython/utils/functions.py | 46 +++++---
binarycpython/utils/grid.py | 106 ++++++++++++------
binarycpython/utils/grid_options_defaults.py | 5 +
4 files changed, 130 insertions(+), 63 deletions(-)
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 92e2ea41a..07ddd447d 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1768,7 +1768,7 @@ def calc_e_integral(
def calc_P_integral(
- options, min_P, integrals_string, interpolator_name, mass_string, verbosity=0
+ options, min_logP, max_logP, integrals_string, interpolator_name, mass_string, verbosity=0
):
"""
Function to calculate the P integral
@@ -1777,7 +1777,6 @@ def calc_P_integral(
"""
global Moecache
- max_logP = 10
# Check if the dict exists
if not Moecache.get(integrals_string, None):
@@ -1788,7 +1787,7 @@ def calc_P_integral(
I = 0
dlogP = 1e-3
- for logP in np.arange(np.log10(min_P), max_logP, dlogP):
+ for logP in np.arange(min_logP, max_logP, dlogP):
# Loop over all the values in the table, between the min and max P
dp_dlogP = Moecache[interpolator_name].interpolate(
[np.log10(options[mass_string]), logP]
@@ -1799,8 +1798,8 @@ def calc_P_integral(
# Set the integral value in the dict
Moecache[integrals_string][options[mass_string]] = I
verbose_print(
- "\tM&S: calc_P_integral: min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
- min_P,
+ "\tM&S: calc_P_integral: min_logP: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
+ min_logP,
integrals_string,
interpolator_name,
mass_string,
@@ -1812,8 +1811,8 @@ def calc_P_integral(
)
else:
verbose_print(
- "\tM&S: calc_P_integral: Found cached value for min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
- min_P,
+ "\tM&S: calc_P_integral: Found cached value for min_logP: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
+ min_logP,
integrals_string,
interpolator_name,
mass_string,
@@ -2000,14 +1999,19 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
options["sep"] = calc_sep_from_period(
options["M_1"], options["M_2"], options["P"]
)
- # TODO: add check for min_P with instant RLOF?
+ # TODO: add check for min_logP with instant RLOF?
# TODO: Actually use the value above.
# Total mass inner binary:
options["M_1+M_2"] = options["M_1"] + options["M_2"]
# Calculate P integral or use cached value
+
+ # get the periods from the Moecahe
+ min_logP = Moecache['logperiods'][0]
+ max_logP = Moecache['logperiods'][-1]
+
calc_P_integral(
- options, 1, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
+ options, min_logP, max_logP, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
)
# Set probabilty for P1
@@ -2125,7 +2129,13 @@ def Moe_di_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["M_1+M_2"], options["mmin"], max_sep)
+ min_P2 = calc_period_from_sep(options["M_1+M_2"],
+ options["mmin"],
+ max_sep)
+ min_logP2 = math.log10(min_P2)
+ #max_logP2 = 10.0
+ #min_logP = Moecache['logperiods'][0]
+ max_logP2 = Moecache['logperiods'][-1]
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
@@ -2152,7 +2162,8 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
# Check for cached value of P integral or calculate
calc_P_integral(
options,
- min_P2,
+ min_logP2,
+ max_logP2,
"P2_integrals",
"rinterpolator_log10P",
"M_1+M_2",
@@ -2243,7 +2254,8 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
# TODO: Make sure we use the correct period idea here.
calc_P_integral(
options,
- min_P2,
+ min_logP2,
+ max_logP2,
"P2_integrals",
"rinterpolator_log10P",
"M_1+M_2",
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 41a051655..dbf90f113 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -410,19 +410,19 @@ def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
def get_moe_di_stefano_dataset(options, verbosity=0):
"""
- Function to get the default moe and di Stefano dataset or accept a user input.
+ Function to get the default Moe and di Stefano dataset or accept a user input.
+
+ Returns a dict containing the (JSON) data.
"""
- if not options.get("file", None):
- verbose_print(
- "Using the default Moe and de Stefano 2017 datafile",
- verbosity,
- 1,
- )
+ json_data = None
- json_data = copy.deepcopy(moe_di_stefano_2017_data.moe_di_stefano_2017_data)
+ if "JSON" in options:
+ # use the JSON data passed in
+ json_data = options["JSON"]
- else:
+ elif "file" in options:
+ # use the file passed in, if provided
if not os.path.isfile(options["file"]):
verbose_print(
"The provided 'file' Moe and de Stefano JSON file does not seem to exist at {}".format(
@@ -444,10 +444,19 @@ def get_moe_di_stefano_dataset(options, verbosity=0):
# Read input data and Clean up the data if there are white spaces around the keys
with open(options["file"], "r") as data_filehandle:
datafile_data = data_filehandle.read()
- datafile_data = datafile_data.replace('" ', '"')
- datafile_data = datafile_data.replace(' "', '"')
- datafile_data = datafile_data.replace(' "', '"')
- json_data = json.loads(datafile_data)
+ datafile_data = datafile_data.replace('" ', '"')
+ datafile_data = datafile_data.replace(' "', '"')
+ datafile_data = datafile_data.replace(' "', '"')
+ json_data = json.loads(datafile_data)
+
+ if not json_data:
+ # no JSON data or filename given, use the default 2017 dataset
+ verbose_print(
+ "Using the default Moe and de Stefano 2017 datafile",
+ verbosity,
+ 1,
+ )
+ json_data = copy.deepcopy(moe_di_stefano_2017_data.moe_di_stefano_2017_data)
return json_data
@@ -1905,10 +1914,15 @@ def update_dicts(dict_1: dict, dict_2: dict) -> dict:
"""
- # Set up new dict
- new_dict = OrderedDict() # TODO: check if this still works
-
+ # Set up new dict of the same type as dict_1
#
+ # Note: setting directly to OrderedDict fails in some cases
+ # so instead we take a (shallow) copy of dict_1 which will
+ # have the same type as dict_1, then clear it. (There must
+ # be a better way to do this...)
+ new_dict = dict_1.copy() # OrderedDict() # TODO: check if this still works
+ new_dict.clear()
+
keys_1 = dict_1.keys()
keys_2 = dict_2.keys()
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index dae0c30d3..d0edf56b4 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -151,8 +151,9 @@ class Population:
# Set the options that are passed at creation of the object
self.set(**kwargs)
- # Load M&s options
+ # Load Moe and di Stefano options
self.grid_options["m&s_options"] = copy.deepcopy(moe_di_stefano_default_options)
+ self.custom_options["m&s_JSON_data"] = None
# Write M&S options to a file. NOTE: not sure why i put this here anymore
os.makedirs(
@@ -3745,30 +3746,33 @@ class Population:
def set_moe_di_stefano_settings(self, options=None):
"""
- Function to set user input configurations for the moe & di Stefano methods
+ Function to set user input configurations for the Moe & di Stefano methods
If nothing is passed then we just use the default 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))
+ if not options:
+ options = {}
- # Write options to a file
- os.makedirs(
- os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
- exist_ok=True,
- )
- with open(
+
+ # 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(options)
+
+ # Write options to a file
+ os.makedirs(
+ os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
+ exist_ok=True,
+ )
+ with open(
os.path.join(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
"moeopts.dat",
),
"w",
- ) as f:
- f.write(json.dumps(self.grid_options["m&s_options"], indent=4))
- f.close()
+ ) as f:
+ f.write(json.dumps(self.grid_options["m&s_options"], indent=4))
+ f.close()
def _load_moe_di_stefano_data(self):
"""
@@ -3776,15 +3780,26 @@ class Population:
"""
# 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_di_stefano_dataset(
- self.grid_options["m&s_options"],
- verbosity=self.grid_options["verbosity"],
- )
+ if (not self.grid_options["_loaded_ms_data"]): # and not Moecache:
- # 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 self.custom_options["m&s_JSON_data"]:
+ # Use the existing (perhaps modified) JSON data
+ json_data = self.custom_options["m&s_JSON_data"]
+
+ else:
+ # Load the JSON data from a file
+ json_data = get_moe_di_stefano_dataset(
+ self.grid_options["m&s_options"],
+ verbosity=self.grid_options["verbosity"],
+ )
+
+ # entry of log10M1 is a list containing 1 dict.
+ # We can take the dict out of the list
+ if isinstance(json_data["log10M1"], list):
+ json_data["log10M1"] = json_data["log10M1"][0]
+
+ # save this data in case we want to modify it later
+ self.custom_options["m&s_JSON_data"] = json_data
# Get all the masses
logmasses = sorted(json_data["log10M1"].keys())
@@ -3815,13 +3830,12 @@ class Population:
logperiods = []
for logmass in logmasses:
if not logperiods:
- logperiods = sorted(json_data["log10M1"][logmass]["logP"])
+ logperiods = sorted(json_data["log10M1"][logmass]["logP"].keys())
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"
+ msg = "Period values are not consistent throughout the dataset\logperiods = " + ' '.join(str(x) for x in logperiods) + "\nCurrent periods = " + ' '.join(str(x) for x in current_logperiods)
verbose_print(
"\tMoe_di_Stefano_2017: {}".format(msg),
self.grid_options["verbosity"],
@@ -3832,7 +3846,6 @@ class Population:
############################################################
# log10period binwidth : of course this assumes a fixed
# binwidth, so we check for this too.
-
for i in range(len(current_logperiods) - 1):
if not dlog10P == (
float(current_logperiods[i + 1]) - float(current_logperiods[i])
@@ -3845,6 +3858,10 @@ class Population:
)
raise ValueError(msg)
+ # save the logperiods list in the cache:
+ # this is used in the renormalization integration
+ Moecache['logperiods'] = logperiods
+
# Write to file
os.makedirs(
os.path.join(self.grid_options["tmp_dir"], "moe_distefano"),
@@ -3884,13 +3901,12 @@ class Population:
# 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
+ # divide by dlog10P
if first:
first = 0
@@ -3923,6 +3939,7 @@ class Population:
/ dlog10P,
]
)
+
Moecache["period_distributions"].append(
[
float(logmass),
@@ -4095,9 +4112,6 @@ class Population:
############################################################
# always require M1, for all systems
#
-
- # TODO: put in option for the time-adaptive grid.
-
# log-spaced m1 with given resolution
self.add_grid_variable(
name="lnm1",
@@ -4478,14 +4492,36 @@ eccentricity3=0
Takes a dictionary as its only argument
"""
+ default_options = {
+ 'apply settings' : True,
+ 'setup grid' : True,
+ 'load data' : True,
+ 'clean cache' : False,
+ 'clean load flag' : False,
+ 'clean all' : False,
+ }
+ if not options:
+ options = {}
+ options = update_dicts(default_options,options)
+
+ # clean cache?
+ if options['clean all'] or options['clean cache']:
+ Moecache.clear()
+
+ if options['clean all'] or options['clean load flag']:
+ self.grid_options["_loaded_ms_data"] = False
+
# Set the user input
- self.set_moe_di_stefano_settings(options=options)
+ if options['apply settings']:
+ self.set_moe_di_stefano_settings(options=options)
# Load the data
- self._load_moe_di_stefano_data()
+ if options['load data']:
+ self._load_moe_di_stefano_data()
# construct the grid here
- self._set_moe_di_stefano_distributions()
+ if options['setup grid']:
+ self._set_moe_di_stefano_distributions()
def _clean_interpolators(self):
"""
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index 3a1ff066d..5467f69ab 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -511,6 +511,11 @@ MIN_MASS_BINARY_C = float(
# Default options for the Moe & di Stefano grid
moe_di_stefano_default_options = {
+ # place holder for the JSON data to be used if a file
+ # isn't specified
+ "JSON" : None,
+
+ # resolution data
"resolutions": {
"M": [
20, # M1
--
GitLab