diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index f5d6db84f0dae9943f62fbed43d3681cd788cafe..cdb941bb4858f25febccdf3720eb401ff18b5a26 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -314,17 +314,16 @@ def three_part_powerlaw(
m0, m1, m2, m_max, p1, p2, p3
)
- #
if m < m0:
- prob = 0 # Below lower bound TODO: make this clear.
- elif m0 < m <= m1:
- prob = three_part_powerlaw_constants[0] * (m ** p1) # Between M0 and M1
- elif m1 < m <= m2:
- prob = three_part_powerlaw_constants[1] * (m ** p2) # Between M1 and M2
- elif m2 < m <= m_max:
- prob = three_part_powerlaw_constants[2] * (m ** p3) # Between M2 and M_MAX
+ prob = 0.0 # Below lower bound
+ elif m <= m1:
+ prob = three_part_powerlaw_constants[0] * (m ** p1) # Between m0 and m1
+ elif m <= m2:
+ prob = three_part_powerlaw_constants[1] * (m ** p2) # Between m1 and m2
+ elif m <= m_max:
+ prob = three_part_powerlaw_constants[2] * (m ** p3) # Between m2 and m_max
else:
- prob = 0 # Above M_MAX
+ prob = 0 # Above m_max
return prob
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 43339aff42c8bfc5494dda276d7ed5110e644f8f..b08016701c1b1e448460650a1be5b5db590ad886 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -50,6 +50,7 @@ import time
import uuid
_count = 0
+from diskcache import Cache
from typing import Union, Any
from collections import (
OrderedDict,
@@ -549,7 +550,7 @@ class Population:
valuerange: Union[list, str],
samplerfunc: str,
probdist: str,
- dphasevol: Union[str, int],
+ dphasevol: Union[str, int] = -1,
gridtype: str = "centred",
branchpoint: int = 0,
branchcode: Union[str, None] = None,
@@ -644,6 +645,11 @@ class Population:
Code added at the very bottom of the block.
"""
+ # check parameters
+ if False and dphasevol != -1.0 and gridtype == 'discrete':
+ print("Error making grid: you have set the phasevol to be not -1 and gridtype to discrete, but a discrete grid has no phasevol calculation. You should only set the gridtype to discrete and not set the phasevol in this case.")
+ sys.exit()
+
# Add grid_variable
grid_variable = {
"name": name,
@@ -666,7 +672,7 @@ class Population:
}
# Check for gridtype input
- if not gridtype in [
+ allowed_gridtypes = [
"edge",
"right",
"right edge",
@@ -675,8 +681,10 @@ class Population:
"centred",
"centre",
"center",
- ]:
- msg = "Unknown gridtype value. Please start another one"
+ "discrete"
+ ]
+ if not gridtype in allowed_gridtypes:
+ msg = "Unknown gridtype {gridtype}. Please choose one of: ".format(gridtype=gridtype) + ','.join(allowed_gridtypes)
raise ValueError(msg)
# Load it into the grid_options
@@ -1367,7 +1375,7 @@ class Population:
# Send closing signal to workers. When they receive this they will terminate
if self.grid_options["verbosity"] >= _LOGGER_VERBOSITY_LEVEL:
- stream_logger.debug(f"Signaling stop to processes") # DEBUG
+ stream_logger.debug(f"Signalling processes to stop") # DEBUG
for _ in range(num_cores):
job_queue.put("STOP")
@@ -2464,7 +2472,10 @@ class Population:
# Setting up the for loop
# Add comment for for loop
self._add_code(
- "# for loop for {}".format(grid_variable["name"]) + "\n",
+ "# for loop for variable {name} gridtype {gridtype}".format(
+ name=grid_variable["name"],
+ gridtype=grid_variable["gridtype"],
+ ) + "\n",
"sampled_values_{} = {}".format(
grid_variable["name"], grid_variable["samplerfunc"]
)
@@ -2472,136 +2483,160 @@ class Population:
if vb:
self._add_code(
- "print('samples','{}',':',sampled_values_{})\n".format(
- grid_variable["name"], grid_variable["name"]
+ "print('samples','{name}',':',sampled_values_{name})\n".format(
+ name=grid_variable["name"],
)
)
- # # Some print statement
- # self._add_code((
- # "print('phasevol_{}:', phasevol_{})".format(grid_variable["name"],
- # grid_variable["name"])
- # + "\n"
- # )
-
- # TODO: make sure this works
- # Adding for loop structure
- # self._add_code((
- # "for {} in sampled_values_{}:".format(
- # grid_variable["name"], grid_variable["name"]
- # )
- # + "\n"
- # )
-
if vb:
self._add_code(
- "print('sample {} from',sampled_values_{})".format(
- grid_variable["name"], grid_variable["name"]
+ "print('sample {name} from',sampled_values_{name})".format(
+ name=grid_variable["name"]
)
+ "\n"
)
-
- # if we're on a centred grid, we have actually
- # one fewer grid point
+ # calculate number of values and starting location
+ #
+ # if we're sampling a continuous variable, we
+ # have one fewer grid point than the length of the
+ # sampled_values list
if (
grid_variable["gridtype"] == "centred"
or grid_variable["gridtype"] == "centre"
or grid_variable["gridtype"] == "center"
+ or grid_variable["gridtype"] == "edge"
+ or grid_variable["gridtype"] == "left edge"
+ or grid_variable["gridtype"] == "left"
+ or grid_variable["gridtype"] == "right"
+ or grid_variable["gridtype"] == "right edge"
):
offset = -1
- else:
+ elif grid_variable["gridtype"] == "discrete":
+ # discrete variables sample all the points
offset = 0
- self._add_code(
- "# for loop for {}".format(grid_variable["name"]) + "\n",
- "nvalues_{} = len(sampled_values_{})+{}".format(
- grid_variable["name"],
- grid_variable["name"],
- offset
- )
- + "\n")
+ start = 0
- # if grid_variable["resolution"] == 0:
- # sample from the sampled_values
+ # for loop over the variable
if vb:
self._add_code(
- "print(\"var {} values \",sampled_values_{},\" len \",len(sampled_values_{})+{},\" gridtype {} offset {}\\n\")\n".format(
- grid_variable["name"],
- grid_variable["name"],
- grid_variable["name"],
- offset,
- grid_variable['gridtype'],
- offset
+ "print(\"var {name} values \",sampled_values_{name},\" len \",len(sampled_values_{name})+{offset},\" gridtype {gridtype} offset {offset}\\n\")\n".format(
+ name=grid_variable["name"],
+ offset=offset,
+ gridtype=grid_variable['gridtype'],
)
)
-
self._add_code(
- "for {}_sample_number in range(max(0,len(sampled_values_{})+{})):".format(
- grid_variable["name"],
- grid_variable["name"],
- offset,
+ "for {name}_sample_number in range({start},len(sampled_values_{name})+{offset}):".format(
+ name=grid_variable["name"],
+ offset=offset,
+ start=start
)
+ "\n"
)
- # old code: use resolution passed in
- # else:
- # # use resolution passed in
- # self._add_code(
- # "for {}_sample_number in range({} if {} != 0 else max(0,len(sampled_values_{})-1)):".format(
- # grid_variable["name"],
- # grid_variable["resolution"],
- # grid_variable["resolution"],
- # grid_variable["name"],
- # )
- # + "\n"
- # )
-
self._increment_indent_depth(+1)
- # {}_this and {}_next are this grid point's index and the next
- # grid point's index, used in many calculations below
- self._add_code("if {}_sample_number == 0:\n".format(grid_variable["name"]))
- self._add_code("{}_this = 0;\n".format(grid_variable["name"]), indent=1)
- self._add_code("else:\n")
- self._add_code(
- "{}_this = {}_sample_number; ".format(
- grid_variable["name"], grid_variable["name"]
- ),
- indent=1,
- )
- self._add_code(
- "\n",
- "{}_next = {}_this + 1".format(
- grid_variable["name"], grid_variable["name"]
- ),
- "\n",
- )
+ # {}_this_index is this grid point's index
+ # {}_prev_index and {}_next_index are the previous and next grid points,
+ # (which can be None if there is no previous or next, or if
+ # previous and next should not be used: this is deliberate)
+ #
- # 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
+ if grid_variable["gridtype"] == "discrete":
+ # discrete grids only care about this,
+ # both prev and next should be None to
+ # force errors where they are used
self._add_code(
- "phasevol_{name} = (sampled_values_{name}[{name}_next]-sampled_values_{name}[{name}_this]) if {name}_next < nvalues_{name} else (sampled_values_{name}[{name}_next-1]-sampled_values_{name}[{name}_this-1])".format(name=grid_variable["name"])
+ "{name}_this_index = {name}_sample_number ".format(
+ name=grid_variable["name"],
+ ),
+ )
+ self._add_code(
+ "\n",
+ "{name}_prev_index = None if {name}_this_index == 0 else ({name}_this_index - 1) ".format(
+ name=grid_variable["name"],
+ ),
+ "\n",
+ )
+ self._add_code(
+ "\n",
+ "{name}_next_index = None if {name}_this_index >= (len(sampled_values_{name})+{offset} - 1) else ({name}_this_index + 1)".format(
+ name=grid_variable["name"],
+ offset=offset
+ ),
+ "\n",
+ )
+
+ elif (grid_variable["gridtype"] == "centred"
+ or grid_variable["gridtype"] == "centre"
+ or grid_variable["gridtype"] == "center"
+ or grid_variable["gridtype"] == "edge"
+ or grid_variable["gridtype"] == "left"
+ or grid_variable["gridtype"] == "left edge"):
+
+ # left and centred grids
+ self._add_code("if {}_sample_number == 0:\n".format(grid_variable["name"]))
+ self._add_code("{}_this_index = 0;\n".format(grid_variable["name"]), indent=1)
+ self._add_code("else:\n")
+ self._add_code("{name}_this_index = {name}_sample_number ".format(name=grid_variable["name"]),indent=1)
+ self._add_code("\n")
+ self._add_code("{name}_prev_index = ({name}_this_index - 1) if {name}_this_index > 0 else None ".format(name=grid_variable["name"]))
+ self._add_code("\n")
+ self._add_code("{name}_next_index = {name}_this_index + 1".format(name=grid_variable["name"]))
+ self._add_code("\n")
+
+ elif(grid_variable["gridtype"] == "right" or
+ grid_variable["gridtype"] == "right edge"):
+
+ # right edged grid
+ self._add_code("if {name}_sample_number == 0:\n".format(name=grid_variable["name"]))
+ self._add_code("{name}_this_index = 1;\n".format(name=grid_variable["name"]),indent=1)
+ self._add_code("else:\n")
+ self._add_code("{name}_this_index = {name}_sample_number + 1 ".format(name=grid_variable["name"],),indent=1)
+ self._add_code("\n")
+ self._add_code("{name}_prev_index = {name}_this_index - 1".format(name=grid_variable["name"]))
+ self._add_code("\n")
+ self._add_code("{name}_next_index = ({name}_this_index + 1) if {name}_this_index < len(sampled_values_{name}) else None".format(name=grid_variable["name"]))
+ self._add_code("\n")
+
+ # calculate phase volume
+ if(grid_variable["dphasevol"] == -1):
+ # no phase volume required so set it to 1.0
+ self._add_code("dphasevol_{name} = 1.0 # 666\n".format(name=grid_variable["name"]))
+
+ elif(grid_variable["gridtype"] == "right" or
+ grid_variable["gridtype"] == "right edge"):
+ # right edges always have this and prev defined
+ self._add_code(
+ "dphasevol_{name} = (sampled_values_{name}[{name}_this_index] - sampled_values_{name}[{name}_prev_index])".format(name=grid_variable["name"])
+ "\n"
)
+ elif grid_variable["gridtype"] == "discrete":
+ # discrete might have next defined, use it if we can,
+ # otherwise use prev
+ self._add_code(
+ "dphasevol_{name} = (sampled_values_{name}[{name}_next_index] - sampled_values_{name}[{name}_this_index]) if {name}_next_index else (sampled_values_{name}[{name}_this_index] - sampled_values_{name}[{name}_prev_index])".format(name=grid_variable["name"])
+ + "\n"
+ )
else:
- self._add_code("phasevol_{} = 1\n".format(grid_variable["name"]))
+ # left and centred always have this and next defined
+ self._add_code(
+ "dphasevol_{name} = (sampled_values_{name}[{name}_next_index] - sampled_values_{name}[{name}_this_index])".format(name=grid_variable["name"])
+ + "\n"
+ )
+
##############
# Add phasevol check:
- self._add_code("if phasevol_{} <= 0:\n".format(grid_variable["name"]))
+ self._add_code("if dphasevol_{name} <= 0:\n".format(name=grid_variable["name"]))
# TODO: We might actually want to add the starcount and probability to the totals regardless.
# n that case we need another local variable which will prevent it from being run but will track those parameters
# Add phasevol check action:
self._add_code(
- 'print("Grid generator: phasevol_{name} <= 0! (this=",{name}_this,"=",sampled_values_{name}[{name}_this],", next=",{name}_next,"=",sampled_values_{name}[{name}_next],") Skipping current sample.")'.format(name=grid_variable["name"])
+ 'print("Grid generator: dphasevol_{name} <= 0! (this=",{name}_this_index,"=",sampled_values_{name}[{name}_this_index],", next=",{name}_next_index,"=",sampled_values_{name}[{name}_next_index],") Skipping current sample.")'.format(name=grid_variable["name"])
+ "\n",
"continue\n",
indent=1,
@@ -2609,34 +2644,22 @@ class Population:
if vb:
self._add_code(
- "print('sample {name} from ',sampled_values_{name},' at this=',{name}_this,', next=',{name}_next)".format(name=grid_variable["name"])
+ "print('sample {name} from ',sampled_values_{name},' at this=',{name}_this_index,', next=',{name}_next_index)".format(name=grid_variable["name"])
+ "\n"
)
- # select point location based on gridtype (left, centre or right)
+ # select sampled 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"
- ):
- self._add_code(
- "{name} = sampled_values_{name}[{name}_this]".format(
- name=grid_variable["name"])
- + "\n"
- )
- if vb:
- self._add_code(
- "print(\"sampled_values_{name}[{name}_this = \",{name}_this,\" {name}_next= \",{name}_next,\"] = \",{name},\"\\\n\")\n".format(
- name=grid_variable["name"],
- ))
- elif (
- grid_variable["gridtype"] == "right"
+ or grid_variable["gridtype"] == "right"
or grid_variable["gridtype"] == "right edge"
+ or grid_variable['gridtype'] == 'discrete'
):
self._add_code(
- + "{name} = sampled_values_{name}[{name}_next]".format(
- name=grid_variable["name"],
- )
+ "{name} = sampled_values_{name}[{name}_this_index]".format(
+ name=grid_variable["name"])
+ "\n"
)
elif (
@@ -2645,17 +2668,17 @@ class Population:
or grid_variable["gridtype"] == "center"
):
self._add_code(
- "{name} = 0.5 * (sampled_values_{name}[{name}_next] + sampled_values_{name}[{name}_this])".format(name=grid_variable["name"])
+ "{name} = 0.5 * (sampled_values_{name}[{name}_next_index] + sampled_values_{name}[{name}_this_index])".format(name=grid_variable["name"])
+ "\n"
)
else:
- msg = "Unknown gridtype value. PLease choose a different one"
+ msg = "Unknown gridtype value {type}.".format(type=grid_variable['gridtype'])
raise ValueError(msg)
if vb:
self._add_code(
- "print('hence {} = ',{})\n".format(
- grid_variable["name"], grid_variable["name"]
+ "print('hence {name} = ',{name})\n".format(
+ name=grid_variable["name"]
)
)
@@ -2666,18 +2689,18 @@ class Population:
if grid_variable["condition"]:
self._add_code(
# Add comment
- "# Condition for {}\n".format(grid_variable["name"]),
+ "# Condition for {name}\n".format(name=grid_variable["name"]),
# Add condition check
- "if not {}:\n".format(grid_variable["condition"]),
+ "if not {condition}:\n".format(condition=grid_variable["condition"]),
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["name"]
+ 'print("Grid generator: Condition for {name} not met!")'.format(
+ name=grid_variable["name"]
)
+ "\n",
"continue" + "\n",
@@ -2691,16 +2714,6 @@ class Population:
# Add some whitespace
self._add_code("\n")
- # # Add debugging:
- # if grid_variable['name']=='q':
- # self._add_code((
- # indent=1,
- # 'print("sampling:", sampled_values_{}, M_1)'.format(
- # grid_variable["name"], grid_variable["name"]
- # )
- # + "\n"
- # )
-
# Add some whitespace
self._add_code("\n")
@@ -2709,8 +2722,8 @@ class Population:
# Add pre-code
if grid_variable["precode"]:
self._add_code(
- "{}".format(
- grid_variable["precode"].replace(
+ "{precode}".format(
+ precode=grid_variable["precode"].replace(
"\n", "\n" + self._indent_block(0)
)
)
@@ -2719,8 +2732,8 @@ class Population:
# Set phasevol
self._add_code(
- "phasevol *= phasevol_{}\n".format(
- grid_variable["name"],
+ "phasevol *= dphasevol_{name}\n".format(
+ name=grid_variable["name"],
)
)
@@ -2730,29 +2743,29 @@ class Population:
self._add_code(
"\n",
"# Setting probabilities\n",
- "d{} = phasevol_{} * ({})".format(
- grid_variable["name"],
- grid_variable["name"],
- grid_variable["probdist"],
+ "d{name} = dphasevol_{name} * ({probdist})".format(
+ name=grid_variable["name"],
+ probdist=grid_variable["probdist"],
)
+ "\n",
# Save probability sum
- "probabilities_sum[{}] += d{}".format(
- grid_variable["grid_variable_number"], grid_variable["name"]
+ "probabilities_sum[{n}] += d{name}".format(
+ n=grid_variable["grid_variable_number"],
+ name=grid_variable["name"]
)
+ "\n",
)
if grid_variable["grid_variable_number"] == 0:
self._add_code(
- "probabilities_list[0] = d{}".format(grid_variable["name"]) + "\n"
+ "probabilities_list[0] = d{name}".format(name=grid_variable["name"]) + "\n"
)
else:
self._add_code(
- "probabilities_list[{}] = probabilities_list[{}] * d{}".format(
- grid_variable["grid_variable_number"],
- grid_variable["grid_variable_number"] - 1,
- grid_variable["name"],
+ "probabilities_list[{this}] = probabilities_list[{prev}] * d{name}".format(
+ this=grid_variable["grid_variable_number"],
+ prev=grid_variable["grid_variable_number"] - 1,
+ name=grid_variable["name"],
)
+ "\n"
)
@@ -2761,8 +2774,8 @@ class Population:
# postcode
if grid_variable["postcode"]:
self._add_code(
- "{}".format(
- grid_variable["postcode"].replace(
+ "{postcode}".format(
+ postcode=grid_variable["postcode"].replace(
"\n", "\n" + self._indent_block(0)
)
)
@@ -2773,14 +2786,14 @@ class Population:
# Increment starcount for this parameter
self._add_code(
"\n",
- "# Increment starcount for {}\n".format(grid_variable["name"]),
- "starcounts[{}] += 1".format(
- grid_variable["grid_variable_number"],
+ "# Increment starcount for {name}\n".format(name=grid_variable["name"]),
+ "starcounts[{n}] += 1".format(
+ n=grid_variable["grid_variable_number"],
)
+ "\n",
# Add value to dict
- 'parameter_dict["{}"] = {}'.format(
- grid_variable["parameter_name"], grid_variable["parameter_name"]
+ 'parameter_dict["{name}"] = {name}'.format(
+ name=grid_variable["parameter_name"]
)
+ "\n",
"\n",
@@ -2826,14 +2839,14 @@ class Population:
self._increment_indent_depth(+1)
self._add_code(
"#" * 40 + "\n",
- "# Code below is for finalising the handling of this iteration of the parameter {}\n".format(
- grid_variable["name"]
+ "# Code below is for finalising the handling of this iteration of the parameter {name}\n".format(
+ name=grid_variable["name"]
),
)
# Set phasevol
# TODO: fix. this isn't supposed to be the value that we give it here. discuss
- self._add_code("phasevol /= phasevol_{}\n\n".format(grid_variable["name"]))
+ self._add_code("phasevol /= dphasevol_{name}\n\n".format(name=grid_variable["name"]))
self._increment_indent_depth(-2)
@@ -2872,7 +2885,7 @@ class Population:
self._increment_indent_depth(+1)
self._add_code("\n", "#" * 40 + "\n", "if print_results:\n")
self._add_code(
- "print('Grid has handled {} stars with a total probability of {:g}'.format(_total_starcount,self.grid_options['_probtot']))\n",
+ "print('Grid has handled {starcount} stars with a total probability of {probtot:g}'.format(starcount=_total_starcount,probtot=self.grid_options['_probtot']))\n",
indent=1,
)
@@ -2895,14 +2908,15 @@ class Population:
# Write to file
gridcode_filename = self._gridcode_filename()
+
self.grid_options["gridcode_filename"] = gridcode_filename
verbose_print(
- "{}Writing grid code to {} [dry_run = {}]{}".format(
- self.ANSI_colours["blue"],
- gridcode_filename,
- dry_run,
- self.ANSI_colours["reset"],
+ "{blue}Writing grid code to {file} [dry_run = {dry}]{reset}".format(
+ blue=self.ANSI_colours["blue"],
+ file=gridcode_filename,
+ dry=dry_run,
+ reset=self.ANSI_colours["reset"],
),
self.grid_options["verbosity"],
1,
@@ -2927,8 +2941,10 @@ class Population:
try:
os.symlink(gridcode_filename, symlink)
verbose_print(
- "{}Symlinked grid code to {} {}".format(
- self.ANSI_colours["blue"], symlink, self.ANSI_colours["reset"]
+ "{blue}Symlinked grid code to {symlink} {reset}".format(
+ blue=self.ANSI_colours["blue"],
+ symlink=symlink,
+ reset=self.ANSI_colours["reset"]
),
self.grid_options["verbosity"],
1,
@@ -2946,12 +2962,12 @@ class Population:
self._add_code("#" * 40 + "\n")
if branchcode:
- self._add_code("# Branch code\nif {}:\n", branchcode)
+ self._add_code("# Branch code\nif {branchcode}:\n".format(branchcode=branchcode))
if branchpoint:
self._add_code(
- "# Code below will get evaluated for every system at this level of multiplicity (last one of that being {})\n".format(
- grid_variable["name"]
+ "# Code below will get evaluated for every system at this level of multiplicity (last one of that being {name})\n".format(
+ name=grid_variable["name"]
)
)
else:
@@ -2963,8 +2979,8 @@ class Population:
self._add_code(
"\n",
"# Weigh the probability by a custom weighting factor\n",
- 'probability = self.grid_options["weight"] * probabilities_list[{}]'.format(
- grid_variable["grid_variable_number"]
+ 'probability = self.grid_options["weight"] * probabilities_list[{n}]'.format(
+ n=grid_variable["grid_variable_number"]
)
+ "\n",
# Take into account the multiplicity fraction:
@@ -2984,8 +3000,8 @@ class Population:
self._add_code(
"_total_starcount += 1\n",
# set probability and phasevol values into the system dict
- 'parameter_dict["{}"] = {}'.format("probability", "probability") + "\n",
- 'parameter_dict["{}"] = {}'.format("phasevol", "phasevol") + "\n",
+ 'parameter_dict["{p}"] = {p}'.format(p="probability") + "\n",
+ 'parameter_dict["{v}"] = {v}'.format(v="phasevol") + "\n",
# Increment total probability
"self._increment_probtot(probability)\n",
indent=1,
@@ -3015,8 +3031,8 @@ class Population:
# Code to load the
verbose_print(
- message="Loading grid code function from {}".format(
- self.grid_options["gridcode_filename"]
+ message="Loading grid code function from {file}".format(
+ file=self.grid_options["gridcode_filename"]
),
verbosity=self.grid_options["verbosity"],
minimal_verbosity=1,
@@ -4165,8 +4181,7 @@ class Population:
condition="({}[int(multiplicity)-1] > 0)".format(
str(self.grid_options["Moe2017_options"]["multiplicity_modulator"])
),
- gridtype="edge",
- dphasevol=-1,
+ gridtype="discrete",
probdist=1,
)
diff --git a/binarycpython/utils/spacing_functions.py b/binarycpython/utils/spacing_functions.py
index 1e4c2ecf72297696b36d304ade77f05baaff531b..b76af0d8b32fe5599e7a86de711ca95b110b7a92 100644
--- a/binarycpython/utils/spacing_functions.py
+++ b/binarycpython/utils/spacing_functions.py
@@ -6,7 +6,10 @@ Tasks:
"""
from typing import Union
+from diskcache import Cache
+import datetime
import functools
+import json
import math
import numpy as np
import py_rinterpolate
@@ -158,25 +161,7 @@ def gaussian_zoom(
@functools.lru_cache(maxsize=16)
-def const_dt(
- self,
- dt=1000.0,
- dlogt=0.1,
- mmin=0.07,
- mmax=100.0,
- nres=1000,
- logspacing=False,
- tmin=3.0, # start at 3Myr
- tmax=None, # use max_evolution_time by default
- mindm=None, # tuple of tuples
- 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,
-):
+def const_dt(self,cachedir=None,usecache=True,**kwargs):
"""
const_dt returns a list of masses spaced at a constant age difference
@@ -197,6 +182,8 @@ def const_dt(
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()
+ usecache: if True (the default) uses cached results if they are saved (in cachedir)
+ cachedir: if None, defaults to grid_options['tmp_dir']/const_dt_cache
Returns:
Array of masses.
@@ -212,222 +199,271 @@ def const_dt(
samplerfunc="const_dt(self,dlogt=0.1,nres=100,mmin=0.07,mmax=80.0,maxdm=((0.07,1.0,0.1),(1.0,10.0,1.0),(10.0,80.0,2.0)),showtable=True,logspacing=True,fsample=1.0/4.0)"
"""
-
- # first, make a stellar lifetime table
- #
- # we should use the bse_options from self
- # so our lifetime_population uses the same physics
- lifetime_population = Population()
- lifetime_population.bse_options = dict(self.bse_options)
-
- # we only want to evolve the star during nuclear burning,
- # we don't want a dry run of the grid
- # we want to use the right number of CPU cores
- lifetime_population.set(
- do_dry_run=False,
- num_cores=self.grid_options["num_cores"],
- max_stellar_type_1=10,
- save_ensemble_chunks=False,
- symlink_latest_gridcode=False,
- modulo=1,
- start_at=0,
- slurm=0,
- condor=0,
- )
-
- # 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)],
- samplerfunc="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",
- precode="M_1=math.exp(lnM_1)",
- condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
- gridtype="edge",
- )
-
- # set up the parse function
- def _parse_function(self, output):
- if output:
- for line in output.splitlines():
- data = line.split()
- if data[0] == "SINGLE_STAR_LIFETIME":
- # append (log10(mass), log10(lifetime)) tuples
- logm = math.log10(float(data[1]))
- logt = math.log10(float(data[2]))
- # print(line)
- # print("logM=",logm,"M=",10.0**logm," -> logt=",logt)
- self.grid_results["interpolation table m->t"][logm] = logt
- self.grid_results["interpolation table t->m"][logt] = logm
-
- lifetime_population.set(
- parse_function=_parse_function,
- )
-
- # run to build the interpolation table
- print("Running population to make lifetime interpolation table, please wait")
- lifetime_population.evolve()
- # print("Data table",lifetime_population.grid_results['interpolation table t->m'])
-
- if not "interpolation table t->m" in lifetime_population.grid_results or len(lifetime_population.grid_results["interpolation table t->m"].keys()) == 0:
- print("\n\n\nError: The t->m lifetime table is empty. One usual cause for this is that the tmax or max_evolution_time option (currently passed in to const_dt as {tmax}) is too short for there to be any entries in the table before the first timestep. Try increasing tmax and max_evolution_time, shorten the timestep or, if using log times, set tstart to be closer to 0.\n".format(tmax=tmax))
- exit()
-
-
- # convert to nested lists for the interpolator
- #
- # make time -> mass table
- data_table_time_mass = []
- times = sorted(lifetime_population.grid_results["interpolation table t->m"].keys())
- for time in times:
- mass = lifetime_population.grid_results["interpolation table t->m"][time]
- # we have to make sure the time is monotonic (not guaranteed at high mass)
- if len(data_table_time_mass) == 0:
- data_table_time_mass.append([time, mass])
- elif mass < data_table_time_mass[-1][1]:
- data_table_time_mass.append([time, mass])
-
- # make mass -> time table
- data_table_mass_time = []
- masses = sorted(lifetime_population.grid_results["interpolation table m->t"].keys())
- for mass in masses:
- time = lifetime_population.grid_results["interpolation table m->t"][mass]
- data_table_mass_time.append([mass, time])
-
- # set up interpolators
- interpolator_time_mass = py_rinterpolate.Rinterpolate(
- table=data_table_time_mass, nparams=1, ndata=1, verbosity=0 # mass # lifetime
- )
- interpolator_mass_time = py_rinterpolate.Rinterpolate(
- table=data_table_mass_time, nparams=1, ndata=1, verbosity=0 # lifetime # mass
- )
-
- # function to get a mass given a time (Myr)
- def _mass_from_time(linear_time):
- return 10.0 ** interpolator_time_mass.interpolate([math.log10(linear_time)])[0]
-
- # function to get a time given a mass (Msun)
- def _time_from_mass(mass):
- return 10.0 ** interpolator_mass_time.interpolate([math.log10(mass)])[0]
-
- # return a unique list
- def _uniq(_list):
- return sorted(list(set(_list)))
-
- # format a whole list like %g
- def _format(_list):
- return [float("{x:g}".format(x=x)) for x in _list]
-
- # construct mass list, always include the min and max
- mass_list = [mmin, mmax]
-
- # first, make sure the stars are separated by only
- # maxdm
- if maxdm:
- for x in maxdm:
- range_min = x[0]
- range_max = x[1]
- dm = x[2]
- 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"]
-
- # set default mass list
- if logspacing:
- logt = math.log10(t)
- logtmin = math.log10(tmin)
- while logt > logtmin:
- m = _mass_from_time(10.0 ** logt)
- mass_list.append(m)
- logt = max(logtmin, logt - dlogt * fsample)
- else:
- while t > tmin:
- m = _mass_from_time(t)
- mass_list.append(m)
- t = max(tmin, t - dt * fsample)
-
- # make mass list unique
- mass_list = _uniq(mass_list)
-
- if mindm:
- for x in mindm:
- range_min = x[0]
- range_max = x[1]
- mindm = x[2]
- # impose a minimum dm: if two masses in the list
- # are separated by < this, remove the second
- for index, mass in enumerate(mass_list):
- if index > 0 and mass >= range_min and mass <= range_max:
- dm = mass_list[index] - mass_list[index - 1]
- if dm < mindm:
- mass_list[index - 1] = 0.0
- mass_list = _uniq(mass_list)
- if mass_list[0] == 0.0:
- mass_list.remove(0.0)
-
- # apply multiplication factor if given
- if factor and factor != 1.0:
- mass_list = [m * factor for m in mass_list]
-
- # reformat numbers
- mass_list = _format(mass_list)
-
- # show the mass<>time table?
- if showtable:
- twas = 0.0
- logtwas = 0.0
- for i, m in enumerate(mass_list):
- t = _time_from_mass(m)
+ if cachedir == None:
+ cachedir = self.grid_options['tmp_dir'] + '/const_dt_cache'
+ cache = Cache(cachedir)
+ @cache.memoize()
+ def _const_dt(cachedir,
+ num_cores,
+ bse_options,
+ dt=1000.0,
+ dlogt=0.1,
+ mmin=0.07,
+ mmax=100.0,
+ nres=1000,
+ logspacing=False,
+ tmin=3.0, # start at 3Myr
+ tmax=None, # use max_evolution_time by default
+ mindm=None, # tuple of tuples
+ 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,
+ usecache=True,
+ ):
+
+ # first thing to do is make a stellar lifetime table
+ #
+ # we should use the bse_options passed in
+ # so our lifetime_population uses the same physics
+ # as the main grid
+
+ # convert bse_options to dict
+ bse_options = json.loads(bse_options)
+
+ # make new population object
+ lifetime_population = Population()
+ lifetime_population.bse_options = bse_options
+
+ # we only want to evolve the star during nuclear burning,
+ # we don't want a dry run of the grid
+ # we want to use the right number of CPU cores
+ lifetime_population.set(
+ do_dry_run=False,
+ num_cores=num_cores,
+ max_stellar_type_1=10,
+ save_ensemble_chunks=False,
+ symlink_latest_gridcode=False,
+ modulo=1,
+ start_at=0,
+ slurm=0,
+ condor=0,
+ )
+
+ # 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)],
+ samplerfunc="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",
+ precode="M_1=math.exp(lnM_1)",
+ condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
+ gridtype="edge",
+ )
+
+ # set up the parse function
+ def _parse_function(self, output):
+ if output:
+ for line in output.splitlines():
+ data = line.split()
+ if data[0] == "SINGLE_STAR_LIFETIME":
+ # append (log10(mass), log10(lifetime)) tuples
+ logm = math.log10(float(data[1]))
+ logt = math.log10(float(data[2]))
+ # print(line)
+ # print("logM=",logm,"M=",10.0**logm," -> logt=",logt)
+ self.grid_results["interpolation table m->t"][logm] = logt
+ self.grid_results["interpolation table t->m"][logt] = logm
+
+ lifetime_population.set(
+ parse_function=_parse_function,
+ )
+
+ # run to build the interpolation table
+ print("Running population to make lifetime interpolation table, please wait")
+ lifetime_population.evolve()
+ # print("Data table",lifetime_population.grid_results['interpolation table t->m'])
+
+ if not "interpolation table t->m" in lifetime_population.grid_results or len(lifetime_population.grid_results["interpolation table t->m"].keys()) == 0:
+ print("\n\n\nError: The t->m lifetime table is empty. One usual cause for this is that the tmax or max_evolution_time option (currently passed in to const_dt as {tmax}) is too short for there to be any entries in the table before the first timestep. Try increasing tmax and max_evolution_time, shorten the timestep or, if using log times, set tstart to be closer to 0.\n".format(tmax=tmax))
+ exit()
+
+
+ # convert to nested lists for the interpolator
+ #
+ # make time -> mass table
+ data_table_time_mass = []
+ times = sorted(lifetime_population.grid_results["interpolation table t->m"].keys())
+ for time in times:
+ mass = lifetime_population.grid_results["interpolation table t->m"][time]
+ # we have to make sure the time is monotonic (not guaranteed at high mass)
+ if len(data_table_time_mass) == 0:
+ data_table_time_mass.append([time, mass])
+ elif mass < data_table_time_mass[-1][1]:
+ data_table_time_mass.append([time, mass])
+
+ # make mass -> time table
+ data_table_mass_time = []
+ masses = sorted(lifetime_population.grid_results["interpolation table m->t"].keys())
+ for mass in masses:
+ time = lifetime_population.grid_results["interpolation table m->t"][mass]
+ data_table_mass_time.append([mass, time])
+
+ # set up interpolators
+ interpolator_time_mass = py_rinterpolate.Rinterpolate(
+ table=data_table_time_mass, nparams=1, ndata=1, verbosity=0 # mass # lifetime
+ )
+ interpolator_mass_time = py_rinterpolate.Rinterpolate(
+ table=data_table_mass_time, nparams=1, ndata=1, verbosity=0 # lifetime # mass
+ )
+
+ # function to get a mass given a time (Myr)
+ def _mass_from_time(linear_time):
+ return 10.0 ** interpolator_time_mass.interpolate([math.log10(linear_time)])[0]
+
+ # function to get a time given a mass (Msun)
+ def _time_from_mass(mass):
+ return 10.0 ** interpolator_mass_time.interpolate([math.log10(mass)])[0]
+
+ # return a unique list
+ def _uniq(_list):
+ return sorted(list(set(_list)))
+
+ # format a whole list like %g
+ def _format(_list):
+ return [float("{x:g}".format(x=x)) for x in _list]
+
+ # construct mass list, always include the min and max
+ mass_list = [mmin, mmax]
+
+ # first, make sure the stars are separated by only
+ # maxdm
+ if maxdm:
+ for x in maxdm:
+ range_min = x[0]
+ range_max = x[1]
+ dm = x[2]
+ 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 bse_options['max_evolution_time']
+
+ # set default mass list
+ if logspacing:
logt = math.log10(t)
- if twas > 0.0:
- print(
- "{i:4d} m={m:13g} t={t:13g} log10(t)={logt:13g} dt={dt:13g} dlog10(t)={dlogt:13g}".format(
- i=i, m=m, t=t, logt=logt, dt=twas - t, dlogt=logtwas - logt
+ logtmin = math.log10(tmin)
+ while logt > logtmin:
+ m = _mass_from_time(10.0 ** logt)
+ mass_list.append(m)
+ logt = max(logtmin, logt - dlogt * fsample)
+ else:
+ while t > tmin:
+ m = _mass_from_time(t)
+ mass_list.append(m)
+ t = max(tmin, t - dt * fsample)
+
+ # make mass list unique
+ mass_list = _uniq(mass_list)
+
+ if mindm:
+ for x in mindm:
+ range_min = x[0]
+ range_max = x[1]
+ mindm = x[2]
+ # impose a minimum dm: if two masses in the list
+ # are separated by < this, remove the second
+ for index, mass in enumerate(mass_list):
+ if index > 0 and mass >= range_min and mass <= range_max:
+ dm = mass_list[index] - mass_list[index - 1]
+ if dm < mindm:
+ mass_list[index - 1] = 0.0
+ mass_list = _uniq(mass_list)
+ if mass_list[0] == 0.0:
+ mass_list.remove(0.0)
+
+ # apply multiplication factor if given
+ if factor and factor != 1.0:
+ mass_list = [m * factor for m in mass_list]
+
+ # reformat numbers
+ mass_list = _format(mass_list)
+
+ # show the mass<>time table?
+ if showtable:
+ twas = 0.0
+ logtwas = 0.0
+ for i, m in enumerate(mass_list):
+ t = _time_from_mass(m)
+ logt = math.log10(t)
+ if twas > 0.0:
+ print(
+ "{i:4d} m={m:13g} t={t:13g} log10(t)={logt:13g} dt={dt:13g} dlog10(t)={dlogt:13g}".format(
+ i=i, m=m, t=t, logt=logt, dt=twas - t, dlogt=logtwas - logt
+ )
)
- )
- else:
- print(
- "{i:4d} m={m:13g} t={t:13g} log10(t)={logt:13g}".format(
- i=i, m=m, t=t, logt=logt
+ else:
+ 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()
+ twas = t
+ logtwas = logt
+ exit()
- # return the mass list as a numpy array
- mass_list = np.unique(np.array(mass_list))
+ # return the mass list as a numpy array
+ 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)
+ # perhaps log the masses
+ if logmasses:
+ mass_list = np.log(mass_list)
+ if log10masses:
+ mass_list = np.log10(mass_list)
- if showlist:
+
+
+ return mass_list
+
+ # call _const_dt and return the mass_list
+ #
+ # Note: because _const_dt is cached to disk, calling it may
+ # use the cached result.
+ #
+ # Note: we send a sorted JSON string instead of the
+ # bse_options dict to make sure the order is preserved
+ mass_list = _const_dt(cachedir,
+ self.grid_options['num_cores'],
+ json.dumps(self.bse_options,sort_keys=True),
+ **kwargs
+ )
+ cache.close()
+
+ if kwargs.get('showlist',True):
print("const_dt mass list ({} masses)\n".format(len(mass_list)), mass_list)
return mass_list
diff --git a/setup.py b/setup.py
index 721dd841d81ef5ec0906bac3b2075e4a4602cb02..077adef0dbf84e2eec2a6aade9ad2a8e7046a6fd 100644
--- a/setup.py
+++ b/setup.py
@@ -265,6 +265,7 @@ setup(
"colorama",
"compress_pickle",
"datasize",
+ "diskcache",
"h5py",
"halo",
"humanize",