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 ea06fc4f567a61e6c46c26e170f70364185151a4..4fc3e4a0ac95be3784158d36fd573749ff3fa5e9 100644 --- a/binarycpython/utils/grid.py +++ b/binarycpython/utils/grid.py @@ -494,7 +494,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, @@ -589,6 +589,11 @@ class Population: Code added at the very bottom of the block. """ + # check parameters + if 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, @@ -611,7 +616,7 @@ class Population: } # Check for gridtype input - if not gridtype in [ + allowed_gridtypes = [ "edge", "right", "right edge", @@ -620,8 +625,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 @@ -1238,7 +1245,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") @@ -2349,144 +2356,170 @@ 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"] ) + "\n") - if vb: + if True: self._add_code( - "print('samples','{}',':',sampled_values_{})\n".format( - grid_variable["name"], grid_variable["name"] + "print('samples','{name}',':',np.exp(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: - offset = 0 + if (grid_variable["gridtype"] == "right" + or grid_variable["gridtype"] == "right edge"): + start = 0 + else: + start = 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") + elif grid_variable["gridtype"] == "discrete": + # discrete variables sample all the points + offset = 0 + 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( + "{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( - "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"]) + "\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["gridtype"] == "discrete" or grid_variable["dphasevol"] == -1): + # no phase volume required for discrete calculations + # 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" ) else: - self._add_code("phasevol_{} = 1\n".format(grid_variable["name"])) + # left or 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, @@ -2494,34 +2527,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 ( @@ -2530,17 +2551,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"] ) ) @@ -2551,18 +2572,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", @@ -2576,16 +2597,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") @@ -2594,8 +2605,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) ) ) @@ -2604,8 +2615,8 @@ class Population: # Set phasevol self._add_code( - "phasevol *= phasevol_{}\n".format( - grid_variable["name"], + "phasevol *= dphasevol_{name}\n".format( + name=grid_variable["name"], ) ) @@ -2615,29 +2626,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" ) @@ -2646,8 +2657,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) ) ) @@ -2658,14 +2669,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", @@ -2711,14 +2722,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) @@ -2757,7 +2768,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, ) @@ -2781,16 +2792,16 @@ class Population: # Write to file gridcode_filename = os.path.join( self.grid_options["tmp_dir"], - "binary_c_grid_{}.py".format(self.grid_options["_population_id"]), + "binary_c_grid_{id}.py".format(id=self.grid_options["_population_id"]), ) 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, @@ -2812,8 +2823,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, @@ -2831,12 +2844,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: @@ -2848,8 +2861,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: @@ -2869,8 +2882,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, @@ -2900,8 +2913,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, @@ -4270,8 +4283,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, )