bjet_mcmc.blazar_utils module

file name: blazar_utils.py

Program Purpose: Implement necessary functions for mcmc (reading in configs

and data, probability functions, etc.)

Note: All file paths are relative to Bjet_MCMC

Parameters are always listed in the following order: [delta, K, n1, n2, gamma_min, gamma_max, gamma_break, B, R]

The additional parameters for EIC are bb_temp, l_nuc, tau, blob_dist, in that order. All parameters are the logarithm of the true value except for delta, n1, and n2.

Parameter

Description

Scale

delta

Doppler factor

Linear

K

Particle density [cm^-3]

Log

n1

alpha_1 (first index)

Linear

n2

alpha_2 (second index)

Linear

gamma_min

Low-energy cutoff

Log

gamma_max

High-energy cutoff

Log

gamma_break

Energy break

Log

B

Magnetic field strength [G]

Log

R

Blob radius [cm]

Log

Additional params for EIC

Parameter

Description

Scale

bb_temp

Black body temp of disk [K]

Log

l_nuc

Nucleus luminosity [ergs/s]

Log

tau

Fraction of luminosity scattered

Log

blob_dist

Distance of blob [cm]

Log
bjet_mcmc.blazar_utils.chi_squared(params, v_data, vFv_data, err_data, name_stem=None, theta=None, redshift=None, min_freq=None, max_freq=None, torus_temp=None, torus_luminosity=None, torus_frac=None, data_folder=None, executable=None, command_params_full=None, command_params_1=None, command_params_2=None, parameter_file=None, prev_files=False, use_param_file=True, verbose=False, eic=False)

Given parameters and data, create a model and then return the chi squared value for it.

Get the chi squared value for a given set of parameters and data.

The purpose of the option to supply command params is to speed up computation– the values for the first portion of the parameters and the second portion will be the same for every run in the MCMC, and it speeds up the code significantly to pass them as arguments instead of re-creating them every time. The entire list of parameters (including the ones that are changed in the MCMC) can be supplied with command_params_full. In this case, params are ignored. Alternatively, only the constant param values are provided with command_params_1 and command_params_2.

