fixing the M&S sampling and adding notebook tests (50971270) · Commits · binary_c / binary_c-python

binarycpython/tests/test_notebooks.py

+9 −6

Original line number	Diff line number	Diff line
		@@ -18,16 +18,19 @@ NOTEBOOKS_DIR = os.path.abspath(
		# define notebooks
		notebooks = [
		"notebook_api_functionality.ipynb",
		"notebook_population.ipynb",
		"notebook_individual_systems.ipynb",
		"notebook_custom_logging.ipynb",
		"notebook_extra_features.ipynb",
		"notebook_luminosity_function_single.ipynb",
		"notebook_luminosity_function_binaries.ipynb",
		"notebook_BHBH.ipynb",
		"notebook_common_envelope_evolution.ipynb",
		"notebook_custom_logging.ipynb",
		"notebook_ensembles.ipynb",
		"notebook_event_based_logging.ipynb",
		"notebook_extra_features.ipynb",
		"notebook_HRD.ipynb",
		"notebook_individual_systems.ipynb",
		"notebook_luminosity_function_binaries.ipynb",
		"notebook_luminosity_function_single.ipynb",
		"notebook_population.ipynb",
		"notebook_solar_system.ipynb",
		"notebook_source_file_sampling.ipynb",
		]

binarycpython/utils/population_extensions/Moe_di_Stefano_2017.py

+4 −4

Original line number	Diff line number	Diff line
		@@ -68,8 +68,8 @@ class Moe_di_Stefano_2017:
		"description": "sampler functions to each of the parameters. NEEDS UPDATING",
		},
		"IMF_distribution": {
		"value": "kroupa",
		"description": "IMF choice for the M&S sampling. Currently only supporting 'kroupa': Kroupa 2001 IMF and 'chabrier': Chabrier 2003 IMF.",
		"value": "kroupa2001",
		"description": "IMF choice for the M&S sampling. Currently only supporting 'kroupa2001': Kroupa 2001 IMF and 'chabrier2003': Chabrier 2003 IMF.",
		},
		"ranges": {
		"value": {
		@@ -158,7 +158,7 @@ defaults to [1,1,0,0] i.e. singles and binaries
		""",
		},
		"q_low_extrapolation_method": {
		"value": "linear",
		"value": "flat",
		"description": """
		q extrapolation (below 0.15) method
		none
		@@ -169,7 +169,7 @@ q extrapolation (below 0.15) method
		""",
		},
		"q_high_extrapolation_method": {
		"value": "linear",
		"value": "flat",
		"description": "Same as q_low_extrapolation_method",
		},
		}

binarycpython/utils/population_extensions/distribution_functions.py

+39 −13

Original line number	Diff line number	Diff line
		@@ -1487,6 +1487,12 @@ class distribution_functions:
		q_interpolator, tmp_table, qdata, verbosity=verbosity
		)

		self.vb_debug(
		"\tMoe and di Stefano 2017: build_q_table: using integration_constant_q: {}".format(
		integration_constant_q
		),
		)

		if integration_constant_q > 0:
		# normalise to 1.0 by dividing the data by 1.0/$I
		q_interpolator.multiply_table_column(
		@@ -1498,6 +1504,12 @@ class distribution_functions:
		q_interpolator, tmp_table, qdata, verbosity=verbosity
		)

		self.vb_debug(
		"\tMoe and di Stefano 2017: build_q_table: using new_integration_constant_q: {}".format(
		new_integration_constant_q
		),
		)

		# fail if error in integral > 1e-6 (should be ~ machine precision)
		if abs(1.0 - new_integration_constant_q) > 1e-6:
		self.vb_debug(
		@@ -1505,6 +1517,7 @@ class distribution_functions:
		integration_constant_q
		),
		)

		# set this new table in the cache
		Moecache["rinterpolator_q_given_{}_log10{}".format(m, p)] = q_interpolator
		self.vb_debug(
		@@ -1594,6 +1607,7 @@ class distribution_functions:
		raise ValueError(msg)

		integration_constant_q += x[0] * dq

		return integration_constant_q

		def fill_data(self, sample_values, data_dict):
		@@ -1824,24 +1838,33 @@ class distribution_functions:
		),
		)

		############################################################
		# always require an IMF for the primary star
		#################
		# Calculate the M1_probability
		# we use an IMF distribution function to calculate this probability
		#
		# NB multiply by M1 to convert dN/dM to dN/dlnM
		# (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 power law now.
		# TODO: is this actually the correct way? putting the M1 in there? Do we sample in log space?
		# TODO: when we go to a higher mass, we want to extrapolate the values of M&S
		# Store M_1 value temporarily
		options["M_1_temp"] = options["M_1"]

		# Check value of M_1. If it exceeds the max value then we cap it
		# NOTE: we choose not to extrapolate the data, we just take the largest available value
		# options["M_1"] = np.max([options["M_1"], 30.])

		# Set the M1_probability
		# TODO: the way we choose the IMF here is a bit manual. I want to generalise this a bit more so any distribution can be used
		if options["IMF_distribution"] == "kroupa":
		# Handle choice for IMF
		if options["IMF_distribution"] == "kroupa2001":
		M1_probability = self.Kroupa2001(options["M_1"]) * options["M_1"]
		elif options["IMF_distribution"] == "chabrier":
		elif options["IMF_distribution"] == "chabrier2003":
		M1_probability = self.imf_chabrier2003(options["M_1"]) * options["M_1"]
		else:
		raise ValueError(
		"IMF_distribution choice ({}) for Moe & diStefano 2017 distributions not supported.".format(
		options["IMF_distribution"]
		)
		)

		#
		# Set and return if zero
		prob_dict["M_1"] = M1_probability
		self.vb_debug(
		"\tMoe_di_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob dict ({})".format(
		@@ -1852,9 +1875,8 @@ class distribution_functions:
		# calc_total_probdens(prob_dict)
		# return prob_dict

		"""
		From here we go through the multiplicities.
		"""
		#################
		# From here we go through the multiplicities.
		if multiplicity >= 2:
		# If the multiplicity is higher than 1, we will need to construct the following tables:
		# - period distribution table
		@@ -2321,4 +2343,8 @@ class distribution_functions:
		prob_dict["total_probdens"],
		),
		)

		# restore temp M-1 value
		options["M_1"] = options["M_1_temp"]

		return prob_dict