Source code for binarycpython.utils.distribution_functions

"""
Module containing the predefined distribution functions

The user can use any of these distribution functions to
generate probability distributions for sampling populations

There are distributions for the following parameters:
    - mass
    - period
    - mass ratio
    - binary fraction

Tasks:
    - TODO: make some things globally present? rob does this in his module..i guess it saves calculations but not sure if im gonna do that now
    - TODO: make global constants stuff
    - TODO: add eccentricity distribution: thermal
    - TODO: Add SFH distributions depending on redshift
    - TODO: Add metallicity distributions depending on redshift
    - TODO: Add initial rotational velocity distributions
"""

import math
import numpy as np
from typing import Optional, Union
from binarycpython.utils.useful_funcs import calc_period_from_sep

###
# File containing probability distributions
# Mostly copied from the perl modules
LOG_LN_CONVERTER = 1.0 / math.log(10.0)
distribution_constants = {}  # To store the constants in


[docs]def prepare_dict(global_dict: dict, list_of_sub_keys: list) -> None: """ Function that makes sure that the global dict is prepared to have a value set there. This dictionary will store values and factors for the distribution functions, so that they dont have to be calculated each time. Args: global_dict: globablly acessible dictionary where factors are stored in list_of_sub_keys: List of keys that must become be(come) present in the global_dict """ internal_dict_value = global_dict # This loop almost mimics a recursive loop into the dictionary. # It checks whether the first key of the list is present, if not; set it with an empty dict. # Then it overrides itself to be that (new) item, and goes on to do that again, until the list exhausted for k in list_of_sub_keys: # If the sub key doesnt exist then make an empty dict if not internal_dict_value.get(k, None): internal_dict_value[k] = {} internal_dict_value = internal_dict_value[k]
[docs]def set_opts(opts: dict, newopts: dict) -> dict: """ Function to take a default dict and override it with newer values. # TODO: consider changing this to just a dict.update Args: opts: dictionary with default values newopts: dictionary with new values Returns: returns an updated dictionary """ if newopts: for opt in newopts.keys(): if opt in opts.keys(): opts[opt] = newopts[opt] return opts
[docs]def flat() -> float: """ Dummy distribution function that returns 1 Returns: a flat uniform distribution: 1 """ return 1.0
[docs]def number(value: Union[int, float]) -> Union[int, float]: """ Dummy distribution function that returns the input Args: value: the value that will be returned by this function. Returns: the value that was provided """ return value
[docs]def const( min_bound: Union[int, float], max_bound: Union[int, float], val: float = None ) -> Union[int, float]: """ a constant distribution function between min=min_bound and max=max_bound. Args: min_bound: lower bound of the range max_bound: upper bound of the range Returns: returns the value of 1/(max_bound-min_bound). If val is provided, it will check whether min_bound < val <= max_bound. if not: returns 0 """ if val: if not min_bound < val <= max_bound: print("out of bounds") prob = 0 return prob prob = 1.0 / (max_bound - min_bound) return prob
[docs]def powerlaw_constant( min_val: Union[int, float], max_val: Union[int, float], k: Union[int, float] ) -> Union[int, float]: """ Function that returns the constant to normalise a powerlaw TODO: what if k is -1? Args: min_val: lower bound of the range max_val: upper bound of the range k: powerlaw slope Returns: constant to normalize the given powerlaw between the min_val and max_val range """ k1 = k + 1.0 # print( # "Powerlaw consts from {} to {}, k={} where k1={}".format( # min_val, max_val, k, k1 # ) # ) powerlaw_const = k1 / (max_val ** k1 - min_val ** k1) return powerlaw_const
[docs]def powerlaw( min_val: Union[int, float], max_val: Union[int, float], k: Union[int, float], x: Union[int, float], ) -> Union[int, float]: """ Single powerlaw with index k at x from min to max Args: min_val: lower bound of the powerlaw max_val: upper bound of the powerlaw k: slope of the power law x: position at which we want to evaluate Returns: `probability` at the given position(x) """ # Handle faulty value if k == -1: print("wrong value for k") raise ValueError if (x < min_val) or (x > max_val): print("input value is out of bounds!") return 0 powerlaw_const = powerlaw_constant(min_val, max_val, k) # powerlaw prob = powerlaw_const * (x ** k) # print( # "Power law from {} to {}: const = {}, y = {}".format( # min_val, max_val, const, y # ) # ) return prob
[docs]def calculate_constants_three_part_powerlaw( m0: Union[int, float], m1: Union[int, float], m2: Union[int, float], m_max: Union[int, float], p1: Union[int, float], p2: Union[int, float], p3: Union[int, float], ) -> Union[int, float]: """ Function to calculate the constants for a three-part powerlaw TODO: use the powerlaw_constant function to calculate all these values Args: m0: lower bound mass m1: second boundary, between the first slope and the second slope m2: third boundary, between the second slope and the third slope m_max: upper bound mass p1: first slope p2: second slope p3: third slope Returns: array of normalisation constants """ # print("Initialising constants for the three-part powerlaw: m0={} m1={} m2={}\ # m_max={} p1={} p2={} p3={}\n".format(m0, m1, m2, m_max, p1, p2, p3)) array_constants_three_part_powerlaw = [0, 0, 0] array_constants_three_part_powerlaw[1] = ( ((m1 ** p2) * (m1 ** (-p1))) * (1.0 / (1.0 + p1)) * (m1 ** (1.0 + p1) - m0 ** (1.0 + p1)) ) array_constants_three_part_powerlaw[1] += ( (m2 ** (1.0 + p2) - m1 ** (1.0 + p2)) ) * (1.0 / (1.0 + p2)) array_constants_three_part_powerlaw[1] += ( ((m2 ** p2) * (m2 ** (-p3))) * (1.0 / (1.0 + p3)) * (m_max ** (1.0 + p3) - m2 ** (1.0 + p3)) ) array_constants_three_part_powerlaw[1] = 1.0 / ( array_constants_three_part_powerlaw[1] + 1e-50 ) array_constants_three_part_powerlaw[0] = array_constants_three_part_powerlaw[1] * ( (m1 ** p2) * (m1 ** (-p1)) ) array_constants_three_part_powerlaw[2] = array_constants_three_part_powerlaw[1] * ( (m2 ** p2) * (m2 ** (-p3)) ) return array_constants_three_part_powerlaw
# $$array[1]=(($m1**$p2)*($m1**(-$p1)))* # (1.0/(1.0+$p1))* # ($m1**(1.0+$p1)-$m0**(1.0+$p1))+ # (($m2**(1.0+$p2)-$m1**(1.0+$p2)))* # (1.0/(1.0+$p2))+ # (($m2**$p2)*($m2**(-$p3)))* # (1.0/(1.0+$p3))* # ($mmax**(1.0+$p3)-$m2**(1.0+$p3)); # $$array[1]=1.0/($$array[1]+1e-50); # $$array[0]=$$array[1]*$m1**$p2*$m1**(-$p1); # $$array[2]=$$array[1]*$m2**$p2*$m2**(-$p3); # #print "ARRAY SET @_ => @$array\n"; # $threepart_powerlaw_consts{"@_"}=[@$array];
[docs]def three_part_powerlaw( m: Union[int, float], m0: Union[int, float], m1: Union[int, float], m2: Union[int, float], m_max: Union[int, float], p1: Union[int, float], p2: Union[int, float], p3: Union[int, float], ) -> Union[int, float]: """ Generalized three-part power law, usually used for mass distributions Args: m: mass at which we want to evaluate the distribution. m0: lower bound mass m1: second boundary, between the first slope and the second slope m2: third boundary, between the second slope and the third slope m_max: upper bound mass p1: first slope p2: second slope p3: third slope Returns: 'probability' at given mass m """ # TODO: add check on whether the values exist three_part_powerlaw_constants = calculate_constants_three_part_powerlaw( m0, m1, m2, m_max, p1, p2, p3 ) # if m < m0: prob = 0 # Below lower bound 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 else: prob = 0 # Above M_MAX return prob
[docs]def gaussian_normalizing_const( mean: Union[int, float], sigma: Union[int, float], gmin: Union[int, float], gmax: Union[int, float], ) -> Union[int, float]: """ Function to calculate the normalisation constant for the gaussian Args: mean: mean of the gaussian sigma: standard deviation of the gaussian gmin: lower bound of the range to calculate the probabilities in gmax: upper bound of the range to calculate the probabilities in Returns: normalisation constant for the gaussian distribution(mean, sigma) between gmin and gmax """ # First time; calculate multipllier for given mean and sigma ptot = 0 resolution = 1000 d = (gmax - gmin) / resolution for i in range(resolution): y = gmin + i * d ptot += d * gaussian_func(y, mean, sigma) # TODO: Set value in global return ptot
[docs]def gaussian_func( x: Union[int, float], mean: Union[int, float], sigma: Union[int, float] ) -> Union[int, float]: """ Function to evaluate a gaussian at a given point, but this time without any boundaries. Args: x: location at which to evaluate the distribution mean: mean of the gaussian sigma: standard deviation of the gaussian Returns: value of the gaussian at x """ gaussian_prefactor = 1.0 / math.sqrt(2.0 * math.pi) r = 1.0 / (sigma) y = (x - mean) * r return gaussian_prefactor * r * math.exp(-0.5 * y ** 2)
[docs]def gaussian( x: Union[int, float], mean: Union[int, float], sigma: Union[int, float], gmin: Union[int, float], gmax: Union[int, float], ) -> Union[int, float]: """ Gaussian distribution function. used for e..g Duquennoy + Mayor 1991 Args: x: location at which to evaluate the distribution mean: mean of the gaussian sigma: standard deviation of the gaussian gmin: lower bound of the range to calculate the probabilities in gmax: upper bound of the range to calculate the probabilities in Returns: 'probability' of the gaussian distribution between the boundaries, evaluated at x """ # # location (X value), mean and sigma, min and max range # my ($x,$mean,$sigma,$gmin,$gmax) = @_; if (x < gmin) or (x > gmax): prob = 0 else: # normalize over given range # TODO: add loading into global var normalisation = gaussian_normalizing_const(mean, sigma, gmin, gmax) prob = normalisation * gaussian_func(x, mean, sigma) return prob
##### # Mass distributions #####
[docs]def Kroupa2001(m: Union[int, float], newopts: dict = None) -> Union[int, float]: """ Probability distribution function for kroupa 2001 IMF, where the default values to the three_part_powerlaw are: default = {"m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 100, "p1": -1.3, "p2": -2.3,"p3": -2.3} Args: m: mass to evaluate the distribution at newopts: optional dict to override the default values. Returns: 'probability' of distribution function evaluated at m """ # Default params and override them default = { "m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 100, "p1": -1.3, "p2": -2.3, "p3": -2.3, } value_dict = default.copy() if newopts: value_dict.update(newopts) return three_part_powerlaw( m, value_dict["m0"], value_dict["m1"], value_dict["m2"], value_dict["mmax"], value_dict["p1"], value_dict["p2"], value_dict["p3"], )
[docs]def ktg93(m: Union[int, float], newopts: dict = None) -> Union[int, float]: """ Probability distribution function for KTG93 IMF, where the default values to the three_part_powerlaw are: default = {"m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 80, "p1": -1.3, "p2": -2.2,"p3": -2.7} Args: m: mass to evaluate the distribution at newopts: optional dict to override the default values. Returns: 'probability' of distribution function evaluated at m """ # TODO: ask rob what this means # if($m eq 'uncertainties') # { # # return (pointer to) the uncertainties hash # return { # m0=>{default=>0.1, # fixed=>1}, # m1=>{default=>0.5, # fixed=>1}, # m2=>{default=>1.0, # fixed=>1}, # mmax=>{default=>80.0, # fixed=>1}, # p1=>{default=>-1.3, # low=>-1.3, # high=>-1.3}, # p2=>{default=>-2.2, # low=>-2.2, # high=>-2.2}, # p3=>{default=>-2.7, # low=>-2.7, # high=>-2.7} # }; # } # set options # opts = set_opts({'m0':0.1, 'm1':0.5, 'm2':1.0, 'mmax':80, 'p1':-1.3, 'p2':-2.2, 'p3':-2.7}, # newopts) defaults = { "m0": 0.1, "m1": 0.5, "m2": 1.0, "mmax": 80, "p1": -1.3, "p2": -2.2, "p3": -2.7, } value_dict = defaults.copy() if newopts: value_dict.update(newopts) return three_part_powerlaw( m, value_dict["m0"], value_dict["m1"], value_dict["m2"], value_dict["mmax"], value_dict["p1"], value_dict["p2"], value_dict["p3"], )
# sub ktg93_lnspace # { # # wrapper for KTG93 on a ln(m) grid # my $m=$_[0]; # return ktg93(@_) * $m; # }
[docs]def imf_tinsley1980(m: Union[int, float]) -> Union[int, float]: """ Probability distribution function for tinsley 1980 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3) Args: m: mass to evaluate the distribution at Returns: 'probability' of distribution function evaluated at m """ return three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3)
[docs]def imf_scalo1986(m: Union[int, float]) -> Union[int, float]: """ Probability distribution function for Scalo 1986 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70) Args: m: mass to evaluate the distribution at Returns: 'probability' of distribution function evaluated at m """ return three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70)
[docs]def imf_scalo1998(m: Union[int, float]) -> Union[int, float]: """ From scalo 1998 Probability distribution function for Scalo 1998 IMF (defined up until 80Msol): three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3) Args: m: mass to evaluate the distribution at Returns: 'probability' of distribution function evaluated at m """ return three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3)
[docs]def imf_chabrier2003(m: Union[int, float]) -> Union[int, float]: """ Probability distribution function for IMF of Chabrier 2003 PASP 115:763-795 Args: m: mass to evaluate the distribution at Returns: 'probability' of distribution function evaluated at m """ chabrier_logmc = math.log10(0.079) chabrier_sigma2 = 0.69 * 0.69 chabrier_a1 = 0.158 chabrier_a2 = 4.43e-2 chabrier_x = -1.3 if m <= 0: print("below bounds") raise ValueError if 0 < m < 1.0: A = 0.158 dm = math.log10(m) - chabrier_logmc prob = chabrier_a1 * math.exp(-(dm ** 2) / (2.0 * chabrier_sigma2)) else: prob = chabrier_a2 * (m ** chabrier_x) prob = prob / (0.1202462 * m * math.log(10)) return prob
######################################################################## # Binary fractions ########################################################################
[docs]def Arenou2010_binary_fraction(m: Union[int, float]) -> Union[int, float]: """ Arenou 2010 function for the binary fraction as f(M1) GAIA-C2-SP-OPM-FA-054 www.rssd.esa.int/doc_fetch.php?id=2969346 Args: m: mass to evaluate the distribution at Returns: binary fraction at m """ return 0.8388 * math.tanh(0.688 * m + 0.079)
# print(Arenou2010_binary_fraction(0.4))
[docs]def raghavan2010_binary_fraction(m: Union[int, float]) -> Union[int, float]: """ Fit to the Raghavan 2010 binary fraction as a function of spectral type (Fig 12). Valid for local stars (Z=Zsolar). The spectral type is converted mass by use of the ZAMS effective temperatures from binary_c/BSE (at Z=0.02) and the new "long_spectral_type" function of binary_c (based on Jaschek+Jaschek's Teff-spectral type table). Rob then fitted the result Args: m: mass to evaluate the distribution at Returns: binary fraction at m """ return min( 1.0, max( (m ** 0.1) * (5.12310e-01) + (-1.02070e-01), (1.10450e00) * (m ** (4.93670e-01)) + (-6.95630e-01), ), )
# print(raghavan2010_binary_fraction(2)) ######################################################################## # Period distributions ########################################################################
[docs]def duquennoy1991(logper: Union[int, float]) -> Union[int, float]: """ Period distribution from Duquennoy + Mayor 1991. Evaluated the function gaussian(logper, 4.8, 2.3, -2, 12) Args: logper: logarithm of period to evaluate the distribution at Returns: 'probability' at gaussian(logper, 4.8, 2.3, -2, 12) """ return gaussian(logper, 4.8, 2.3, -2, 12)
[docs]def sana12( M1: Union[int, float], M2: Union[int, float], a: Union[int, float], P: Union[int, float], amin: Union[int, float], amax: Union[int, float], x0: Union[int, float], x1: Union[int, float], p: Union[int, float], ) -> Union[int, float]: """ distribution of initial orbital periods as found by Sana et al. (2012) which is a flat distribution in ln(a) and ln(P) respectively for stars * less massive than 15Msun (no O-stars) * mass ratio q=M2/M1<0.1 * log(P)<0.15=x0 and log(P)>3.5=x1 and is be given by dp/dlogP ~ (logP)^p for all other binary configurations (default p=-0.55) arguments are M1, M2, a, Period P, amin, amax, x0=log P0, x1=log P1, p example args: 10, 5, sep(M1, M2, P), sep, ?, -2, 12, -0.55 # TODO: Fix this function! Args: M1: Mass of primary M2: Mass of secondary a: separation of binary P: period of binary amin: minimum separation of the distribution (lower bound of the range) amax: maximum separation of the distribution (upper bound of the range) x0: log of minimum period of the distribution (lower bound of the range) x1: log of maximum period of the distribution (upper bound of the range) p: slope of the distributoon Returns: 'probability' of orbital period P given the other parameters """ res = 0 if (M1 < 15) or (M2 / M1 < 0.1): res = 1.0 / (math.log(amax) - math.log(amin)) else: p1 = 1.0 + p # For more details see the LyX document of binary_c for this distribution # where the variables and normalizations are given # we use the notation x=log(P), xmin=log(Pmin), x0=log(P0), ... to determine the x = LOG_LN_CONVERTER * math.log(P) xmin = LOG_LN_CONVERTER * math.log(calc_period_from_sep(M1, M2, amin)) xmax = LOG_LN_CONVERTER * math.log(calc_period_from_sep(M1, M2, amax)) # print("M1 M2 amin amax P x xmin xmax") # print(M1, M2, amin, amax, P, x, xmin, xmax) # my $x0 = 0.15; # my $x1 = 3.5; A1 = 1.0 / ( x0 ** p * (x0 - xmin) + (x1 ** p1 - x0 ** p1) / p1 + x1 ** p * (xmax - x1) ) A0 = A1 * x0 ** p A2 = A1 * x1 ** p if x < x0: res = 3.0 / 2.0 * LOG_LN_CONVERTER * A0 elif x > x1: res = 3.0 / 2.0 * LOG_LN_CONVERTER * A2 else: res = 3.0 / 2.0 * LOG_LN_CONVERTER * A1 * x ** p return res
# print(sana12(10, 2, 10, 100, 1, 1000, math.log(10), math.log(1000), 6))
[docs]def interpolate_in_mass_izzard2012( M: Union[int, float], high: Union[int, float], low: Union[int, float] ) -> Union[int, float]: """ Function to interpolate in mass TODO: fix this function. TODO: describe the args high: at M=16.3 low: at 1.15 Args: M: mass high: low: Returns: """ log_interpolation = False if log_interpolation: return (high - low) / (math.log10(16.3) - math.log10(1.15)) * ( math.log10(M) - math.log10(1.15) ) + low else: return (high - low) / (16.3 - 1.15) * (M - 1.15) + low
[docs]def Izzard2012_period_distribution( P: Union[int, float], M1: Union[int, float], log10Pmin: Union[int, float] = 1 ) -> Union[int, float]: """ period distribution which interpolates between Duquennoy and Mayor 1991 at low mass (G/K spectral type <~1.15Msun) and Sana et al 2012 at high mass (O spectral type >~16.3Msun) This gives dN/dlogP, i.e. DM/Raghavan's Gaussian in log10P at low mass and Sana's power law (as a function of logP) at high mass TODO: fix this function Args: P: period M1: Primary star mass log10Pmin: minimum period in base log10 (optional) Returns: 'probability' of interpolated distribution function at P and M1 """ # Check if there is input and force it to be at least 1 log10Pmin //= -1.0 log10Pmin = max(-1.0, log10Pmin) # save mass input and limit mass used (M1 from now on) to fitted range Mwas = M1 M1 = max(1.15, min(16.3, M1)) # print("Izzard2012 called for M={} (trunc'd to {}), P={}\n".format(Mwas, M1, P)) # Calculate the normalisations # need to normalize the distribution for this mass # (and perhaps secondary mass) prepare_dict(distribution_constants, ["Izzard2012", M1]) if not distribution_constants["Izzard2012"][M1].get(log10Pmin): distribution_constants["Izzard2012"][M1][ log10Pmin ] = 1 # To prevent this loop from going recursive N = 200.0 # Resolution for normalisation. I hope 1000 is enough dlP = (10.0 - log10Pmin) / N C = 0 # normalisation const. for lP in np.arange(log10Pmin, 10, dlP): C += dlP * Izzard2012_period_distribution(10 ** lP, M1, log10Pmin) distribution_constants["Izzard2012"][M1][log10Pmin] = 1.0 / C # print( # "Normalization constant for Izzard2012 M={} (log10Pmin={}) is\ # {}\n".format( # M1, log10Pmin, distribution_constants["Izzard2012"][M1][log10Pmin] # ) # ) lP = math.log10(P) # log period # # fits mu = interpolate_in_mass_izzard2012(M1, -17.8, 5.03) sigma = interpolate_in_mass_izzard2012(M1, 9.18, 2.28) K = interpolate_in_mass_izzard2012(M1, 6.93e-2, 0.0) nu = interpolate_in_mass_izzard2012(M1, 0.3, -1) g = 1.0 / (1.0 + 1e-30 ** (lP - nu)) lPmu = lP - mu # print( # "M={} ({}) P={} : mu={} sigma={} K={} nu={} norm=%g\n".format( # Mwas, M1, P, mu, sigma, K, nu # ) # ) # print "FUNC $distdata{Izzard2012}{$M}{$log10Pmin} * (exp(- (x-$mu)**2/(2.0*$sigma*$sigma) ) + $K/MAX(0.1,$lP)) * $g;\n"; if (lP < log10Pmin) or (lP > 10.0): return 0 else: return ( distribution_constants["Izzard2012"][M1][log10Pmin] * (math.exp(-lPmu * lPmu / (2.0 * sigma * sigma)) + K / max(0.1, lP)) * g )
######################################################################## # Mass ratio distributions ########################################################################
[docs]def flatsections(x: float, opts: dict) -> Union[float, int]: """ Function to generate flat distributions, possibly in multiple sections Args: x: mass ratio value opts: list containing the flat sections. Which are themselves dictionaries, with keys "max": upper bound, "min": lower bound and "height": value Returns: probability of that mass ratio. """ c = 0 y = 0 for opt in opts: dc = (opt["max"] - opt["min"]) * opt["height"] # print("added flatsection ({}-{})*{} = {}\n".format( # opt['max'], opt['min'], opt['height'], dc)) c += dc if opt["min"] <= x <= opt["max"]: y = opt["height"] # print("Use this\n") c = 1.0 / c y = y * c # print("flatsections gives C={}: y={}\n",c,y) return y
# print(flatsections(1, [{'min': 0, 'max': 2, 'height': 3}])) ######################################################################## # Eccentricity distributions ######################################################################## ######################################################################## # Star formation histories ########################################################################
[docs]def cosmic_SFH_madau_dickinson2014(z): """ Cosmic star formation history distribution from Madau & Dickonson 2014 (https://arxiv.org/pdf/1403.0007.pdf) Args: z: redshift Returns: Cosmic star formation rate in Solarmass year^-1 megaparsec^-3 """ CSFH = 0.015 * ((1+z)**2.7)/(1+(((1+z)/2.9)**5.6)) return CSFH
######################################################################## # Metallicity distributions ########################################################################