From 042dc03d76a7d0c64004487edde735c912600e1d Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Wed, 5 May 2021 21:01:08 +0100
Subject: [PATCH] Fixed bugs in the moe&destefano grid
---
binarycpython/utils/distribution_functions.py | 45 ++++++++++-------
binarycpython/utils/grid.py | 49 +++++++++++++++----
setup.py | 2 +-
3 files changed, 67 insertions(+), 29 deletions(-)
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 79431dadd..79e0b3c89 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1162,7 +1162,7 @@ def build_q_table(options, m, p):
# qdata contains the table that we modify: we get
# the original data by interpolating Moe's table
- qdata = []
+ qdata = {}
can_renormalize = 1
qeps = 1e-8 # small number but such that qeps+1 != 1
@@ -1220,10 +1220,15 @@ def build_q_table(options, m, p):
# as a function of M1 and orbital period,
# to obtain the q distribution data.
- for q in np.arange(0.15, np.min(0.950001, qmax + 0.0001), 0.1):
- qdata[q] = Moecache["rinterpolator_q"].interpolate(
+ for q in np.arange(0.15, np.min([0.950001, qmax + 0.0001]), 0.1):
+ val = Moecache["rinterpolator_q"].interpolate(
[np.log10(options[m]), np.log10(options[p]), q]
)[0]
+ # print("val: ", val)
+ # print("q: ", q)
+ # print("qdata: ", qdata)
+ # print("type(qdata): ", type(qdata))
+ qdata[q] = val
# just below qmin, if qmin>qeps, we want nothing
if qmin - 0.15 > qeps:
@@ -1267,10 +1272,10 @@ def build_q_table(options, m, p):
method = options.get("q_{}_extrapolation_method", None)
qlimit = qmin if pre == "log" else qmax
- print("Q: {} method: {}".format(pre, method))
- print("indices: {}".format(indices))
- print("End index: {}".format(end_index))
- print("QS: {}".format(str(indices)))
+ # print("Q: {} method: {}".format(pre, method))
+ # print("indices: {}".format(indices))
+ # print("End index: {}".format(end_index))
+ # print("QS: {}".format(str(indices)))
# truncate the distribution
qdata[max(0.0, min(1.0, qlimit + sign * qeps))] = 0
@@ -1284,12 +1289,12 @@ def build_q_table(options, m, p):
qdata[qlimit] = qdata[qs[end_index]]
elif method == "linear":
# linear extrapolation
- print("Linear 2 {}".format(pre))
+ # print("Linear 2 {}".format(pre))
dq = qs[indices[1]] - qs[indices[0]]
if dq == 0:
# No change
- print("dq = 0")
+ # print("dq = 0")
qdata[qlimit] = qs[end_index]
else:
slope = (
@@ -1299,11 +1304,11 @@ def build_q_table(options, m, p):
qdata[qs[indices[0]]] - slope * qs[indices[0]]
)
qdata[qlimit] = max(0.0, slope * qlimit + intercept)
- print(
- "Slope: {} intercept: {} dn/dq({}) = {}".format(
- slope, intercept, qlimit, qdata[qlimit]
- )
- )
+ # print(
+ # "Slope: {} intercept: {} dn/dq({}) = {}".format(
+ # slope, intercept, qlimit, qdata[qlimit]
+ # )
+ # )
elif method == "plaw2":
newq = 0.05
# use a power-law extrapolation down to q=0.05, if possible
@@ -1394,7 +1399,7 @@ def build_q_table(options, m, p):
)
# TODO: build a check in here to see if the interpolator build was successful
- print("Can renormalize?: {}".format(can_renormalize))
+ # print("Can renormalize?: {}".format(can_renormalize))
if can_renormalize:
# now we integrate and renormalize (if the table is not all zero)
#
@@ -1412,7 +1417,7 @@ def build_q_table(options, m, p):
raise ValueError
else:
I += x[0] * dq
- print("dn/dq ({}) = {}".format())
+ # print("dn/dq ({}) = {} I -> = {}".format(q, x[0], I))
if I > 0:
# normalize to 1.0 by dividing the data by 1.0/$I
@@ -1422,7 +1427,7 @@ def build_q_table(options, m, p):
I = 0
for q in np.arange(0, 1 + 2e-6, dq):
I += q_interpolator.interpolate([q])[0] * dq
- print("Q integral: {}, {}".format(I, q_interpolator))
+ # print("Q integral: {}, {}".format(I, q_interpolator))
# fail if error in integral > 1e-6 (should be ~ machine precision)
if abs(1.0 - I) > 1e-6:
@@ -1430,6 +1435,8 @@ def build_q_table(options, m, p):
# set this new table in the cache
Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)] = q_interpolator
+ if not Moecache.get("rinterpolator_q_metadata", None):
+ Moecache["rinterpolator_q_metadata"] = {}
Moecache["rinterpolator_q_metadata"][m] = options[m]
Moecache["rinterpolator_q_metadata"][p] = options[p]
@@ -1756,7 +1763,7 @@ def Moe_de_Stefano_2017_pdf(options):
options["M1"], len(prob), str(prob), prob_dens
)
)
- elif multipliticty == 2:
+ elif multiplicity == 2:
print(
"M1={} q={} log10P={} ({}): {} -> {}\n".format(
options["M1"],
@@ -1781,7 +1788,7 @@ def Moe_de_Stefano_2017_pdf(options):
prob_dens,
)
)
- elif multipliticty == 4:
+ elif multiplicity == 4:
print(
"M1={} q={} log10P={} ({}) : M3={} P2={}} ecc2={} : M4={} P3={} ecc3={} : {} -> {}".format(
options["M1"],
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 6c2b67cdd..77fbb94c5 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -1130,7 +1130,7 @@ class Population:
verbose_print(
"Process {}: Warning! Ensemble output is empty. ".format(ID),
self.grid_options["verbosity"],
- 2,
+ 1,
)
#
@@ -1145,6 +1145,9 @@ class Population:
ensemble_raw_output
) # Load this into a dict so that we can combine it later
+ #
+ ensemble_json["population_id"] = self.grid_options["_population_id"]
+
else:
# If we do not allow this, automatically we will export this to the data_dir, in some formatted way
output_file = os.path.join(
@@ -1565,6 +1568,7 @@ class Population:
# TODO: Make clear that the phasevol only works good
# TODO: add option to ignore this phasevol calculation and set it to 1
# if you sample linearly in that thing.
+ # if phasevol is <= 0 then we SKIP that whole loop. Its not working then.
if (
not grid_variable["dphasevol"] == -1
): # A method to turn off this calculation and allow a phasevol = 1
@@ -1631,6 +1635,27 @@ class Population:
# Add some whiteline
code_string += indent * (depth + 1) + "\n"
+ ##############3
+ # Add phasevol check:
+ code_string += (
+ indent * (depth + 1)
+ + "if phasevol_{} <= 0:".format(grid_variable["name"])
+ + "\n"
+ )
+
+ # Add phasevol check action:
+ code_string += (
+ indent * (depth + 2)
+ + 'print("phasevol_{} <= 0!")'.format(
+ grid_variable["name"]
+ )
+ + "\n"
+ )
+ code_string += indent * (depth + 2) + "continue" + "\n"
+
+ # Add some whiteline
+ code_string += indent * (depth + 1) + "\n"
+
#########################
# Setting up pre-code and value in some cases
# Add pre-code
@@ -1901,14 +1926,16 @@ class Population:
if not dry_run:
# 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) + "yield(parameter_dict)\n"
- # The below solution might be a good one to add things to specific queues
- # $self->queue_evolution_code_run($self->{_flexigrid}->{thread_q},
- # $system);
-
# 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"
@@ -3227,9 +3254,11 @@ qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
options["ranges"]["q"][0]
if 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)
else "qmax",
+
],
resolution=options["resolutions"]["M"][1],
probdist=1,
@@ -3241,11 +3270,13 @@ 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)
- else "options.get('Mmin', 0.07)/M_1",
+ 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)
+ if options.get("ranges", {}).get("q", [None, None])[1]
else "qmax",
+
options["resolutions"]["M"][1],
),
)
@@ -3470,4 +3501,4 @@ eccentricity3=0
)
# and finally the probability calculator
- # self.grid_options['_grid_variables'][self.last_grid_variable()]['probdist'] = probdist_addition
+ self.grid_options['_grid_variables'][self.last_grid_variable()]['probdist'] = probdist_addition
diff --git a/setup.py b/setup.py
index 502a2067a..966c01581 100644
--- a/setup.py
+++ b/setup.py
@@ -210,7 +210,7 @@ BINARY_C_PYTHON_API_MODULE = Extension(
extra_compile_args=[],
language="C",
)
-
+headers=['src/includes/header.h']
############################################################
# Making the extension function
############################################################
--
GitLab