diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 92e2ea41a128a46860452d37508f225e9fd5f62f..07ddd447de6868239badd50da8a3202fc408244d 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 41a0516557707bbf48820b5077465b31fc029b25..dbf90f1139e80be0a615ea21c6d9a199635f5b41 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 dae0c30d3882f58885ef16988f846dcabfd23808..d0edf56b43f0f0df5a87f5c3997c36a895711ae6 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 3a1ff066d763ad352190a615d9e210534e8eebc2..5467f69abc61128b187f4dc3cccbb750eff78d5c 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