From 0b4202572759be1eb811eac0cad825b84ca10052 Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Sun, 6 Jun 2021 18:35:20 +0100
Subject: [PATCH] changing the structure of the pdf output: it will be a
dictionary instead of a float
---
binarycpython/utils/distribution_functions.py | 130 ++++++++++++++++--
binarycpython/utils/grid.py | 31 ++++-
binarycpython/utils/grid_options_defaults.py | 4 +-
3 files changed, 153 insertions(+), 12 deletions(-)
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 68952f5a5..673157013 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1198,6 +1198,15 @@ def build_q_table(options, m, p, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ else:
+ verbose_print(
+ "\tM&S: build_q_table: Cached values for different m={} p={}. Freeing current table and making new table".format(
+ options[m], options[p]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)].destroy()
#
if not incache:
@@ -1618,8 +1627,6 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
mmax => maximum allowed stellar mass (default 80.0)
"""
- verbosity_shift = -2
-
verbose_print(
"\tM&S: Moe_de_Stefano_2017_pdf with options:\n\t\t{}".format(
json.dumps(options)
@@ -1629,6 +1636,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
)
prob = [] # Value that we will return
+ prob_dict = {} # Dictionary containing all the pdf values for the different parameters
# Get the multiplicity from the options, and if its not there, calculate it based on the
# TODO: the function below makes no sense. We NEED to pass the multiplicity in the
@@ -1667,6 +1675,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ prob_dict['multiplicity'] = prob[-1]
############################################################
# always require an IMF for the primary star
@@ -1684,11 +1693,113 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ prob_dict['M1'] = prob[-1]
"""
From here we go through the multiplicities.
"""
if multiplicity >= 2:
+ # If the multiplicity is higher than 1, we will need to construct the following tables:
+ # - period distribution table
+ # - q distribution table
+ # - eccentricity distribution table
+
+ # Set up the interpolator for the periods
+ if not Moecache.get("rinterpolator_log10P", None):
+ Moecache["rinterpolator_log10P"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["period_distributions"], # Contains the table of data
+ nparams=2, # log10M, log10P
+ ndata=2, # binary, triple
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new period interpolator: {}".format(
+ Moecache["rinterpolator_q"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ # Make a table storing Moe's data for q distributions
+ if not Moecache.get("rinterpolator_q", None):
+ Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["q_distributions"], # Contains the table of data
+ nparams=3, # log10M, log10P, q
+ ndata=1, #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
+ Moecache["rinterpolator_q"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ # Make a table storing Moe's data for q distributions
+ if not Moecache.get("rinterpolator_e", None):
+ Moecache["rinterpolator_e"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["ecc_distributions"], # Contains the table of data
+ nparams=3, # log10M, log10P, e
+ ndata=1, #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new e interpolator: {}".format(
+ Moecache["rinterpolator_e"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ if not Moecache.get("rinterpolator_log10P2", None):
+ Moecache["rinterpolator_log10P2"] = py_rinterpolate.Rinterpolate(
+ table=Moecache[
+ "period_distributions"
+ ], # Contains the table of data
+ nparams=2, # log10(M1+M2), log10P # TODO: Really? Is the first column of that table M1+M2?
+ ndata=2, # binary, triple
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+
+ if not Moecache.get("rinterpolator_q", None):
+ Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["q_distributions"], # Contains the table of data
+ nparams=3, # log10(M1+M2), log10P2, q
+ ndata=1, #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
+ Moecache["rinterpolator_q"]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ if not Moecache.get("rinterpolator_log10P2", None):
+ Moecache[
+ "rinterpolator_log10P2"
+ ] = py_rinterpolate.Rinterpolate(
+ table=Moecache[
+ "period_distributions"
+ ], # Contains the table of data
+ nparams=2, # log10(M1+M2), log10P3
+ ndata=2, # binary, triple
+ verbosity=verbosity
+ - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+
+ if not Moecache.get("rinterpolator_q", None):
+ Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
+ table=Moecache[
+ "q_distributions"
+ ], # Contains the table of data
+ nparams=3, # log10(M1+M2), log10P2, q
+ ndata=1, #
+ verbosity=verbosity
+ - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
###############
# Calculation for period of the binary
@@ -1733,6 +1844,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
verbosity,
_MS_VERBOSITY_LEVEL,
)
+ prob_dict['P'] = prob[-1]
############################################################
# mass ratio (0 < q = M2/M1 < qmax)
@@ -1765,7 +1877,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# Construct the q table
build_q_table(options, m_label, p_label, verbosity=verbosity)
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ".format(),
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ({}) for m={} p={}".format(Moecache[
+ "rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
+ ], options[m_label], options[p_label]),
verbosity,
_MS_VERBOSITY_LEVEL,
)
@@ -1794,15 +1908,15 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
# Make a table storing Moe's data for q distributions
if not Moecache.get("rinterpolator_e", None):
- Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
- table=Moecache["q_distributions"], # Contains the table of data
- nparams=3, # log10M, log10P, q
+ Moecache["rinterpolator_e"] = py_rinterpolate.Rinterpolate(
+ table=Moecache["ecc_distributions"], # Contains the table of data
+ nparams=3, # log10M, log10P, e
ndata=1, #
verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
)
verbose_print(
- "\tM&S: Moe_de_Stefano_2017_pdf: Created new q interpolator: {}".format(
- Moecache["rinterpolator_q"]
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created new e interpolator: {}".format(
+ Moecache["rinterpolator_e"]
),
verbosity,
_MS_VERBOSITY_LEVEL,
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 7bd10cc05..26c52d937 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -81,6 +81,10 @@ from binarycpython.utils.distribution_functions import (
fill_data,
)
+from binarycpython.utils.spacing_functions import (
+ const,
+)
+
from binarycpython import _binary_c_bindings
@@ -3481,12 +3485,12 @@ class Population:
),
precode="""orbital_period = 10.0**log10per
qmin={}/M_1
-qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
+qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
""".format(
options.get("Mmin", 0.07)
),
- )
-
+ ) # TODO: change the maximum_mass_ratio_for_RLOF
+
# binaries: mass ratio
self.add_grid_variable(
name="q",
@@ -3752,6 +3756,9 @@ eccentricity3=0
_MS_VERBOSITY_LEVEL,
)
+ # Signal that the M&S grid has been set
+ self.grid_options['_set_ms_grid'] = True
+
def _clean_interpolators(self):
"""
Function to clean up the interpolators after a run
@@ -3768,3 +3775,21 @@ eccentricity3=0
for key in interpolator_keys:
Moecache[key].destroy()
del Moecache[key]
+
+ def plot_MS_grid(self, options, mass):
+ """
+ Function to plot the period distribution along with its probability for a given mass.
+ """
+
+ import matplotlib.pyplot as plt
+
+ plt.title("Period distribution for M={}".format(mass))
+ plt.plot()
+ plt.show()
+
+ sampling_range = const(
+ options["ranges"]["logP"][0],
+ options["ranges"]["logP"][1],
+ options["resolutions"]["logP"][0])
+
+
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index 971a13ecd..ae1f87736 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -40,10 +40,11 @@ grid_options_defaults_dict = {
"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
##########################
# Execution log:
##########################
- "verbosity": 0, # Level of verbosity of the simulation. 0=INFO,
+ "verbosity": 0, # Level of verbosity of the simulation
"log_file": os.path.join(
temp_dir(), "binary_c_python.log"
), # Set to None to not log to file. The directory will be created
@@ -478,6 +479,7 @@ grid_options_descriptions = {
"_total_probability_weighted_mass_run": "To count the total mass * probability for each system that thread/process has ran",
"_actually_evolve_system": "Whether to actually evolve the systems of just act as if. for testing. used in _process_run_population_grid",
"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",
}
###
--
GitLab