diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 807f938bfc1d1073147a6b69cd70276fd4929d37..e8e0176a5fefd1185b4dc200cc368ec5c8b3de92 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -2007,8 +2007,8 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
# Calculate P integral or use cached value
# get the periods from the Moecahe
- min_logP = Moecache['logperiods'][0]
- max_logP = Moecache['logperiods'][-1]
+ min_logP = float(Moecache['logperiods'][0])
+ max_logP = float(Moecache['logperiods'][-1])
calc_P_integral(
options, min_logP, max_logP, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
@@ -2135,7 +2135,7 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
min_logP2 = math.log10(min_P2)
#max_logP2 = 10.0
#min_logP = Moecache['logperiods'][0]
- max_logP2 = Moecache['logperiods'][-1]
+ max_logP2 = float(Moecache['logperiods'][-1])
if options["P2"] < min_P2:
# period is too short : system is not hierarchical
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index d615cdab21e1138e78ab1c968027e890a527dd5c..75998803213d3d0daca6963adb8bcf5dc218598d 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -39,7 +39,7 @@ import strip_ansi
import sys
import time
import uuid
-
+_count = 0
from typing import Union, Any
from collections import (
OrderedDict,
@@ -475,19 +475,22 @@ class Population:
def add_grid_variable(
self,
name: str,
+ parameter_name: str,
longname: str,
valuerange: Union[list, str],
resolution: str,
spacingfunc: str,
probdist: str,
dphasevol: Union[str, int],
- parameter_name: str,
gridtype: str = "centred",
branchpoint: int = 0,
branchcode: Union[str, None] = None,
precode: Union[str, None] = None,
postcode: Union[str, None] = None,
+ topcode: Union[str, None] = None,
+ bottomcode: Union[str, None] = None,
condition: Union[str, None] = None,
+
) -> None:
"""
Function to add grid variables to the grid_options.
@@ -505,11 +508,23 @@ class Population:
Args:
name:
- name of parameter. This is evaluated as a parameter and you can use it throughout
+ name of parameter used in the grid Python code.
+ This is evaluated as a parameter and you can use it throughout
the rest of the function
Examples:
name = 'lnm1'
+
+ parameter_name:
+ name of the parameter in binary_c
+
+ This name must correspond to a Python variable of the same name,
+ which is automatic if parameter_name == name.
+
+ Note: if parameter_name != name, you must set a
+ variable in "precode" or "postcode" to define a Python variable
+ called parameter_name
+
longname:
Long name of parameter
@@ -561,11 +576,18 @@ class Population:
Method on how the value range is sampled. Can be either 'edge' (steps starting at
the lower edge of the value range) or 'centred'
(steps starting at lower edge + 0.5 * stepsize).
+
+ topcode:
+ Code added at the very top of the block.
+
+ bottomcode:
+ Code added at the very bottom of the block.
"""
# Add grid_variable
grid_variable = {
"name": name,
+ "parameter_name":parameter_name,
"longname": longname,
"valuerange": valuerange,
"resolution": resolution,
@@ -574,11 +596,12 @@ class Population:
"postcode": postcode,
"probdist": probdist,
"dphasevol": dphasevol,
- "parameter_name": parameter_name,
"condition": condition,
"gridtype": gridtype,
"branchpoint": branchpoint,
"branchcode": branchcode,
+ "topcode" : topcode,
+ "bottomcode" : bottomcode,
"grid_variable_number": len(self.grid_options["_grid_variables"]),
}
@@ -589,6 +612,7 @@ class Population:
# Load it into the grid_options
self.grid_options["_grid_variables"][grid_variable["name"]] = grid_variable
+
verbose_print(
"Added grid variable: {}".format(json.dumps(grid_variable, indent=4)),
self.grid_options["verbosity"],
@@ -1365,6 +1389,7 @@ class Population:
# Have some user-defined function do stuff with the data.
if self.grid_options["parse_function"]:
+ self.custom_options["parameter_dict"] = full_system_dict
self.grid_options["parse_function"](self, out)
def _process_run_population_grid(self, job_queue, result_queue, ID):
@@ -1625,12 +1650,12 @@ class Population:
# Debug line: logging all the lines
if self.grid_options["log_runtime_systems"] == 1:
with open(
- os.path.join(
- self.grid_options["tmp_dir"],
- "runtime_systems",
- "process_{}.txt".format(self.process_ID),
- ),
- "a+",
+ os.path.join(
+ self.grid_options["tmp_dir"],
+ "runtime_systems",
+ "process_{}.txt".format(self.process_ID),
+ ),
+ "a+",
) as f:
binary_cmdline_string = self._return_argline(full_system_dict)
f.write(
@@ -1664,12 +1689,12 @@ class Population:
# Set status to finishing
with open(
- os.path.join(
- self.grid_options["tmp_dir"],
- "process_status",
- "process_{}.txt".format(self.process_ID),
- ),
- "w",
+ os.path.join(
+ self.grid_options["tmp_dir"],
+ "process_status",
+ "process_{}.txt".format(self.process_ID),
+ ),
+ "w",
) as f:
f.write("FINISHING")
f.close()
@@ -1756,23 +1781,23 @@ class Population:
colour = 'cyan on black'
verbose_print(
self._boxed(
- "{colour}Process {ID} finished:\ngenerator started at {start}\ngenerator finished at {end}\ntotal: {timesecs}\nof which {binary_c_secs} with binary_c\nRan {nsystems} systems\nwith a total probability of {psystems:g}\n{failcolour}This thread had {nfail} failing systems{colour}\n{failcolour}with a total failed probability of {pfail}{colour}\n{zerocolour}Skipped a total of {nzero} zero-probability systems{zeroreset}\n".format(
- colour=self.ANSI_colours[colour],
- ID=ID,
- start=start_process_time.isoformat(),
- end=end_process_time.isoformat(),
- timesecs=timedelta((end_process_time - start_process_time).total_seconds()),
- binary_c_secs=timedelta(total_time_calling_binary_c),
- nsystems=number_of_systems_run,
- psystems=probability_of_systems_run,
- failcolour=self.ANSI_colours['red'] if self.grid_options["_failed_count"]>0 else '',
- failreset=self.ANSI_colours[colour] if self.grid_options["_failed_count"]>0 else '',
- nfail=self.grid_options["_failed_count"],
- pfail=self.grid_options["_failed_prob"],
- nzero=zero_prob_stars_skipped,
- zerocolour=self.ANSI_colours["yellow"] if zero_prob_stars_skipped>0 else '',
- zeroreset=self.ANSI_colours[colour] if zero_prob_stars_skipped>0 else ''
- ),colour=colour),
+ "{colour}Process {ID} finished:\ngenerator started at {start}\ngenerator finished at {end}\ntotal: {timesecs}\nof which {binary_c_secs} with binary_c\nRan {nsystems} systems\nwith a total probability of {psystems:g}\n{failcolour}This thread had {nfail} failing systems{colour}\n{failcolour}with a total failed probability of {pfail}{colour}\n{zerocolour}Skipped a total of {nzero} zero-probability systems{zeroreset}\n".format(
+ colour=self.ANSI_colours[colour],
+ ID=ID,
+ start=start_process_time.isoformat(),
+ end=end_process_time.isoformat(),
+ timesecs=timedelta((end_process_time - start_process_time).total_seconds()),
+ binary_c_secs=timedelta(total_time_calling_binary_c),
+ nsystems=number_of_systems_run,
+ psystems=probability_of_systems_run,
+ failcolour=self.ANSI_colours['red'] if self.grid_options["_failed_count"]>0 else '',
+ failreset=self.ANSI_colours[colour] if self.grid_options["_failed_count"]>0 else '',
+ nfail=self.grid_options["_failed_count"],
+ pfail=self.grid_options["_failed_prob"],
+ nzero=zero_prob_stars_skipped,
+ zerocolour=self.ANSI_colours["yellow"] if zero_prob_stars_skipped>0 else '',
+ zeroreset=self.ANSI_colours[colour] if zero_prob_stars_skipped>0 else ''
+ ),colour=colour),
self.grid_options["verbosity"],
1,
)
@@ -1794,24 +1819,24 @@ class Population:
"zero_prob_stars_skipped": zero_prob_stars_skipped,
}
with open(
- os.path.join(
- self.grid_options["tmp_dir"],
- "process_summary",
- "process_{}.json".format(self.process_ID),
- ),
- "w",
+ os.path.join(
+ self.grid_options["tmp_dir"],
+ "process_summary",
+ "process_{}.json".format(self.process_ID),
+ ),
+ "w",
) as f:
f.write(json.dumps(summary_dict, indent=4))
f.close()
# Set status to running
with open(
- os.path.join(
- self.grid_options["tmp_dir"],
- "process_status",
- "process_{}.txt".format(self.process_ID),
- ),
- "w",
+ os.path.join(
+ self.grid_options["tmp_dir"],
+ "process_status",
+ "process_{}.txt".format(self.process_ID),
+ ),
+ "w",
) as f:
f.write("FINISHED")
f.close()
@@ -1936,7 +1961,7 @@ class Population:
raise ValueError
if not any(
- [key.startswith("ensemble_filter_") for key in self.bse_options]
+ [key.startswith("ensemble_filter_") for key in self.bse_options]
):
verbose_print(
"Warning: Running the ensemble without any filter requires a lot of available RAM",
@@ -2092,7 +2117,7 @@ class Population:
# a variable grid
###################################################
- def _add_code(self, *args, indent=0, mindent=-1):
+ def _add_code(self, *args, indent=0):
"""
Function to add code to the grid code string
@@ -2104,21 +2129,9 @@ class Population:
don't use both!
"""
- # if we have multiple arguments and have forgotten
- # to set mindent = 0, do it here
- if indent == 0 and mindent == -1 and len(args) > 0:
- mindent = 0
-
- # if mindent isn't set, use indent
- if mindent == -1:
- self.code_string += self._indent_block(indent)
-
- # put everything in the code_string, indenting
- # each line if required
+ indent_block = self._indent_block(indent)
for thing in args:
- if mindent != -1 or len(args) > 1:
- self.code_string += self._indent_block(mindent)
- self.code_string += thing
+ self.code_string += indent_block + thing
def _indent_block(self, n=0):
"""
@@ -2192,7 +2205,6 @@ class Population:
"probabilities_sum = [0 for i in range({})]\n".format(total_grid_variables),
"parameter_dict = {}\n",
"phasevol = 1\n",
- "\n",
)
# Set up the system parameters
@@ -2213,13 +2225,13 @@ class Population:
)
for grid_variable_el in sorted(
- self.grid_options["_grid_variables"].items(),
- key=lambda x: x[1]["grid_variable_number"],
+ self.grid_options["_grid_variables"].items(),
+ key=lambda x: x[1]["grid_variable_number"],
):
# Make probabilities dict
grid_variable = grid_variable_el[1]
self._add_code(
- 'probabilities["{}"] = 0\n'.format(grid_variable["parameter_name"])
+ 'probabilities["{}"] = 0\n'.format(grid_variable["name"])
)
#################################################################################
@@ -2233,10 +2245,10 @@ class Population:
# Generate code
print("Generating grid code")
for loopnr, grid_variable_el in enumerate(
- sorted(
- self.grid_options["_grid_variables"].items(),
- key=lambda x: x[1]["grid_variable_number"],
- )
+ sorted(
+ self.grid_options["_grid_variables"].items(),
+ key=lambda x: x[1]["grid_variable_number"],
+ )
):
verbose_print(
"Constructing/adding: {}".format(grid_variable_el[0]),
@@ -2244,11 +2256,18 @@ class Population:
2)
grid_variable = grid_variable_el[1]
+ ####################
+ # top code
+ if grid_variable["topcode"]:
+ self._add_code(
+ grid_variable["topcode"]
+ )
+
#########################
# Setting up the for loop
# Add comment for for loop
self._add_code(
- "# for loop for {}".format(grid_variable["parameter_name"]) + "\n",
+ "# for loop for {}".format(grid_variable["name"]) + "\n",
"sampled_values_{} = {}".format(
grid_variable["name"], grid_variable["spacingfunc"]
)
@@ -2312,14 +2331,16 @@ class Population:
self._add_code("else:\n")
self._add_code(
"{}_this = {}_sample_number; ".format(
- grid_variable["name"], grid_variable["name"]
+ grid_variable["name"],
+ grid_variable["name"]
),
indent=1,
)
self._add_code(
"\n",
"{}_next = {}_this + 1".format(
- grid_variable["name"], grid_variable["name"]
+ grid_variable["name"],
+ grid_variable["name"]
),
"\n",
)
@@ -2329,7 +2350,7 @@ class Population:
# 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
+ not grid_variable["dphasevol"] == -1
): # A method to turn off this calculation and allow a phasevol = 1
self._add_code(
"phasevol_{} = sampled_values_{}[{}_next]-sampled_values_{}[{}_this]".format(
@@ -2356,8 +2377,8 @@ class Population:
grid_variable["name"], grid_variable["name"], grid_variable["name"]
)
+ "\n",
- "continue" + "\n",
- mindent=1,
+ "continue\n",
+ indent=1,
)
if vb:
@@ -2373,9 +2394,9 @@ class Population:
# select point location based on gridtype (left, centre or right)
if (
- grid_variable["gridtype"] == "edge"
- or grid_variable["gridtype"] == "left"
- or grid_variable["gridtype"] == "left edge"
+ grid_variable["gridtype"] == "edge"
+ or grid_variable["gridtype"] == "left"
+ or grid_variable["gridtype"] == "left edge"
):
self._add_code(
"{} = sampled_values_{}[{}_this]".format(
@@ -2386,8 +2407,8 @@ class Population:
+ "\n"
)
elif (
- grid_variable["gridtype"] == "right"
- or grid_variable["gridtype"] == "right edge"
+ grid_variable["gridtype"] == "right"
+ or grid_variable["gridtype"] == "right edge"
):
self._add_code(
+"{} = sampled_values_{}[{}_next]".format(
@@ -2398,9 +2419,9 @@ class Population:
+ "\n"
)
elif (
- grid_variable["gridtype"] == "centred"
- or grid_variable["gridtype"] == "centre"
- or grid_variable["gridtype"] == "center"
+ grid_variable["gridtype"] == "centred"
+ or grid_variable["gridtype"] == "centre"
+ or grid_variable["gridtype"] == "center"
):
self._add_code(
"{} = 0.5 * (sampled_values_{}[{}_next] + sampled_values_{}[{}_this])".format(
@@ -2431,28 +2452,28 @@ class Population:
self._add_code(
# Add comment
- "# Condition for {}\n".format(grid_variable["parameter_name"]),
+ "# Condition for {}\n".format(grid_variable["name"]),
# Add condition check
"if not {}:\n".format(grid_variable["condition"]),
- mindent=0,
+ indent=0,
)
# Add condition failed action:
if self.grid_options['verbosity'] >= 3:
self._add_code(
'print("Grid generator: Condition for {} not met!")'.format(
- grid_variable["parameter_name"]
+ grid_variable["name"]
)
+ "\n",
"continue" + "\n",
- mindent=1,
+ indent=1,
)
else:
self._add_code(
"continue" + "\n",
- mindent=1,
+ indent=1,
)
- # Add some whitespace
+ # Add some whitespace
self._add_code("\n")
# # Add debugging:
@@ -2539,7 +2560,7 @@ class Population:
self._add_code(
"\n",
"# Increment starcount for {}\n".format(
- grid_variable["parameter_name"]
+ grid_variable["name"]
),
"starcounts[{}] += 1".format(
grid_variable["grid_variable_number"],
@@ -2547,7 +2568,8 @@ class Population:
+ "\n",
# Add value to dict
'parameter_dict["{}"] = {}'.format(
- grid_variable["parameter_name"], grid_variable["parameter_name"]
+ grid_variable["parameter_name"],
+ grid_variable["parameter_name"]
)
+ "\n",
"\n",
@@ -2570,6 +2592,13 @@ class Population:
# increment indent_depth
self._increment_indent_depth(+1)
+ ####################
+ # bottom code
+ if grid_variable["bottomcode"]:
+ self._add_code(
+ grid_variable["bottomcode"]
+ )
+
self._increment_indent_depth(-1)
self._add_code("\n")
@@ -2609,7 +2638,7 @@ class Population:
self._add_code(
# Add comment
"# Condition for branchpoint at {}".format(
- reverse_sorted_grid_variables[loopnr + 1][1]["parameter_name"]
+ reverse_sorted_grid_variables[loopnr + 1][1]["name"]
)
+ "\n",
# # Add condition check
@@ -2635,7 +2664,7 @@ class Population:
self._add_code("\n", "#" * 40 + "\n", "if print_results:\n")
self._add_code(
"print('Grid has handled {} stars with a total probability of {}'.format(_total_starcount,self.grid_options['_probtot']))\n",
- mindent=1,
+ indent=1,
)
################
@@ -2676,8 +2705,10 @@ class Population:
# perhaps create symlink
if self.grid_options['symlink latest gridcode']:
+ global _count
symlink = os.path.join(self.grid_options['tmp_dir'],
- 'binary_c_grid-latest')
+ 'binary_c_grid-latest' + str(_count))
+ _count += 1
if os.path.exists(symlink):
os.unlink(symlink)
@@ -2746,7 +2777,7 @@ class Population:
'parameter_dict["{}"] = {}'.format("phasevol", "phasevol") + "\n",
# Increment total probability
"self._increment_probtot(probability)\n",
- mindent=1,
+ indent=1,
)
if not dry_run:
@@ -3401,7 +3432,7 @@ class Population:
def write_ensemble(self,output_file,json_data):
"""
- write_ensemble : Write ensemble results to a file.
+ write_ensemble : Write ensemble results to a file.
Args:
output_file : the output filename.
@@ -3472,10 +3503,10 @@ class Population:
############################################################
def write_binary_c_calls_to_file(
- self,
- output_dir: Union[str, None] = None,
- output_filename: Union[str, None] = None,
- include_defaults: bool = False,
+ self,
+ output_dir: Union[str, None] = None,
+ output_filename: Union[str, None] = None,
+ include_defaults: bool = False,
) -> None:
"""
Function that loops over the grid code and writes the generated parameters to a file.
@@ -3531,7 +3562,7 @@ class Population:
# Check if there is an output dir configured
if self.custom_options.get("data_dir", None):
binary_c_calls_output_dir = self.custom_options["data_dir"]
- # otherwise check if there's one passed to the function
+ # otherwise check if there's one passed to the function
else:
if not output_dir:
print(
@@ -3543,7 +3574,7 @@ class Population:
# check if there's a filename passed to the function
if output_filename:
binary_c_calls_filename = output_filename
- # otherwise use default value
+ # otherwise use default value
else:
binary_c_calls_filename = "binary_c_calls.txt"
diff --git a/binarycpython/utils/spacing_functions.py b/binarycpython/utils/spacing_functions.py
index 563ca859e4e373b00f3d54ab1c880bfbceb5f016..1df93cb4316c0ab55997a916c3919d98c396b004 100644
--- a/binarycpython/utils/spacing_functions.py
+++ b/binarycpython/utils/spacing_functions.py
@@ -144,6 +144,9 @@ def const_dt(self,
maxdm=((0.07,1.0,0.1),(1.0,300.0,1.0)), # tuple of tuples
fsample=1.0,
factor=1.0,
+ logmasses=False,
+ log10masses=False,
+ showlist=False,
showtable=False):
"""
const_dt returns a list of masses spaced at a constant age difference
@@ -157,11 +160,14 @@ def const_dt(self,
logspacing: whether to use log-spaced time, in which case dt is actually d(log10(t))
tmin: the minimum time to consider (Myr, default 3.0 Myr)
tmax: the maximum time to consider (Myr, default None which means we use the grid option 'max_evolution_time')
- mindm: a tuple of tuples containing a mass range and minimum mass spacing in that range. The default is ((0.07,1.0,0.1),(1.0,300.0,1.0)) allocated a minimum dm of 0.1Msun in the mass range 0.07 to 1.0 Msun and 1.0Msun in the range 1.0 to 300.0 Msun. Anything you set overrides this.
- maxdm: a list of tuples similar to mindm but specifying a maximum mass spacing. (None)
+ mindm: a tuple of tuples containing a mass range and minimum mass spacing in that range. The default is ((0.07,1.0,0.1),(1.0,300.0,1.0)) allocated a minimum dm of 0.1Msun in the mass range 0.07 to 1.0 Msun and 1.0Msun in the range 1.0 to 300.0 Msun. Anything you set overrides this. Note, if you use only one tuple, you must set it with a trailing comma, thus, e.g. ((0.07,1.0,0.1),). (default None)
+ maxdm: a list of tuples similar to mindm but specifying a maximum mass spacing. In the case of maxdm, if the third option in each tuple is negative it is treated as a log step (its absolute value is used as the step). (default None)
fsample: a global sampling (Shannon-like) factor (<1) to improve resolution (default 1.0, set to smaller to improve resolution)
factor: all masses generated are multiplied by this after generation
showtable: if True, the mass list and times are shown to stdout after generation
+ showlist: if True, show the mass list once generated
+ logmasses: if True, the masses are logged with math.log()
+ log10masses: if True, the masses are logged with math.log10()
Returns:
Array of masses.
@@ -197,13 +203,13 @@ def const_dt(self,
# make a grid in M1
lifetime_population.add_grid_variable(
name="lnM_1",
+ parameter_name="M_1",
longname="log Primary mass", # == single-star mass
valuerange=[math.log(mmin),math.log(mmax)],
resolution=0,
spacingfunc="const(math.log({mmin}),math.log({mmax}),{nres})".format(mmin=mmin,mmax=mmax,nres=nres),
probdist="1", # dprob/dm1 : we don't care, so just set it to 1
dphasevol="dlnM_1",
- parameter_name="M_1",
precode="M_1=math.exp(lnM_1)",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
gridtype="edge"
@@ -290,11 +296,24 @@ def const_dt(self,
range_min = x[0]
range_max = x[1]
dm = x[2]
- m = mmin
- while m <= mmax:
- if m>=range_min and m<=range_max:
- mass_list.append(m)
- m += dm
+ if dm < 0.0:
+ # use log scale
+ dlogm = -dm
+ logm = math.log(mmin)
+ logmmax = math.log(mmax)
+ logrange_min = math.log(range_min)
+ logrange_max = math.log(range_max)
+ while logm <= logmmax:
+ if logm>=logrange_min and logm<=logrange_max:
+ mass_list.append(math.exp(logm))
+ logm += dlogm
+ else:
+ # use linear scale
+ m = mmin
+ while m <= mmax:
+ if m>=range_min and m<=range_max:
+ mass_list.append(m)
+ m += dm
# start time loop at tmax or max_evolution_time
t = tmax if tmax else self.bse_options['max_evolution_time']
@@ -333,7 +352,7 @@ def const_dt(self,
mass_list.remove(0.0)
# apply multiplication factor if given
- if factor:
+ if factor and factor!=1.0:
mass_list = [m * factor for m in mass_list]
# reformat numbers
@@ -352,6 +371,18 @@ def const_dt(self,
print("{i:4d} m={m:13g} t={t:13g} log10(t)={logt:13g}".format(i=i,m=m,t=t,logt=logt))
twas = t
logtwas = logt
+ exit()
# return the mass list as a numpy array
- return np.unique(np.array(mass_list))
+ mass_list = np.unique(np.array(mass_list))
+
+ # perhaps log the masses
+ if logmasses:
+ mass_list = np.log(mass_list)
+ if log10masses:
+ mass_list = np.log10(mass_list)
+
+ if showlist:
+ print("const_dt mass list\n",mass_list)
+
+ return mass_list