Parameters:

  • params (Any) – The parameters of the model.

  • v_data (List or ndarray) – The observed values of the dependent variable (v). Data values for v (NOT LOG)

  • vFv_data (List or ndarray) – The observed values of the dependent variable (vFv). Data values for vFv (NOT LOG)

  • err_data (List or ndarray) – The observed values of the dependent variable errors.

  • name_stem (str) – The name stem for the output and intermediate files. (optional) Default is none; will be then set to default.

  • theta (float) – The angle of the jet with respect to the line of sight in degrees. (optional) Default is none, and it will be set to the default value of 0.57.

  • redshift (float) – The redshift of the object. (optional). Default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • min_freq (float) – The minimum frequency for the model calculations. Minimum frequency for the SED model. Default is none, where it will be set to the default value of 5.0e+7 in blazar_model.process_model.(optional)

  • max_freq (float) – The maximum frequency for the model calculations. Maximum frequency for the SED model. Default is none, where it will be set to the default value of 1.0e+29 in blazar_model.process_model.(optional)

  • torus_temp (float) – The temperature of the torus. (optional) Default is none, and it will be set to the default of 2.0e+4

  • torus_luminosity (float) – The luminosity of the torus. (optional) Default is none, and it will be set to the default of 5.5e+20

  • torus_frac (float) – The fraction of the total luminosity from the torus. (optional) Value for the fraction of the torus luminosity reprocessed isotropically. Default is none, and it will be set to the default of 9.0e-5

  • data_folder (str) – Relative path to the folder where data will be saved to. (optional)

  • executable (str) – Where the bjet executable is located. (optional)

  • command_params_full (str) – Full set of parameters to pass to the bjet exec–see the README for information. Should be length 22, 23, 28, or 29.

  • command_params_1 (str) – The settings and transformation parameters: [prev files flag, data folder, model type, redshift, hubble constant, theta]

  • command_params_2 (str) – The constant and numerical parameters: [length of the emitting region, absorption by EBL, blob distance, # of spectral points, min freq, max freq, file name prefix]

  • parameter_file (str) – Relative path of the parameter file (the file where parameters are written to when modeling, will be overwritten) Default is <PARAMETER_FOLDER>/params.txt (PARAMETER_FOLDER is in home directory)

  • prev_files (bool) – Whether bjet should create _prev files after each run; default is False

  • use_param_file (bool) – Whether bjet should be called with a parameter file or with command line args; default is False

  • verbose (bool) – Whether information on the model should be displayed; default is False

  • eic (bool) – states whether the run is eic or std; default is false (std)

Returns:

The chi-squared value for the model.

Return type:

float

bjet_mcmc.blazar_utils.chi_squared_Limit_to_err(P, func_nu_data, vFv_Limit_data)

Method to include upper and lower limits in Chi2 calculation. This will return the equivalent 1sigma error to be included in a standard Chi2 formula in a way that the resulting Chi2 will be equal to -2*ln(P). P being the UL/LL likelihood (above and below)

Parameters:

  • P (float) – likelihood value associated with a data point. For ULs/LLs P is usually set either at 0.05 or 0.95 (probability above and below UL/LL)

  • func_nu_data (numpy.array) – Model flux value a nu = nu_data

  • vFv_Limit_data (numpy.array) – measure flux UL/LL.

Returns:

The error limit calculated. standard deviation to be included in a Chi2 formula

Return type:

numpy.array

bjet_mcmc.blazar_utils.chi_squared_from_model(model_results, v_data, vFv_data, err_data)

Calculates the chi-squared value based on the given model_results, v_data, vFv_data, and err_data. Consider asymmetric error bars in the dataset

Parameters:

  • model_results (list) – The results of the model calculation. Should contain logv_all and logvFv_all. Model results (logv, logvFv, v, vFv)

  • v_data (numpy array) – The observed v values. Data values (NOT LOG)

  • vFv_data (numpy array) – The observed vFv values. Data values (NOT LOG)

  • err_data (numpy array) – The error values for vFv data. Should contain two values: the lower limit and the upper limit. [err_data_down, err_data_up] Data values (NOT LOG)

Returns:

The calculated chi-squared value.

Return type:

float

bjet_mcmc.blazar_utils.e_to_v(val)
bjet_mcmc.blazar_utils.get_random_parameters(param_min_vals=None, param_max_vals=None, alpha2_limits=None, redshift=None, tau_var=None, use_variability=True, eic=False, fixed_params=None)

Get a random array of parameters in the parameter space.

Parameters:

  • param_min_vals (numpy.array or None) – minimum values for the params in the normal order; default is None, and the values are set to the defaults in blazar_utils.min_max_parameters If not provided, default_min_max is used.

  • param_max_vals (numpy.array or None) – maximum values for the params in the normal order; default is None, and the values are set to the defaults in blazar_utils.min_max_parameters. If not provided, default_min_max is used.

  • alpha2_limits (numpy.array or None) – values for alpha2 limits; default is None, and the alpha2 limits will be set to the defaults

  • redshift (float or None) – redshift value; default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • tau_var (float or None) – time in hours for tau variability; default is None, so the log_prior function will use the default, which is 24 hours

  • use_variability (bool) – whether variability should be taken into account; default is True

  • eic (bool) – states whether the run is eic or std; default is false (std)

  • fixed_params (dict or None) – Fixed parameters. If not provided, default_min_max is used.

Returns:

Randomly generated parameters.

Return type:

numpy.array

Note

random parameters within the min/max bounds. They will be valid (gamma min < gamma break < gamma max, etc.)

The function first checks if param_min_vals and param_max_vals are provided. If not, it uses default_min_max. If param_min_vals is not provided, it assigns default_min_max[0] to param_min_vals. If param_max_vals is not provided, it assigns default_min_max[1] to param_max_vals. Then it calculates the parameter_size and adds param_min_vals with parameter_size * random numbers between 0 and 1. It checks if the generated parameters are finite using log_prior function. If not, it generates new parameters until they are valid. Finally, it returns the valid parameters.

bjet_mcmc.blazar_utils.log_prior(params, param_min_vals=None, param_max_vals=None, redshift=None, tau_var=None, use_variability=True, alpha2_limits=None, eic=False, fixed_params=None)

Using a uniform prior distribution. Return whether input params are valid. list parameters with eic: [“delta”, “K”, “n_1”, “n_2”, “gamma_min”, “gamma_max”, “gamma_break”, “B”, “R”, “bb_temp”, “l_nuc”, “tau”, “blob_dist”]

Parameters:

  • params (numpy.ndarray) – The parameter values.

  • param_min_vals (numpy.ndarray) – The minimum allowed values for each parameter. If None, the function will use default values based on the specified constraints.

  • param_max_vals (numpy.ndarray) – The maximum allowed values for each parameter. If None, the function will use default values based on the specified constraints.

  • redshift (float) – redshift value; default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • tau_var (float) – time in hours for tau variability; default is None, so the log_prior function will use the default, which is 24 hours

  • use_variability (bool) – whether variability should be taken into account; default is True

  • alpha2_limits (tuple) – The limits for the alpha2 parameter. If None, the function will use default values based on the specified constraints.

  • eic (bool) – A flag indicating whether to consider EIC constraints. Default is False.

  • fixed_params (numpy.ndarray) – The fixed parameter values. If provided, the function will replace the corresponding parameters in the ‘params’ array with the fixed values.

Returns:

0 if all parameters are valid: n1 > n2, g_min < g_max, all parameters are in range and otherwise -np.inf

Return type:

float

bjet_mcmc.blazar_utils.log_prob_from_model(model_results, v_data, vFv_data, err_data)

This is the log_prob for the modeling (bigger = better fit). It returns -0.5 * the chi squared value for the v and vFv values from the given model.

Parameters:

  • model_results (Any) – The results of the model calculation. model results (logv, logvFv, v, vFv)

  • v_data (Any) – Data for frequency

  • vFv_data (Any) – Data for energy flux

  • err_data (Any) – Data for vFv error

Returns:

The log probability. -0.5 * the chi squared value

Return type:

float

bjet_mcmc.blazar_utils.log_probability(params, v_data, vFv_data, err_data, name_stem=None, param_min_vals=None, param_max_vals=None, theta=None, redshift=None, tau_var=None, use_variability=True, min_freq=None, max_freq=None, torus_temp=None, torus_luminosity=None, torus_frac=None, data_folder=None, executable=None, command_params_full=None, command_params_1=None, command_params_2=None, unique_name=False, parameter_file=None, prev_files=False, use_param_file=True, verbose=False, eic=False, fixed_params=None)

This is the log_prob for the modeling (bigger = better fit). It returns -0.5 * the chi squared value for the v and vFv values from the model with the given parameters.

The purpose of the option to supply command params is to speed up computation– the values for the first portion of the parameters and the second portion will be the same for every run in the MCMC, and it speeds up the code significantly to pass them as arguments instead of re-creating them every time. The entire list of parameters (including the ones that are changed in the MCMC) can be supplied with command_params_full. In this case, params are ignored. Alternatively, only the constant param values are provided with command_params_1 and command_params_2.

Parameters:

  • params (numpy array) – The set of parameters.

  • v_data (numpy array) – The observed frequency data. (NOT LOG)

  • vFv_data (numpy array) – The observed frequency times flux data. (NOT LOG)

  • err_data (numpy array) – The observed error data. (NOT LOG)

  • name_stem (str) – The name stem used in file naming. (optional)

  • param_min_vals (numpy array) – The minimum values for each parameter. minimum values for the params in the normal order; default is None, and the values are set to the defaults in blazar_utils.min_max_parameters

  • param_max_vals (numpy array) – The maximum values for the params in the normal order; default is None, and the values are set to the defaults in blazar_utils.min_max_parameters

  • theta (float) – The theta value. (optional). Angle from the line of sight. Default is none, and it will be set to the default value of 0.57.

  • redshift (float) – redshift value; default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • tau_var (float) – The tau_var value. (optional). Time in hours for tau variability; default is None, so the log_prior function will use the default, which is 24 hours

  • use_variability (bool) – True to use variability, False otherwise. (default is True)

  • min_freq (float) – The minimum frequency value. (optional) Default is none, where it will be set to the default value of 5.0e+7 in blazar_model.process_model.

  • max_freq (float) – The maximum frequency value. (optional) Default is none, where it will be set to the default value of 1.0e+29 in blazar_model.process_model.

  • torus_temp (float) – The torus temperature value. (optional) Default is none, and it will be set to the default of 2.0e+4

  • torus_luminosity (float) – The torus luminosity value. (optional) Default is none, and it will be set to the default of 5.5e+20

  • torus_frac (float) – The torus fraction value. (optional) Value for the fraction of the torus luminosity reprocessed isotropically. Default is none, and it will be set to the default of 9.0e-5

  • data_folder (str) – The data folder path. (optional)

  • executable (str) – The executable file path. (optional)

  • command_params_full (list) – Full set of parameters to pass to the bjet exec–see the README for information. Should be length 22, 23, 28, or 29.

  • command_params_1 (list) – The settings and transformation parameters: [prev files flag, data folder, model type, redshift, hubble constant, theta]

  • command_params_2 (list) – The constant and numerical parameters: [length of the emitting region, absorption by EBL, blob distance, # of spectral points, min freq, max freq, file name prefix]

  • unique_name (bool) – Specifies if the name stem should be created to be unique. This uses a random number, creating a very low risk of conflicts; default is

  • parameter_file (str) – Name of parameter file. This will be created from name_stem if not provided.

  • prev_files (bool) – Whether bjet should create _prev files after each run; default is False

  • use_param_file (bool) – Whether bjet should be called with a parameter file or with command line args; default is False

  • verbose (bool) – True to enable verbose mode, False otherwise. (default is False)

  • eic (bool) – True to enable EIC mode, False otherwise. (default is False)

  • fixed_params (numpy array) – The fixed parameters. (optional)

Returns:

The log-likelihood value. -0.5 * chi squared value

Return type:

float

bjet_mcmc.blazar_utils.min_max_parameters(alpha2_limits=None, eic=False, fixed_params=None)

Get the default minimum and maximum values for all the parameters.

Parameters:

  • alpha2_limits (tuple) – A tuple of two values representing the lower and upper limits for the alpha2 parameter. Defaults to (1.5, 7.5).

  • eic (bool) – A boolean value indicating whether to include extra parameters. Defaults to False.

  • fixed_params (list or None) – A list of fixed parameter values to be removed from the parameter list. Defaults to None.

Returns:

A tuple containing two arrays representing the minimum and maximum parameter values.

Return type:

tuple

Note

both returns are in the standard order:

  • [delta, K, n1, n2, gamma_min, gamma_max, gamma_break, B, R]

  • with eic: ["delta", "K", "n_1", "n_2", "gamma_min", "gamma_max", "gamma_break", "B", "R", "bb_temp", "l_nuc", "tau", "blob_dist"]

bjet_mcmc.blazar_utils.random_defaults(walkers, param_min_vals=None, param_max_vals=None, alpha2_limits=None, redshift=None, tau_var=None, use_variability=True, eic=False, fixed_params=None)

Get the values used for the initial values for the MCMC. The defaults are random values in the acceptable range that satisfy the log_prior criteria.

Parameters:

  • walkers (int) – number of walkers (specifies how many defaults to generate)

  • param_min_vals (numpy.array or None) – Minimum values for the parameters. If None, default values will be used. (in the standard order)

  • param_max_vals (numpy.array or None) – Maximum values for the parameters. If None, default values will be used. (in the standard order)

  • alpha2_limits (tuple or None) – Limits for the alpha2 parameter. If None, default values will be used.

  • redshift (float or None) – Redshift value. Default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • tau_var (float or None) – Variability value for tau. default is None, so the log_prior function will use the default, which is 24 hours

  • use_variability (bool) – Flag to indicate whether variability should be used. Default is True.

  • eic (bool) – Flag to indicate whether EIC (Extended Inverted Chirality) should be used. Default is False.

  • fixed_params (numpy.array or None) – Values for fixed parameters. If None, default values will be used.

Returns:

A numpy array with the random default parameter values for each walker. 2D np array of floats with NUM_DIM rows (the NUM_DIM parameters) and <walkers> rows. default values for all NUM_DIM parameters for each walker

Return type:

numpy.array

bjet_mcmc.blazar_utils.random_eic_from_std(std_values, walkers, param_min_vals=None, param_max_vals=None, redshift=None, tau_var=None, use_variability=True)

Given a current state of a chain for a non-EIC run (with 9 free parameters), fill in the last 4 parameters with random defaults. This is used when an EIC run that starts with non-EIC values approximated is desired.

Parameters:

  • std_values (list of lists) – 2D np array with # of walkers rows and 9 columns the current state for the std defaults

  • walkers (int) – number of walkers (specifies how many defaults to generate)

  • param_min_vals (list, optional) – 1D np array of NUM_DIM floats minimum values (in the standard order)

  • param_max_vals (list, optional) – 1D np array of NUM_DIM floats maximum values (in the standard order)

  • redshift (float, optional) – redshift value; default is None, so the log_prior function will use the default, which is 0.143 (the value for J1010)

  • tau_var (float, optional) – time in hours for tau variability; default is None, so the log_prior function will use the default, which is 24 hours.

  • use_variability (bool, optional) – Determines whether to use variability in generating random parameters. Defaults to True.

Returns:

2D np array of floats with NUM_DIM rows (the NUM_DIM parameters) and <walkers> rows default values for all NUM_DIM parameters for each walker

Return type:

numpy.ndarray

Raises:

  • ValueError – If the length of std_values is not equal to the number of walkers or if the length of std_values[0] is not 9.

  • ValueError – If the length of param_min_vals is not 13.

bjet_mcmc.blazar_utils.read_configs(config_file=None, config_string=None, verbose=False)

Given a relative path to a file, read the mcmc configs. If custom alpha 2 limits are not specified, the defaults of 1.5 and 7.5 will be used.

See README for format of the configs file.

This can also be used to parse a string of dictionary values of configs

Parameters:

  • config_file (str or None) – Absolute path to file with configurations; default is None, using the relative path to “mcmc_config.txt”

  • config_string (str or None) – This is a string of the form {'key': value, ...}. If a config_string is given, it will be parsed instead of reading from the file. This is useful when the values from a previous configuration dictionary are read from a file.

  • verbose (bool) – Controls whether values are shown; defaults to False

Returns:

The configurations as a dictionary.

Return type:

dict

bjet_mcmc.blazar_utils.read_data(data_file, cols=(0, 1, 4), use_E=True, verbose=False, instrument=False)

Read the data for the SED (energy data, vFv, flux). By default, the data file satisfies the following: - Frequency or energy data in column 0, vFv data in column 1, and vFv error data in column 4 - Space-separated (not changeable) - First row is a header (data not read from first row) - The first (#0) column has energy data in eV (which is then converted into frequency)

Parameters:

  • data_file (str) – relative path of the file w/ the data

  • cols (tuple[int]) – columns with v or E data (v_data * h), vFv_data, and err_data; default is (0, 1, 4)

  • use_E (bool) – specifies if the first (#0) column has v data or E data; default is True

  • verbose (bool) – Flag indicating whether to print additional information during execution. Default is False.

  • instrument (bool) – return the instrument use for each data point in addidion to a better data format for plotting the SED; default is False

Returns:

tuple of 3 1D np arrays of floats v_data, vFv_data, err_data

Return type:

tuple

bjet_mcmc.blazar_utils.v_to_e(val)