diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 9c2e39510db99b7ca2fd0fdbe9284d5283232375..68952f5a5fa3aad6c07802dd9abcc77a217a5da2 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -21,9 +21,18 @@ Tasks:
"""
import math
+import json
+
import numpy as np
+
from typing import Optional, Union
+
from binarycpython.utils.useful_funcs import calc_period_from_sep, calc_sep_from_period
+from binarycpython.utils.functions import verbose_print
+from binarycpython.utils.grid_options_defaults import (
+ _MS_VERBOSITY_LEVEL,
+ _MS_VERBOSITY_INTERPOLATOR_LEVEL,
+)
###
# File containing probability distributions
@@ -943,9 +952,10 @@ import py_rinterpolate
Moecache = {}
-def poisson(lambda_val, n, nmax=None):
+def poisson(lambda_val, n, nmax=None, verbosity=0):
"""
- Function that calculates the poisson value and normalizes TODO: improve the description
+ Function that calculates the poisson value and normalizes
+ TODO: improve the description
"""
cachekey = "{} {} {}".format(lambda_val, n, nmax)
@@ -953,11 +963,15 @@ def poisson(lambda_val, n, nmax=None):
if distribution_constants.get("poisson_cache", None):
if distribution_constants["poisson_cache"].get(cachekey, None):
p_val = distribution_constants["poisson_cache"][cachekey]
- # print(
- # "Found cached value: Poisson ({}, {}, {}) = {}\n".format(
- # lambda_val, n, nmax, p_val
- # )
- # )
+
+ verbose_print(
+ "\tM&S: found cached value for poisson({}, {}, {}): {}".format(
+ lambda_val, n, nmax, p_val
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
return p_val
# Poisson distribution : note, n can be zero
@@ -977,7 +991,11 @@ def poisson(lambda_val, n, nmax=None):
distribution_constants["poisson_cache"] = {}
distribution_constants["poisson_cache"][cachekey] = p_val
- # print("Poisson ({}, {}, {}) = {}\n".format(lambda_val, n, nmax, p_val))
+ verbose_print(
+ "\tM&S: Poisson({}, {}, {}): {}".format(lambda_val, n, nmax, p_val),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
return p_val
@@ -989,11 +1007,15 @@ def _poisson(lambda_val, n):
return (lambda_val ** n) * np.exp(-lambda_val) / (1.0 * math.factorial(n))
-def Moe_de_Stefano_2017_multiplicity_fractions(options):
+def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
+ """
+ Function that returns a list of multiplicity fractions for a given list
+ TODO: make an extrapolation functionality in this. log10(1.6e1) is low, we can probably go a bit further
+ """
# Returns a list of multiplicity fractions for a given input list of masses
- global Moecache
- # TODO: make an extrapolation functionality in this. log10(1.6e1) is low, we can probably go a bit further
+ # Use the global Moecache
+ global Moecache
result = {}
@@ -1020,6 +1042,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options):
table=Moecache["multiplicity_table"], # Contains the table of data
nparams=1, # logM1
ndata=4, # The amount of datapoints (the parameters that we want to interpolate)
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
)
if options["multiplicity_model"] == "Poisson":
@@ -1029,7 +1052,8 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options):
for n in range(4):
result[n] = (
- options["multiplicity_modulator"][n] * poisson(multiplicity, n, 3)
+ options["multiplicity_modulator"][n]
+ * poisson(multiplicity, n, 3, verbosity)
if options["multiplicity_modulator"][n] > 0
else 0
)
@@ -1053,11 +1077,13 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options):
if options["normalize_multiplicities"] == "raw":
# do nothing : use raw Poisson predictions
pass
+
elif options["normalize_multiplicities"] == "norm":
# simply normalize so all multiplicities add to 1
sum_result = sum([result[key] for key in result.keys()])
for key in result.keys():
result[key] = result[key] / sum_result
+
elif options["normalize_multiplicities"] == "merge":
# if multiplicity == 1, merge binaries, triples and quads into singles
# (de facto this is the same as "norm")
@@ -1067,11 +1093,16 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options):
# TODO: ask rob about this part. Not sure if i understand it.
- multiplicitty = 0
+ multiplicity = 0
for n in range(4):
if options["multiplicity_modulator"][n] > 0:
multiplicity = n + 1
- print("Multiplicity: {}: {}".format(multiplicity, str(result)))
+
+ verbose_print(
+ "\tM&S: Multiplicity: {}: {}".format(multiplicity, str(result)),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
if multiplicity == 1:
result[0] = 1.0
@@ -1098,24 +1129,31 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options):
# so do nothing
pass
else:
- print(
- "Error: in the Moe distribution, we seem to want no stars at all (multiplicity == {}).\n".format(
- multiplicity
- )
+ msg = "\tM&S: Error: in the Moe distribution, we seem to want no stars at all (multiplicity == {})".format(
+ multiplicity
+ )
+ verbose_print(
+ msg,
+ verbosity,
+ 0,
)
+ raise ValueError(msg)
sum_result = sum([result[key] for key in result.keys()])
for key in result.keys():
result[key] = result[key] / sum_result
- # print("Multiplicity array: {}".format(str(result)))
+ verbose_print(
+ "\tM&S: Multiplicity array: {}".format(str(result)),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# return array reference
return result
-# @profile
-def build_q_table(options, m, p):
+def build_q_table(options, m, p, verbosity=0):
############################################################
#
# Build an interpolation table for q, given a mass and
@@ -1132,24 +1170,34 @@ def build_q_table(options, m, p):
# mass is in $opts->{$m}
# period is $opts->{$p}
#
- ############################################################
-
# Since the information from the table for M&S is independent of any choice we make,
- # we need to take into account that for example our choice of minimum mass leads to a minimum q_min that is not the same as in the table
- # We should ignore those parts of the table and renormalize. If we are below the lowest value of qmin in the table we need to extrapolate the data
+ # we need to take into account that for example our choice of minimum mass leads to
+ # a minimum q_min that is not the same as in the table
+ # We should ignore those parts of the table and renormalize.
+ # If we are below the lowest value of qmin in the table we need to extrapolate the data
+ #
+ # Anyway, the goal of this function is to provide some extrapolated values for q when we should sample outside of the boundaries
+ ############################################################
# We can check if we have a cached value for this already:
# TODO: fix this cache check.
incache = False
if Moecache.get("rinterpolator_q_metadata", None):
- if (Moecache["rinterpolator_q_metadata"][m]) and (
- Moecache["rinterpolator_q_metadata"][p]
+ if (Moecache["rinterpolator_q_metadata"].get(m, None)) and (
+ Moecache["rinterpolator_q_metadata"].get(p, None)
):
if (Moecache["rinterpolator_q_metadata"][m] == options[m]) and (
Moecache["rinterpolator_q_metadata"][p] == options[p]
):
incache = True
- # print("INCACHE: {}".format(incache))
+
+ verbose_print(
+ "\tM&S: build_q_table: Found cached values for m={} p={}".format(
+ options[m], options[p]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
#
if not incache:
@@ -1159,7 +1207,8 @@ def build_q_table(options, m, p):
# the companions are planets
qmin = options["ranges"]["M"][
0
- ] # TODO: this lower range must not be lower than Mmin. However, since the q_min is used as a sample range for the M2, it should use this value I think. discuss
+ ] # TODO: this lower range must not be lower than Mmin.
+ # However, since the q_min is used as a sample range for the M2, it should use this value I think. discuss
# TODO: fix that this works. Should depend on metallicity too I think. iirc this function does not work now.
# qmax = maximum_mass_ratio_for_RLOF(options[m], options[p])
@@ -1230,10 +1279,6 @@ def build_q_table(options, m, p):
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
@@ -1267,261 +1312,361 @@ def build_q_table(options, m, p):
qs[1] = 1.0
qdata[qs[1]] = qs[0]
+ # We actually should do the extrapolation now.
else:
+ # Loop over both the lower end and the upper end
for pre in ["low", "high"]:
if require_extrapolation[pre] == 0:
continue
else:
sign = -1 if pre == "low" else 1
end_index = 0 if pre == "low" else len(qs) - 1
- indices = [0, 1] if pre == "low" else [len(qs) - 1, len(qs) - 2]
- 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)))
+ indices = (
+ [0, 1] if pre == "low" else [len(qs) - 1, len(qs) - 2]
+ ) # Based on whether we do the high or low end we need to use two different indices
+ method = options.get(
+ "q_{}_extrapolation_method".format(pre), None
+ )
+ qlimit = qmin if pre == "low" else qmax
+
+ verbose_print(
+ "\tM&S: build_q_table: Extrapolating: Q: {} method: {}, indices: {} End index: {}".format(
+ pre, method, indices, end_index
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# truncate the distribution
qdata[max(0.0, min(1.0, qlimit + sign * qeps))] = 0
if method == None:
# no extrapolation : just interpolate between 0.10 and 0.95
+ verbose_print(
+ "\tM&S: build_q_table: using no extrapolations".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
continue
elif method == "flat":
# use the end value value and extrapolate it
# with zero slope
qdata[qlimit] = qdata[qs[end_index]]
+ verbose_print(
+ "\tM&S: build_q_table: using constant extrapolation".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
elif method == "linear":
- # linear extrapolation
- # print("Linear 2 {}".format(pre))
- dq = qs[indices[1]] - qs[indices[0]]
-
- if dq == 0:
- # No change
- # print("dq = 0")
- qdata[qlimit] = qs[end_index]
- else:
- slope = (
- qdata[qs[indices[1]]] - qdata[qs[indices[0]]]
- ) / dq
- intercept = (
- 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]
- # )
- # )
+ qdata[qlimit] = linear_extrapolation_q(
+ qs=qs,
+ indices=indices,
+ qlimit=qlimit,
+ qdata=qdata,
+ end_index=end_index,
+ verbosity=verbosity,
+ )
+
+ verbose_print(
+ "\tM&S: build_q_table: using linear extrapolation".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
elif method == "plaw2":
- newq = 0.05
- # use a power-law extrapolation down to q=0.05, if possible
- if (qdata[qs[indices[0]]] == 0) and (
- qdata[qs[indices[1]]] == 0.0
- ):
- # not possible
- qdata[newq] = 0
- else:
- slope = (
- np.log10(qdata[qs[indices[1]]])
- - np.log10(qdata[qs[indices[0]]])
- ) / (
- np.log10(qs[indices[1]]) - np.log10(qs[indices[0]])
- )
- intercept = np.log10(
- qdata[qs[indices[0]]]
- ) - slope * log10(qs[indices[0]])
- qdata[newq] = slope * newq + intercept
-
+ qdata[qlimit] = powerlaw_extrapolation_q(
+ qs=qs, indices=indices, qdata=qdata, verbosity=verbosity
+ )
+
+ verbose_print(
+ "\tM&S: build_q_table: using powerlaw extrapolation".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
elif method == "nolowq":
newq = 0.05
qdata[newq] = 0
+ verbose_print(
+ "\tM&S: build_q_table: setting lowq to 0".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
else:
- print("No other methods available")
+ verbose_print(
+ "\tM&S: build_q_table: Error no other methods available. The chosen method ({}) does not exist!".format(
+ method
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
raise ValueError
- # TODO: consider implementing this
- # elsif($method =~ /^(log)?poly(\d+)/)
- # {
- # ############################################################
- # # NOT WORKING / TESTED
- # ############################################################
-
- # # fit a polynomial of degree n in q or log10(q)
- # my $dolog = defined $1 ? 1 : 0;
- # my $n = $2;
-
- # # make xdata : list of qs
- # my @xdata;
- # if($dolog)
- # {
- # @xdata = @qs; # linear
- # }
- # else
- # {
- # @xdata = map{log10(MAX(1e-20,$_))}@qs; # log
- # }
-
- # # make ydata : from JSON
- # my @ydata = map{$qdata->{$_}}@qs; # y data
-
- # # make parameters: these are 1/nn
- # my @parameters; # parameters
- # $#parameters = $n - 1; # set length of parameters array
- # {
- # my $nn = $n;
- # @parameters = map{1.0/$nn--}@parameters; # set all to 1.0/parameter number
- # }
-
- # # make polynomial with parameters
- # my $formula = Math::Polynomial->new(@parameters); # use a polynomial fit
- # my $variable ='x';
- # my $max_iter = 100;
- # my $square_residual = Algorithm::CurveFit->curve_fit(
- # $formula,
- # \@parameters,
- # $variable,
- # \@xdata,
- # \@ydata,
- # $max_iter);
-
- # # evaluate at q=0.05
- # my $newq = 0.05;
- # $qdata->{$newq} = $formula->evaluate($newq);
-
- # print Data::Dumper::Dumper(\@parameters);
- # exit;
- # }
-
- # regenerate qs in new table
+ # TODO: consider implementing the poly method (see perl version)
+
+ # regenerate qs in new table. This is now the updated list of qs where we have some extrapolated numbers
tmp_table = []
for q in sorted(qdata.keys()):
tmp_table.append([q, qdata[q]])
- if Moecache.get("rinterpolator_q_given_{}_log10{}".format(m, p), None):
- print(
- "Present interpolator: {}".format(
- Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)]
- )
- )
- print("Destroying present interpolator:")
- interpolator = Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)]
- print(interpolator)
- print(type(interpolator))
- print(dir(interpolator))
- x = Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)].interpolate(
- [0.5]
- )
- print("Interpolated a value q=0.5: {}".format(x))
- Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)].destroy()
- print(interpolator)
- print(type(interpolator))
- print(dir(interpolator))
- print(
- "Present interpolator: {}".format(
- Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)]
- )
- )
- x = Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)].interpolate(
- [0.5]
- )
- print("Interpolated a value q=0.5: {}".format(x))
- # del Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)]
-
- print("CREATING A NEW TABLE Q table")
# Make an interpolation table to contain our modified data
q_interpolator = py_rinterpolate.Rinterpolate(
- table=tmp_table, nparams=1, ndata=1 # Contains the table of data # q #
+ table=tmp_table,
+ nparams=1,
+ ndata=1, # Contains the table of data # q #
+ verbosity=verbosity - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
+ )
+ verbose_print(
+ "\tM&S: build_q_table: Created a new Q table",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
)
- print("CREATed A NEW TABLE Q table")
- # TODO: build a check in here to see if the interpolator build was successful
- # print("Can renormalize?: {}".format(can_renormalize))
if can_renormalize:
+ verbose_print(
+ "\tM&S: build_q_table: Renormalizing table",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
# now we integrate and renormalize (if the table is not all zero)
- #
- dq = 1e-3 # resolution of the integration/renormalization
- I = 0
-
- # integrate: note that the value of the integral is
- # meaningless to within a factor (which depends on $dq)
- for q in np.arange(0, 1 + 2e-6, dq):
- x = q_interpolator.interpolate([q])
- if len(x) == 0:
- print("Q interpolator table interpolation failed")
- print("tmp_table = {}".format(str(tmp_table)))
- print("q_data = {}".format(str(qdata)))
- raise ValueError
- else:
- I += x[0] * dq
- # print("dn/dq ({}) = {} I -> = {}".format(q, x[0], I))
+ I = get_integration_constant_q(q_interpolator, verbosity=verbosity)
if I > 0:
# normalize to 1.0 by dividing the data by 1.0/$I
q_interpolator.multiply_table_column(1, 1.0 / I)
# test this
- 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))
+ new_I = get_integration_constant_q(q_interpolator, verbosity=verbosity)
# fail if error in integral > 1e-6 (should be ~ machine precision)
- if abs(1.0 - I) > 1e-6:
- print("Error: > 1e-6 in q probability integral: {}".format(I))
-
+ if abs(1.0 - new_I) > 1e-6:
+ verbose_print(
+ "\tM&S: build_q_table: Error: > 1e-6 in q probability integral: {}".format(
+ I
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# set this new table in the cache
- print(
- "STORING Q INTERPOLATOR AS {}".format(
- "rinterpolator_q_given_{}_log10{}".format(m, p)
- )
- )
Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)] = q_interpolator
- print(
- "STORed Q INTERPOLATOR AS {}".format(
+ verbose_print(
+ "\tM&S: build_q_table: stored q_interpolater as {}".format(
"rinterpolator_q_given_{}_log10{}".format(m, p)
- )
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
)
+ # Store the values for which this table was set up in the dict
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]
-def Moe_de_Stefano_2017_pdf(options):
- # Moe and de Stefano probability density function
- #
- # takes a dctionary as input (in options) with options:
- #
- # M1, M2, M3, M4 => masses (Msun) [M1 required, rest optional]
- # P, P2, P3 => periods (days) [number: none=binary, 2=triple, 3=quadruple]
- # ecc, ecc2, ecc3 => eccentricities [numbering as for P above]
- #
- # mmin => minimum allowed stellar mass (default 0.07)
- # mmax => maximum allowed stellar mass (default 80.0)
- #
+def powerlaw_extrapolation_q(qdata, qs, indices, verbosity=0):
+ """
+ Function to do the powerlaw extrapolation at the lower end of the q range
+ """
+ newq = 0.05
+
+ # use a power-law extrapolation down to q=0.05, if possible
+ if (qdata[qs[indices[0]]] == 0.0) and (qdata[qs[indices[1]]] == 0.0):
+ # not possible
+ return 0
+
+ else:
+ slope = (np.log10(qdata[qs[indices[1]]]) - np.log10(qdata[qs[indices[0]]])) / (
+ np.log10(qs[indices[1]]) - np.log10(qs[indices[0]])
+ )
+ intercept = np.log10(qdata[qs[indices[0]]]) - slope * np.log10(qs[indices[0]])
+
+ return slope * newq + intercept
+
+
+def linear_extrapolation_q(qs, indices, qlimit, qdata, end_index, verbosity=0):
+ """
+ Function to do the linear extrapolation for q.
+ """
+
+ # linear extrapolation
+ dq = qs[indices[1]] - qs[indices[0]]
+
+ if dq == 0:
+ verbose_print(
+ "\tM&S: build_q_table: linear dq=0".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ # No change
+ return qs[end_index]
+ else:
+ slope = (qdata[qs[indices[1]]] - qdata[qs[indices[0]]]) / dq
+
+ intercept = qdata[qs[indices[0]]] - slope * qs[indices[0]]
+ qdata[qlimit] = max(0.0, slope * qlimit + intercept)
+ verbose_print(
+ "\tM&S: build_q_table: linear Slope: {} intercept: {} dn/dq({}) = {}".format(
+ slope, intercept, qlimit, qdata[qlimit]
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ return max(0.0, slope * qlimit + intercept)
+
+
+def get_integration_constant_q(q_interpolator, verbosity=0):
+ """
+ Function to integrate the q interpolator and return the integration constant
+ """
+
+ dq = 1e-3 # resolution of the integration/renormalization
+ I = 0
+
+ # integrate: note that the value of the integral is
+ # meaningless to within a factor (which depends on $dq)
+ for q in np.arange(0, 1 + 2e-6, dq):
+ x = q_interpolator.interpolate([q])
+ if len(x) == 0:
+ verbose_print(
+ "\tM&S: build_q_table: Q interpolator table interpolation failed.\n\t\ttmp_table = {}\n\t\tq_data = {}".format(
+ str(tmp_table), str(qdata)
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ raise ValueError
+ else:
+ I += x[0] * dq
+ # verbose_print(
+ # "\tM&S: build_q_table: dn/dq ({}) = {} I -> = {}".format(q, x[0], I),
+ # verbosity,
+ # _MS_VERBOSITY_LEVEL,
+ # )
+ return I
+
+
+def fill_data(sample_values, data_dict):
+ """
+ Function that returns the normalized array of values for given logmass and logperiod
+ used for the e and q values
+
+ TODO: make sure we do the correct thing with the dstep
+ """
+
+ data = {}
+ I = 0
+
+ dstep = float(sample_values[1]) - float(sample_values[0])
+
+ # Read out the data
+ for sample_value in sample_values:
+ val = data_dict[sample_value]
+ data[sample_value] = val
+ I += val
+
+ # Normalize the data
+ for sample_value in sample_values:
+ data[sample_value] = data[sample_value] / I
+
+ return data
+
+
+def calc_P_integral(
+ options, min_P, integrals_string, interpolator_name, mass_string, verbosity=0
+):
+ """
+ Function to calculate the P integral
+
+ We need to renormalize this because min_per > 0, and not all periods should be included
+ """
+
+ global Moecache
+ max_logP = 10
+
+ I = 0
+ dlogP = 1e-3
+
+ for logP in np.arange(np.log10(min_P), 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]
+ )[0]
+
+ I += dp_dlogP * dlogP
+ Moecache[integrals_string][options[mass_string]] = I
+ verbose_print(
+ "\tM&S: calc_P_integral: min_P: {} integrals_string: {} interpolator_name: {} mass_string: {} I: {}".format(
+ min_P, integrals_string, interpolator_name, mass_string, I
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+
+def Moe_de_Stefano_2017_pdf(options, verbosity=0):
+ """
+ Moe & distefano function to calculate the probability density.
+
+
+ takes a dctionary as input (in options) with options:
+
+ M1, M2, M3, M4 => masses (Msun) [M1 required, rest optional]
+ P, P2, P3 => periods (days) [number: none=binary, 2=triple, 3=quadruple]
+ ecc, ecc2, ecc3 => eccentricities [numbering as for P above]
+
+ mmin => minimum allowed stellar mass (default 0.07)
+ 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)
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
prob = [] # Value that we will return
- multiplicity = options["multiplicity"]
+ # 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
if not options.get("multiplicity", None):
- multiplicity = 1
- for n in range(2, 5):
- multiplicity += 1 if options.get("M{}".format(n), None) else 0
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Did not find a multiplicity value in the options dictionary",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ raise ValueError
+ # multiplicity = 1
+ # for n in range(2, 5):
+ # multiplicity += 1 if options.get("M{}".format(n), None) else 0
else:
multiplicity = options["multiplicity"]
- # immediately return 0 if the multiplicity modulator is 0
+ # Immediately return 0 if the multiplicity modulator is 0
if options["multiplicity_modulator"][multiplicity - 1] == 0:
- print("_pdf ret 0 because of mult mod\n")
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: returning 0 because of the multiplicity modulator being 0",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
return 0
############################################################
# multiplicity fraction
- prob.append(Moe_de_Stefano_2017_multiplicity_fractions(options)[multiplicity - 1])
+ # Calculate the probabilty, or rather, fraction, of stars that belong to this mass
+ prob.append(
+ Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity)[multiplicity - 1]
+ )
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob list ({})".format(
+ options["M1"], prob[-1], prob
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
############################################################
# always require an IMF for the primary star
@@ -1530,59 +1675,63 @@ def Moe_de_Stefano_2017_pdf(options):
# (dlnM = dM/M, so 1/dlnM = M/dM)
# TODO: Create an n-part-powerlaw method that can have breakpoints and slopes. I'm using a three-part powerlaw now.
+ # TODO: is this actually the correct way? putting the M1 in there? Do we sample in logspace?
prob.append(Kroupa2001(options["M1"]) * options["M1"])
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob list ({})".format(
+ options["M1"], prob[-1], prob
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ """
+ From here we go through the multiplicities.
+ """
if multiplicity >= 2:
+ ###############
+ # Calculation for period of the binary
+
# Separation of the inner binary
options["sep"] = calc_sep_from_period(
options["M1"], options["M2"], options["P"]
)
+ # TODO: add check for min_P with instant RLOF?
# Total mass inner binary:
options["M1+M2"] = options["M1"] + options["M2"]
- # binary, triple or quadruple system
+ # 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),
)
- # TODO: fix this function
- # test the period distribution integrates to 1.0
- # if(0){
- # if(!defined $Moecache{'P_integrals'}->{$opts->{'M1'}})
- # {
- # my $I = 0.0;
- # my $dlogP = 1e-3;
- # for(my $logP = 0.0; $logP < 10.0; $logP += $dlogP)
- # {
- # my $dp_dlogP = # dp / dlogP
- # $Moecache{'rinterpolator_log10P'}->interpolate(
- # [
- # log10($opts->{'M1'}),
- # $logP
- # ])->[0];
- # $I += $dp_dlogP * $dlogP;
- # #printf "logM1 = %g, logP = %g -> dp/logP = %g\n",
- # # log10($opts->{'M1'}),
- # # $logP,
- # # $dp_dlogP;
- # }
- # $Moecache{'P_integrals'}->{$opts->{'M1'}} = $I;
-
- # #printf "M1=%g : P integral %g\n",
- # # $opts->{'M1'},
- # # $Moecache{'P_integrals'}->{$opts->{'M1'}};
- # }
- # }
+ if not Moecache.get("P_integrals", None):
+ Moecache["P_integrals"] = {}
- prob.append(
- Moecache["rinterpolator_log10P"].interpolate(
- [np.log10(options["M1"]), np.log10(options["P"])]
- )[0]
+ # Calculate normalisation contant for P1
+ if not Moecache["P_integrals"].get(options["M1"], None):
+ calc_P_integral(
+ options, 1, "P_integrals", "rinterpolator_log10P", "M1", verbosity
+ )
+
+ # Set probabilty for P1
+ p_val = Moecache["rinterpolator_log10P"].interpolate(
+ [np.log10(options["M1"]), np.log10(options["P"])]
+ )[0]
+ p_val = p_val / Moecache["P_integrals"][options["M1"]]
+ prob.append(p_val)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended period (m={}, P={}) probability ({}) to the prob list ({})".format(
+ options["M1"], options["P"], prob[-1], prob
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
)
############################################################
@@ -1591,32 +1740,78 @@ def Moe_de_Stefano_2017_pdf(options):
# we need to construct the q table for the given M1
# subject to qmin = Mmin/M1
- # Make a table storing Moe's data for q distributions
+ # Build the table for q
+ primary_mass = options["M1"]
+ secondary_mass = options["M2"]
+ m_label = "M1"
+ p_label = "P"
+ # 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,
)
- # Build the table for q
- primary_mass = options["M1"]
- secondary_mass = options["M2"]
- m_label = "M1"
- p_label = "P"
+ # 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(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
- build_q_table(options, m_label, p_label)
+ # Add probability for the mass ratio
prob.append(
Moecache[
"rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
].interpolate([secondary_mass / primary_mass])[0]
)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M1"],
+ options["P"],
+ options["M2"] / options["M1"],
+ prob[-1],
+ prob,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ ############################################################
+ # Eccentricity
# TODO add eccentricity calculations
+ # 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
+ 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,
+ )
+
+ # Calculations for when multiplicity is bigger than 3
+ # BEWARE: binary_c does not evolve these systems actually and the code below should be revised for when binary_c actually evolves triples.
+ # For that reason, I would not advise to use things with multiplicity > 3
if multiplicity >= 3:
- # triple or quadruple system
############################################################
# orbital period 2 =
@@ -1629,23 +1824,31 @@ def Moe_de_Stefano_2017_pdf(options):
# stars 1 and 2
# TODO: Is this a correct assumption?
- max_sep = (
- 10.0 * options["sep"] * (1.0 + options["ecc"])
- ) # TODO: isnt this the minimal separation?
+ max_sep = 10.0 * options["sep"] * (1.0 + options["ecc"])
min_P2 = calc_period_from_sep(options["M1+M2"], options["mmin"], max_sep)
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
prob.append(0)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: period2 is too short: {} < {}, system is not hierarchichal. Added 0 to probability list".format(
+ options["P1"], min_P2
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ # TODO: consider directly returning here
+
else:
# period is long enough that the system is hierarchical
-
# hence the separation between the outer star
# and inner binary
options["sep2"] = calc_sep_from_period(
options["M3"], options["M1+M2"], options["P2"]
)
+ # TODO: Consider actually removing these double interpolators.
+ # Not necessary to have these around, its the same dataset thats loaded
if not Moecache.get("rinterpolator_log10P2", None):
Moecache["rinterpolator_log10P2"] = py_rinterpolate.Rinterpolate(
table=Moecache[
@@ -1653,29 +1856,37 @@ def Moe_de_Stefano_2017_pdf(options):
], # 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),
)
+ # Create the dictionary where we store the area of the P2 integrals
if not Moecache.get("P2_integrals", None):
Moecache["P2_integrals"] = {}
+ # Calculate the P2 integrals if they don't exist
if not Moecache["P2_integrals"].get(options["M1+M2"], None):
- # normalize because $min_per > 0 so not all the periods
- # should be included
- I_p2 = 0
- dlogP2 = 1e-3
- for logP2 in np.arange(np.log10(min_P2), 10, dlogP2):
- # dp_dlogP2 = dp / dlogP
- dp_dlogP2 = Moecache["rinterpolator_log10P2"].interpolate(
- [np.log10(options["M1"]), logP2]
- )[0]
- I_p2 += dp_dlogP2 * dlogP2
- Moecache["P2_integrals"][options["M1+M2"]] = I_p2
+ calc_P_integral(
+ options,
+ min_P2,
+ "P2_integrals",
+ "rinterpolator_log10P2",
+ "M1+M2",
+ verbosity,
+ )
+ # Add the probability
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M1+M2"]), np.log10(options["P2"])]
)[0]
- p_val = p_val / Moecache["P2_integrals"]
+ p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
prob.append(p_val)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
+ options["M1"], options["M2"], options["P2"], prob[-1], prob
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
############################################################
# mass ratio 2 = q2 = M3 / (M1+M2)
@@ -1684,32 +1895,59 @@ def Moe_de_Stefano_2017_pdf(options):
# subject to qmin = Mmin/(M1+M2)
#
# Make a table storing Moe's data for q distributions
- # TODO: Check if this is correct.
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,
)
- # Build the table for q2
+ # Set the variables for the masses and their names
primary_mass = options["M1+M2"]
secondary_mass = options["M3"]
m_label = "M1+M2"
p_label = "P2"
- build_q_table(options, m_label, p_label)
+ # Build q table
+ build_q_table(options, m_label, p_label, verbosity=verbosity)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Called build_q_table",
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ # Add the probability
prob.append(
Moecache[
"rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
].interpolate([secondary_mass / primary_mass])[0]
)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M1+M2"],
+ options["M3"],
+ options["P"],
+ secondary_mass / primary_mass,
+ prob[-1],
+ prob,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
- # TODO: ecc2
+ # TODO: Implement ecc2 calculation
if multiplicity == 4:
# quadruple system.
- # TODO: Ask Rob about the strructure of the quadruple. Is htis only double binary quadrupples?
+ # TODO: Ask Rob about the strructure of the quadruple. Is this only double binary quadrupples?
############################################################
# orbital period 3
@@ -1718,7 +1956,8 @@ def Moe_de_Stefano_2017_pdf(options):
# as for any other stars but Pmax must be such that
# sep3 < sep2 * 0.2
- max_sep3 = 0.2 * options["sep2"] * (1.0 + opts["ecc2"])
+ # TODO: fix this here
+ max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
max_per3 = calc_period_from_sep(
options["M1+M2"], options["mmin"], max_sep3
)
@@ -1732,32 +1971,37 @@ def Moe_de_Stefano_2017_pdf(options):
], # 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("P2_integrals", None):
Moecache["P2_integrals"] = {}
if not Moecache["P2_integrals"].get(options["M1+M2"], None):
- # normalize because $min_per > 0 so not all the periods
- # should be included
- I_p2 = 0
- dlogP2 = 1e-3
-
- for logP2 in np.arange(np.log10(min_P2), 10, dlogP2):
- # dp_dlogP2 = dp / dlogP
- dp_dlogP2 = Moecache["rinterpolator_log10P2"].interpolate(
- [np.log10(options["M1"]), logP2]
- )[0]
- I_p2 += dp_dlogP2 * dlogP2
- Moecache["P2_integrals"][options["M1+M2"]] = I_p2
+ calc_P_integral(
+ options,
+ min_P2,
+ "P2_integrals",
+ "rinterpolator_log10P2",
+ "M1+M2",
+ verbosity,
+ )
# TODO: should this really be the log10P interpolator?
+ #
p_val = Moecache["rinterpolator_log10P"].interpolate(
[np.log10(options["M1+M2"]), np.log10(options["P2"])]
)[0]
-
- p_val = p_val / Moecache["P2_integrals"]
+ p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
prob.append(p_val)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
+ options["M1"], options["M2"], options["P2"], prob[-1], prob
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
############################################################
# mass ratio 2
@@ -1773,6 +2017,8 @@ def Moe_de_Stefano_2017_pdf(options):
], # Contains the table of data
nparams=3, # log10(M1+M2), log10P2, q
ndata=1, #
+ verbosity=verbosity
+ - (_MS_VERBOSITY_INTERPOLATOR_LEVEL - 1),
)
# Build the table for q2
@@ -1781,79 +2027,101 @@ def Moe_de_Stefano_2017_pdf(options):
m_label = "M1+M2"
p_label = "P2"
- build_q_table(options, m_label, p_label)
+ build_q_table(options, m_label, p_label, verbosity=verbosity)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Created q_table ".format(),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ # Add the probability
prob.append(
Moecache[
"rinterpolator_q_given_{}_log10{}".format(m_label, p_label)
].interpolate([secondary_mass / primary_mass])[0]
)
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+ options["M1+M2"],
+ options["M3"],
+ options["P"],
+ secondary_mass / primary_mass,
+ prob[-1],
+ prob,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# todo ecc 3
elif multiplicity not in range(1, 5):
- print("Unknown multiplicity {}\n".format(multiplicity))
-
- # TODO: Translate this to get some better output
- # printf "PDFPROBS (M1=%g, logP=%g q=%g) : @prob\n",
- # $opts->{'M1'},
- # $multiplicity == 2 ? ($opts->{'M2'}/$opts->{'M1'}) : -999,
- # $multiplicity == 2 ? log10($opts->{'P'}) : -999;
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: Unknown multiplicity {}".format(
+ multiplicity
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
# the final probability density is the product of all the
# probability density functions
- # TODO: Where do we take into account the stepsize? probdens is not a probability if we dont take a stepsize
prob_dens = 1 * np.prod(prob)
- print_info = 1
- if print_info:
- # print("Probability density")
-
- if multiplicity == 1:
- print(
- "M1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
- options["M1"], len(prob), str(prob), prob_dens
- )
- )
- elif multiplicity == 2:
- print(
- "M1={} q={} log10P={} ({}): {} -> {}\n".format(
- options["M1"],
- options["M2"] / options["M1"],
- np.log10(options["P"]),
- len(prob),
- str(prob),
- prob_dens,
- )
- )
- elif multiplicity == 3:
- print(
- "M1={} q={} log10P={} ({}): M3={}} P2={} ecc2={} : {} - {}".format(
- options["M1"],
- options["M2"] / options["M1"],
- np.log10(options["P"]),
- len(prob),
- options["M3"],
- np.log10("P2"),
- np.log10("ecc2"),
- str(prob),
- prob_dens,
- )
- )
- elif multiplicity == 4:
- print(
- "M1={} q={} log10P={} ({}) : M3={} P2={}} ecc2={} : M4={} P3={} ecc3={} : {} -> {}".format(
- options["M1"],
- options["M2"] / options["M1"],
- np.log10(options["P"]),
- len(prob),
- options["M3"],
- np.log10("P2"),
- np.log10(options["ecc2"]),
- options["M4"],
- np.log10("P3"),
- np.log10(options["ecc3"]),
- str(prob),
- prob_dens,
- )
- )
+ if multiplicity == 1:
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
+ options["M1"], len(prob), str(prob), prob_dens
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ elif multiplicity == 2:
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}): {} -> {}\n".format(
+ options["M1"],
+ options["M2"] / options["M1"],
+ np.log10(options["P"]),
+ len(prob),
+ str(prob),
+ prob_dens,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ elif multiplicity == 3:
+ verbose_print(
+ "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}): M3={}} P2={} ecc2={} : {} - {}".format(
+ options["M1"],
+ options["M2"] / options["M1"],
+ np.log10(options["P"]),
+ len(prob),
+ options["M3"],
+ np.log10("P2"),
+ np.log10("ecc2"),
+ str(prob),
+ prob_dens,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
+ elif multiplicity == 4:
+ verbose_print(
+ "M&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ({}) : M3={} P2={}} ecc2={} : M4={} P3={} ecc3={} : {} -> {}".format(
+ options["M1"],
+ options["M2"] / options["M1"],
+ np.log10(options["P"]),
+ len(prob),
+ options["M3"],
+ np.log10("P2"),
+ np.log10(options["ecc2"]),
+ options["M4"],
+ np.log10("P3"),
+ np.log10(options["ecc3"]),
+ str(prob),
+ prob_dens,
+ ),
+ verbosity,
+ _MS_VERBOSITY_LEVEL,
+ )
return prob_dens
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 3a054a3a0431b5a62613ba68b11bc8dc01ceea7e..0df3abe35e29c35b7ab96d72932c55bf22ffa44d 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -26,6 +26,8 @@ from collections import defaultdict
from binarycpython import _binary_c_bindings
import binarycpython.utils.moe_distefano_data as moe_distefano_data
+import py_rinterpolate
+
########################################################
# Unsorted
########################################################
@@ -1705,6 +1707,8 @@ def binaryc_json_serializer(obj: Any) -> Any:
if inspect.isfunction(obj):
return str(obj)
+ elif isinstance(obj, py_rinterpolate.Rinterpolate):
+ return str(obj)
else:
return obj
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 9911206b5fbe667ee471bb29409897abc641c9ad..7bd10cc05e9b8e9fbc28785b5838d350a85513e7 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -43,6 +43,7 @@ from pathos.pools import _ProcessPool as Pool
from binarycpython.utils.grid_options_defaults import (
grid_options_defaults_dict,
moe_distefano_default_options,
+ _MS_VERBOSITY_LEVEL,
)
from binarycpython.utils.custom_logging_functions import (
@@ -64,6 +65,7 @@ from binarycpython.utils.functions import (
get_moe_distefano_dataset,
)
+import py_rinterpolate
# from binarycpython.utils.hpc_functions import (
# get_condor_version,
@@ -76,6 +78,7 @@ from binarycpython.utils.functions import (
from binarycpython.utils.distribution_functions import (
Moecache,
LOG_LN_CONVERTER,
+ fill_data,
)
from binarycpython import _binary_c_bindings
@@ -743,11 +746,8 @@ class Population:
],
}
- # print(Moecache)
- print(Moecache.keys())
-
##
- # Clean up code: remove files, unset values. This is placed in the general evolve function,
+ # Clean up code: remove files, unset values, unload interpolators etc. This is placed in the general evolve function,
# because that makes for easier control
self._cleanup()
@@ -928,7 +928,7 @@ class Population:
# pathos_multiprocess
manager = multiprocessing.Manager()
- job_queue = manager.Queue(maxsize=self.grid_options['max_queue_size'])
+ job_queue = manager.Queue(maxsize=self.grid_options["max_queue_size"])
result_queue = manager.Queue(maxsize=self.grid_options["amt_cores"])
# Create process instances
@@ -958,7 +958,6 @@ class Population:
sentinel = object()
for output_dict in iter(result_queue.get, sentinel):
- print(output_dict)
combined_output_dict = merge_dicts(combined_output_dict, output_dict)
if result_queue.empty():
break
@@ -1358,6 +1357,16 @@ class Population:
result_queue.put(output_dict)
stream_logger.debug(f"Process-{self.process_ID} is finished.")
+ # Clean up the interpolators if they exist
+ # TODO: make a cleanup function for the individual threads
+ # TODO: make sure this is necessary. Actually its probably not, because we have a centralized queue
+ verbose_print(
+ "Process {}: Cleaning up interpolators".format(ID),
+ self.grid_options["verbosity"],
+ 1,
+ )
+ self._clean_interpolators()
+
return
# Single system
@@ -1597,6 +1606,9 @@ class Population:
# Unload/free custom_logging_code
# TODO: cleanup custom logging code.
+ # Clean up the interpolator functions
+ self._clean_interpolators()
+
###################################################
# Gridcode functions
#
@@ -1800,7 +1812,7 @@ class Population:
# Add condition failed action:
code_string += (
indent * (depth + 2)
- + 'print("Condition for {} not met!")'.format(
+ + 'print("Grid generator: Condition for {} not met!")'.format(
grid_variable["parameter_name"]
)
+ "\n"
@@ -1821,7 +1833,9 @@ class Population:
# Add phasevol check action:
code_string += (
indent * (depth + 2)
- + 'print("phasevol_{} <= 0!")'.format(grid_variable["name"])
+ + 'print("Grid generator: phasevol_{} <= 0!")'.format(
+ grid_variable["name"]
+ )
+ "\n"
)
code_string += indent * (depth + 2) + "continue" + "\n"
@@ -3079,8 +3093,17 @@ class Population:
# Require input options
if not options:
- print("Please provide a options dictionary.")
- raise ValueError
+ 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)
@@ -3089,7 +3112,7 @@ class Population:
options = update_dicts(default_options, options)
# Write options to a file
- with open("/tmp/moeopts.dat", "w") as f:
+ with open(os.path.join(ms_tmp_dir, "moeopts.dat"), "w") as f:
f.write(json.dumps(options, indent=4))
json_data = get_moe_distefano_dataset(options)
@@ -3100,9 +3123,16 @@ class Population:
# Get all the masses
logmasses = sorted(json_data["log10M1"].keys())
if not logmasses:
- print("The table does not contain masses")
- raise ValueError
- with open("/tmp/moe.log", "w") as logfile:
+ 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))
# Get all the periods and see if they are all consistently present
@@ -3115,8 +3145,13 @@ class Population:
current_logperiods = sorted(json_data["log10M1"][logmass]["logP"])
if not (logperiods == current_logperiods):
- print("Period values are not consistent throughout the dataset")
- raise ValueError
+ msg = "Period values are not consistent throughout the dataset"
+ verbose_print(
+ "\tMoe_de_Stefano_2017: {}".format(msg),
+ self.grid_options["verbosity"],
+ 0,
+ )
+ raise ValueError(msg)
############################################################
# log10period binwidth : of course this assumes a fixed
@@ -3126,16 +3161,22 @@ class Population:
if not dlog10P == (
float(current_logperiods[i + 1]) - float(current_logperiods[i])
):
- print("Period spacing is not consistent throughout the dataset")
- raise ValueError
+ 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)
- with open("/tmp/moe.log", "a") as logfile:
+ # 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:
- # print("Logmass: {}".format(logmass))
-
# Create the multiplicity table
if not Moecache.get("multiplicity_table", None):
Moecache["multiplicity_table"] = []
@@ -3151,7 +3192,6 @@ class Population:
]
)
- print("Size multiplicity table: {}", len(Moecache["multiplicity_table"]))
############################################################
# a small log10period which we can shift just outside the
# table to force integration out there to zero
@@ -3216,10 +3256,6 @@ class Population:
/ dlog10P,
]
)
- print(
- "Size period_distributions table: {}",
- len(Moecache["period_distributions"]),
- )
############################################################
# distributions as a function of mass, period, q
@@ -3227,18 +3263,11 @@ class Population:
# First, get a list of the qs given by Moe
#
qs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["q"])
- qdata = {}
- I_q = 0.0
- # get the q distribution from M&S, which ranges 0.15 to 0.95
- for q in qs:
- val = json_data["log10M1"][logmass]["logP"][logperiod]["q"][q]
- qdata[q] = val
- I_q += val
-
- # hence normalize to 1.0
- for q in qs:
- qdata[q] = qdata[q] / I_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):
@@ -3252,18 +3281,11 @@ class Population:
############################################################
# eccentricity distributions as a function of mass, period, ecc
eccs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["e"])
- ecc_data = {}
- I_e = 0
-
- # Read out the data
- for ecc in eccs:
- val = json_data["log10M1"][logmass]["logP"][logperiod]["e"][ecc]
- ecc_data[ecc] = val
- I_e += val
- # Normalize the data
- for ecc in eccs:
- ecc_data[ecc] = ecc_data[ecc] / I_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):
@@ -3279,7 +3301,7 @@ class Population:
Moecache["period_distributions"].append(
[
float(logmass),
- float(logperiods[-1]) + 0.5 * dlog10P,
+ float(logperiods[-1]) + 0.5 * dlog10P, # TODO: why this shift?
json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_bin_frac_p_dist"
]
@@ -3298,14 +3320,38 @@ class Population:
0.0,
]
)
- print(
- "Size period_distributions table: {}".format(
- len(Moecache["period_distributions"])
- )
- )
+
+ 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("/tmp/moecache.json", "w") as cache_filehandle:
+ with open(os.path.join(ms_tmp_dir, "moecache.json"), "w") as cache_filehandle:
cache_filehandle.write(json.dumps(Moecache, indent=4))
# TODO: remove this and make the plotting of the multiplicity fractions more modular
@@ -3349,7 +3395,6 @@ class Population:
############################################################
# make a version of the options which can be passed to
# the probability density function in the gridcode
- # TODO: check if that works, if we just pass the options as a dict
############################################################
# first, the multiplicity, this is 1,2,3,4, ...
@@ -3357,11 +3402,13 @@ class Population:
max_multiplicity = 0
for i in range(1, 5):
mod = options["multiplicity_modulator"][i - 1]
- print("Check mod {} = {}".format(i, mod))
if mod > 0:
max_multiplicity = i
- print("Max multiplicity = {}".format(max_multiplicity))
-
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Max multiplicity = {}".format(max_multiplicity),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
+ )
######
# Setting up the grid variables
@@ -3408,7 +3455,7 @@ class Population:
gridtype="centred",
dphasevol="dlnm1",
precode='M_1 = np.exp(lnm1); options["M1"]=M_1',
- probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}) if multiplicity == 1 else 1".format(
+ probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity']) if multiplicity == 1 else 1".format(
str(options)[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
),
)
@@ -3677,7 +3724,7 @@ 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({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}})".format(
+ probdist_addition = "Moe_de_Stefano_2017_pdf({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}, verbosity=self.grid_options['verbosity'])".format(
str(options)[1:-1],
'"multiplicity": multiplicity',
'"M1": M_1',
@@ -3696,3 +3743,28 @@ eccentricity3=0
self.grid_options["_grid_variables"][self.last_grid_variable()][
"probdist"
] = probdist_addition
+
+ verbose_print(
+ "\tMoe_de_Stefano_2017: Added final call to the pdf function".format(
+ max_multiplicity
+ ),
+ self.grid_options["verbosity"],
+ _MS_VERBOSITY_LEVEL,
+ )
+
+ def _clean_interpolators(self):
+ """
+ Function to clean up the interpolators after a run
+
+ We look in the MoeCache global variable for items that are interpolators.
+ Should be called by the general cleanup function AND the thread cleanup function
+ """
+
+ interpolator_keys = []
+ for key in Moecache.keys():
+ if isinstance(Moecache[key], py_rinterpolate.Rinterpolate):
+ interpolator_keys.append(key)
+
+ for key in interpolator_keys:
+ Moecache[key].destroy()
+ del Moecache[key]
diff --git a/binarycpython/utils/grid_options_defaults.py b/binarycpython/utils/grid_options_defaults.py
index 98bf8b72ed63f2bae3a660e6bb5b3b8d4554f9d0..971a13ecd55f5fd6f242f4cb1615e81f5cb1c2fd 100644
--- a/binarycpython/utils/grid_options_defaults.py
+++ b/binarycpython/utils/grid_options_defaults.py
@@ -19,6 +19,11 @@ import os
from binarycpython.utils.custom_logging_functions import temp_dir
+_MS_VERBOSITY_LEVEL = 5
+_MS_VERBOSITY_INTERPOLATOR_LEVEL = 6
+_MS_VERBOSITY_INTERPOLATOR_EXTRA_LEVEL = 7
+
+
# Options dict
grid_options_defaults_dict = {
##########################
@@ -34,7 +39,7 @@ grid_options_defaults_dict = {
"_commandline_input": "",
"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!
+ "max_queue_size": 1000, # Maximum size of the system call queue. Can't be too big!
##########################
# Execution log:
##########################
@@ -472,7 +477,7 @@ grid_options_descriptions = {
"_total_mass_run": "To count the total mass that thread/process has ran",
"_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"
+ "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",
}
###
@@ -579,8 +584,8 @@ moe_distefano_default_options = {
# linear2
# plaw2
# nolowq
- "q_low_extrapolation_method": "linear2",
- "q_high_extrapolation_method": "linear2",
+ "q_low_extrapolation_method": "linear",
+ "q_high_extrapolation_method": "linear",
}
moe_distefano_default_options_description = {