processing.py

import numpy as np
import pandas as pn
import random
import json
import poststat


class Process:
    fit_parameters = {'theta_0': {'log10_prior': True},
                      'E_K_iso': {'log10_prior': True},
                      'n_ref': {'log10_prior': True},
                      'theta_obs_frac': {'log10_prior': False},
                      'p': {'log10_prior': False},
                      'epsilon_B': {'log10_prior': True},
                      'epsilon_e_bar': {'log10_prior': True},
                      'a': {'log10_prior': True},
                      'l1': {'log10_prior': True},
                      'l2': {'log10_prior': True},
                      'w': {'log10_prior': True},
                      'A_V': {'log10_prior': False}}

    parameters = {'theta_0': {'log10': True},
                  'E_K_iso': {'log10': True},
                  'n_ref': {'log10': True},
                  'theta_obs_frac': {'log10': False},
                  'p': {'log10': False},
                  'epsilon_B': {'log10': True},
                  'epsilon_e_bar': {'log10': True},
                  'A_V': {'log10': False},
                  'z': {'log10': True},
                  'E_gamma_iso': {'log10': True},
                  'E_K_true': {'log10': True},
                  'epsilon_gamma': {'log10': True},
                  'M_BH': {'log10': True},
                  'epsilon_e': {'log10': True},
                  'E_gamma_true': {'log10': True},
                  'E_true': {'log10': True},
                  't_90': {'log10': True},
                  'E_K_iso_over_n_ref': {'log10': True},
                  'v_*': {'log10': False}}

    def __init__(self, ism_fits_path='ism/', wind_fits_path='wind/',
                 analysis_filename='analysis.json', credible_level='68'):
        self.ism_fits_path = ism_fits_path
        self.wind_fits_path = wind_fits_path
        self.analysis_filename = analysis_filename
        self.credible_level = credible_level

        self.prompt_features = pn.read_csv('features.txt', delimiter=',', comment='#', header=None,
                                           names=['Name', 'z', 'E_gamma_iso', 'dE_gamma_iso', 't_90', 'dt_90',
                                                  'alpha', 'dalpha', 'E_peak', 'dE_peak', 'Instrument'])

        self.lgrb = ['970508', '980703', '990510', '991208', '991216', '000301C', '000418', '000926', '010222',
                     '030329', '050820A', '050904', '060418', '080319B', '090323', '090328', '090423', '090902B',
                     '090926A', '120521C',
                     '130427A', '130702A', '130907A', '140304A', '140311A']
        self.sgrb = ['051221A', '130603B', '140903A', '200522A']
        self.allgrbs = self.lgrb + self.sgrb
        self.bad_fits = ['080319B', '090323', '140311A']
        self.wind_grb_names = []
        self.ism_grb_names = []
        self.not_sure = []

        with open(self.ism_fits_path + self.analysis_filename) as json_data_file:
            self.ism_analysis_dict = json.load(json_data_file)
        with open(self.wind_fits_path + self.analysis_filename) as json_data_file:
            self.wind_analysis_dict = json.load(json_data_file)
        self.combined_analysis_dict = self.ism_analysis_dict

        f_evidence = open('evidence_table.txt', 'w')
        for grb_name in self.allgrbs:
            ism_evidence = self.ism_analysis_dict[grb_name]['global evidence']
            ism_evidence_error = self.ism_analysis_dict[grb_name]['global evidence error']
            wind_evidence = self.wind_analysis_dict[grb_name]['global evidence']
            wind_evidence_error = self.wind_analysis_dict[grb_name]['global evidence error']
            ratio = np.exp(max(ism_evidence, wind_evidence) - min(ism_evidence, wind_evidence))
            f_evidence.write(r'{:s}, ${:.2f} \pm {:.2f}$ & ${:.2f} \pm {:.2f}$, ${:.2f}$'.format(grb_name,
                                                                                                 ism_evidence,
                                                                                                 ism_evidence_error,
                                                                                                 wind_evidence,
                                                                                                 wind_evidence_error,
                                                                                                 ratio)
                             + '\n')

            wind = True
            if ism_evidence > wind_evidence:
                wind = False
                self.ism_grb_names.append(grb_name)

            if grb_name in self.sgrb:
                wind = False
                if grb_name not in self.ism_grb_names:
                    self.ism_grb_names.append(grb_name)
            else:
                if wind:
                    self.wind_grb_names.append(grb_name)

            evidence_difference = 100 * (wind_evidence - ism_evidence) if wind else (
                    ism_evidence - wind_evidence)
            if evidence_difference < 2:
                self.not_sure.append(grb_name)
            if not wind:
                self.combined_analysis_dict[grb_name] = self.ism_analysis_dict[grb_name]
            else:
                self.combined_analysis_dict[grb_name] = self.wind_analysis_dict[grb_name]
        print('Long-GRBs ({:d}):'.format(len(self.lgrb)), self.lgrb)
        print('Short-GRBs ({:d}):'.format(len(self.sgrb)), self.sgrb)
        print('ISM type GRBs ({:d}):'.format(len(self.ism_grb_names)), self.ism_grb_names)
        print('Wind type GRBs ({:d}):'.format(len(self.wind_grb_names)), self.wind_grb_names)
        print('{:.3f}% ISM; {:.3f}% Wind'.format(
            100. * float(len(self.ism_grb_names) /
                         (len(self.ism_grb_names) + len(self.wind_grb_names))),
            100. * float(len(self.wind_grb_names) /
                         (len(self.ism_grb_names) + len(self.wind_grb_names)))))
        f_evidence.close()
        self.wlgrbs = [x for x in self.wind_grb_names if x not in (self.bad_fits + self.sgrb)]
        self.hlgrbs = [x for x in self.ism_grb_names if x not in (self.bad_fits + self.sgrb)]
        pass

    def physical_parameter_value(self, grb_name: str, parameter_name: str, posterior: bool, n_sample: int, clean=False):
        mu = None
        std_lo = None
        std_up = None
        if posterior and parameter_name in self.fit_parameters:
            n_posterior = len(self.combined_analysis_dict[grb_name]['theta_0']['posterior'])
            if n_sample > 0:
                indices = random.sample(range(n_posterior), n_sample)
            else:
                indices = range(n_posterior)
            std_lo = np.zeros(len(indices))
            std_up = np.zeros(len(indices))
            if parameter_name in self.fit_parameters:
                mu = np.array(self.combined_analysis_dict[grb_name][parameter_name]['posterior'])[indices]
            if n_sample == 1:
                mu = mu[0]
                std_up = 0
                std_lo = 0

        else:
            if parameter_name in self.fit_parameters:
                islog10 = self.combined_analysis_dict[grb_name][parameter_name]['isLog10']
                # if not getdist:
                if islog10:
                    posterior = np.log10(
                        np.array(self.combined_analysis_dict[grb_name][parameter_name]['posterior'])).reshape(
                        (-1, 1))
                else:
                    posterior = np.array(
                        self.combined_analysis_dict[grb_name][parameter_name]['posterior']).reshape(
                        (-1, 1))
                analyzer = poststat.Analyzer(posterior, parameter_names=[parameter_name])
                credible_level = float(self.credible_level) / 100.
                result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                mu_local = result['mode']
                std_lo_local = np.abs(result['uncertainty_neg'])
                std_up_local = np.abs(result['uncertainty_pos'])
                # else:
                #     mu_local = self.combined_analysis_dict[grb_name][parameter_name]['mean']
                #     std_lo_local = self.combined_analysis_dict[grb_name][parameter_name][
                #         self.credible_level + '% lower']
                #     std_up_local = self.combined_analysis_dict[grb_name][parameter_name][
                #         self.credible_level + '% upper']
                mu = mu_local if not islog10 else 10 ** mu_local
                std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
            elif parameter_name in self.prompt_features:
                mask = self.prompt_features['Name'] == grb_name
                mu = self.prompt_features[parameter_name][mask].values[0]
                if parameter_name == 'z':
                    mu = mu + 1
                    std_lo = 0
                    std_up = 0
                else:
                    std_lo = self.prompt_features['d' + parameter_name][mask].values[0]
                    std_up = self.prompt_features['d' + parameter_name][mask].values[0]
            else:
                if parameter_name == 'E_K_true':
                    if posterior:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                        post_theta_0, _, _ = self.physical_parameter_value(grb_name, 'theta_0',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                    else:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                        post_theta_0, _, _ = self.physical_parameter_value(grb_name, 'theta_0',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                    post = post_e_k_iso * (1 - np.cos(post_theta_0))
                    if posterior:
                        if hasattr(post, '__len__'):
                            return post, np.zeros(len(post)), np.zeros(len(post))
                        else:
                            return post, 0, 0
                    islog10 = self.parameters[parameter_name]['log10']
                    post = np.log10(post) if islog10 else post
                    analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                    credible_level = float(self.credible_level) / 100.
                    result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                    mu_local = result['mode']
                    std_lo_local = np.abs(result['uncertainty_neg'])
                    std_up_local = np.abs(result['uncertainty_pos'])
                    mu = mu_local if not islog10 else 10 ** mu_local
                    std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                    std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
                elif parameter_name == 'E_gamma_true':
                    post_e_gamma_iso, _, _ = self.physical_parameter_value(grb_name, 'E_gamma_iso',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                    post_theta_0, _, _ = self.physical_parameter_value(grb_name, 'theta_0',
                                                                       posterior=True, n_sample=-1,
                                                                       clean=clean)
                    post = post_e_gamma_iso * (1 - np.cos(post_theta_0))
                    islog10 = self.parameters[parameter_name]['log10']
                    post = np.log10(post) if islog10 else post
                    analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                    credible_level = float(self.credible_level) / 100.
                    result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                    mu_local = result['mode']
                    std_lo_local = np.abs(result['uncertainty_neg'])
                    std_up_local = np.abs(result['uncertainty_pos'])
                    mu = mu_local if not islog10 else 10 ** mu_local
                    std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                    std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)

                elif parameter_name == 'E_true':
                    if posterior:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                        post_e_gamma_iso, _, _ = self.physical_parameter_value(grb_name, 'E_gamma_iso',
                                                                               posterior=True, n_sample=n_sample,
                                                                               clean=clean)
                        post_theta_0, _, _ = self.physical_parameter_value(grb_name, 'theta_0',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                    else:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                        post_e_gamma_iso, _, _ = self.physical_parameter_value(grb_name, 'E_gamma_iso',
                                                                               posterior=True, n_sample=-1,
                                                                               clean=clean)
                        post_theta_0, _, _ = self.physical_parameter_value(grb_name, 'theta_0',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                    post = (post_e_k_iso + post_e_gamma_iso) * (1 - np.cos(post_theta_0))
                    if posterior:
                        if hasattr(post, '__len__'):
                            return post, np.zeros(len(post)), np.zeros(len(post))
                        else:
                            return post, 0, 0
                    islog10 = self.parameters[parameter_name]['log10']
                    post = np.log10(post) if islog10 else post
                    analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                    credible_level = float(self.credible_level) / 100.
                    result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                    mu_local = result['mode']
                    std_lo_local = np.abs(result['uncertainty_neg'])
                    std_up_local = np.abs(result['uncertainty_pos'])
                    mu = mu_local if not islog10 else 10 ** mu_local
                    # std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                    # std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
                    std_lo = std_lo_local if not islog10 else mu - 10 ** (mu_local - std_lo_local)
                    std_up = std_up_local if not islog10 else 10 ** (mu_local + std_up_local) - mu
                elif parameter_name == 'epsilon_gamma':
                    if posterior:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                        post_e_g_iso, _, _ = self.physical_parameter_value(grb_name, 'E_gamma_iso',
                                                                           posterior=True, n_sample=n_sample,
                                                                           clean=clean)
                    else:
                        post_e_k_iso, _, _ = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                        post_e_g_iso, _, _ = self.physical_parameter_value(grb_name, 'E_gamma_iso',
                                                                           posterior=True, n_sample=-1,
                                                                           clean=clean)
                    post_e_total_iso = post_e_k_iso + post_e_g_iso
                    post = post_e_g_iso / post_e_total_iso
                    if posterior:
                        if hasattr(post, '__len__'):
                            return post, np.zeros(len(post)), np.zeros(len(post))
                        else:
                            return post, 0, 0
                    islog10 = self.parameters[parameter_name]['log10']
                    post = np.log10(post) if islog10 else post
                    analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                    credible_level = float(self.credible_level) / 100.
                    result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                    mu_local = result['mode']
                    std_lo_local = np.abs(result['uncertainty_neg'])
                    std_up_local = np.abs(result['uncertainty_pos'])
                    mu = mu_local if not islog10 else 10 ** mu_local
                    std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                    std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
                elif parameter_name == 'epsilon_e':
                    if posterior:
                        post_epsilon_e_bar, _, _ = self.physical_parameter_value(grb_name, 'epsilon_e_bar',
                                                                                 posterior=True, n_sample=n_sample,
                                                                                 clean=clean)
                        post_p, _, _ = self.physical_parameter_value(grb_name, 'p',
                                                                     posterior=True, n_sample=n_sample,
                                                                     clean=clean)
                    else:
                        post_epsilon_e_bar, _, _ = self.physical_parameter_value(grb_name, 'epsilon_e_bar',
                                                                                 posterior=True, n_sample=-1,
                                                                                 clean=clean)
                        post_p, _, _ = self.physical_parameter_value(grb_name, 'p',
                                                                     posterior=True, n_sample=-1,
                                                                     clean=clean)
                    mu_p, std_lo_p, std_up_p = self.physical_parameter_value(grb_name, 'p',
                                                                             posterior=False, n_sample=-1,
                                                                             clean=clean)
                    mask = post_p > 2
                    if mu_p > 2:
                        post = (post_p[mask] - 1) / (post_p[mask] - 2) * post_epsilon_e_bar[mask]
                        if posterior:
                            if hasattr(post, '__len__'):
                                return post, np.zeros(len(post)), np.zeros(len(post))
                            else:
                                return post, 0, 0
                        islog10 = self.parameters[parameter_name]['log10']
                        post = np.log10(post) if islog10 else post
                        analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                        credible_level = float(self.credible_level) / 100.
                        result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                        mu_local = result['mode']
                        std_lo_local = np.abs(result['uncertainty_neg'])
                        std_up_local = np.abs(result['uncertainty_pos'])
                        mu = mu_local if not islog10 else 10 ** mu_local
                        std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                        std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
                    else:
                        mu = np.nan
                        std_lo = 0
                        std_up = 0
                elif parameter_name == 'M_BH':
                    post_e_true, _, _ = self.physical_parameter_value(grb_name, 'E_true',
                                                                             posterior=True, n_sample=-1,
                                                                             clean=clean)
                    a = 0.9
                    nu = 0.1
                    post_e_rot = post_e_true / nu
                    term = 2 - ((1 + (1 - a ** 2) ** 0.5) ** 2 + a ** 2) ** 0.5
                    c = 2.99792458e10
                    m_solar = 1.989e+33
                    post = 2 * post_e_rot / term / c ** 2 / m_solar
                    if posterior:
                        if hasattr(post, '__len__'):
                            return post, np.zeros(len(post)), np.zeros(len(post))
                        else:
                            return post, 0, 0
                    islog10 = self.parameters[parameter_name]['log10']
                    post = np.log10(post) if islog10 else post
                    analyzer = poststat.Analyzer(post.reshape(-1, 1), parameter_names=[parameter_name])
                    credible_level = float(self.credible_level) / 100.
                    result = analyzer.marginal_statistics(0, probability=credible_level, clean=clean)
                    mu_local = result['mode']
                    std_lo_local = np.abs(result['uncertainty_neg'])
                    std_up_local = np.abs(result['uncertainty_pos'])
                    mu = mu_local if not islog10 else 10 ** mu_local
                    std_lo = std_lo_local if not islog10 else std_lo_local * mu / np.log10(np.e)
                    std_up = std_up_local if not islog10 else std_up_local * mu / np.log10(np.e)
                elif parameter_name == 'E_K_iso_over_n_ref':
                    mu_e_k_iso, std_lo_e_k_iso, std_up_e_k_iso = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                                               posterior, n_sample,
                                                                                               clean=clean)
                    mu_n_ref, std_lo_n_ref, std_up_n_ref = self.physical_parameter_value(grb_name, 'n_ref',
                                                                                         posterior, n_sample,
                                                                                         clean=clean)
                    if grb_name in self.wind_grb_names:
                        k = 2
                    else:
                        k = 0
                    mu = (mu_e_k_iso / (mu_n_ref * 1e52)) ** (1. / (3. - k))
                    if k == 2:
                        std_lo = 1. / (mu_n_ref * 1e52) * std_lo_e_k_iso + mu_e_k_iso / (
                                mu_n_ref * 1e52) ** 2 * std_up_n_ref
                        std_up = 1. / (mu_n_ref * 1e52) * std_up_e_k_iso + mu_e_k_iso / (
                                mu_n_ref * 1e52) ** 2 * std_lo_n_ref
                    else:
                        std_lo = 1. / 3. * (mu_n_ref * 1e52) ** (-1. / 3.) * mu_e_k_iso ** (-2. / 3.) * std_lo_e_k_iso + \
                                 1. / 3. * (mu_e_k_iso / 1e52) ** (1. / 3.) * mu_n_ref ** (-4. / 3.) * std_up_n_ref
                        std_up = 1. / 3. * (mu_n_ref * 1e52) ** (-1. / 3.) * mu_e_k_iso ** (-2. / 3.) * std_up_e_k_iso + \
                                 1. / 3. * (mu_e_k_iso / 1e52) ** (1. / 3.) * mu_n_ref ** (-4. / 3.) * std_lo_n_ref
                elif parameter_name == 'v_*':
                    mu_e_k_iso, std_lo_e_k_iso, std_up_e_k_iso = self.physical_parameter_value(grb_name, 'E_K_iso',
                                                                                               posterior, n_sample,
                                                                                               clean=clean)
                    mu_e_k_iso_52 = mu_e_k_iso / 1e52
                    std_lo_e_k_iso_52 = std_lo_e_k_iso / 1e52
                    std_up_e_k_iso_52 = std_up_e_k_iso / 1e52
                    mu_n_ref, std_lo_n_ref, std_up_n_ref = self.physical_parameter_value(grb_name, 'n_ref',
                                                                                         posterior, n_sample,
                                                                                         clean=clean)

                    mu_t_90, std_lo_t_90, std_up_t_90 = self.physical_parameter_value(grb_name, 't_90',
                                                                                      posterior, n_sample,
                                                                                      clean=clean)
                    mu_z, std_lo_z, std_up_z = self.physical_parameter_value(grb_name, 'z',
                                                                             posterior, n_sample,
                                                                             clean=clean)
                    m_dot = 1e-5
                    v_wind = 1000
                    gamma_0 = 300

                    m_dot_cgs = m_dot * 1.989e33 / 3.154e7
                    v_wind_cgs = v_wind * 1e5
                    a_0 = m_dot_cgs / (4 * np.pi * v_wind_cgs)
                    c = 2.99792458e10
                    m_p = 1.6726219e-24
                    a_star = 0.01
                    a = a_star * a_0
                    mu_rho_ref = mu_n_ref * m_p

                    t_dec = mu_e_k_iso / (4 * np.pi * a * gamma_0 ** 4 * c ** 3)
                    r_dec = 2 * gamma_0 ** 2 * c * t_dec
                    # r_other = mu_t_90 * c
                    ratio = 10  # wind / star
                    r_standoff = np.sqrt(a / mu_rho_ref) * ratio
                    distance = r_standoff - r_dec
                    term = 1.1e17 * ((1 + mu_z) / 2) ** -0.5 * (mu_e_k_iso_52 / a_star) ** 0.5
                    t_cross = (distance / term) ** 2 + t_dec / 86400
                    # mu = 4 * np.pi * a ** (3. / 2.) * gamma_0 ** 4 * c ** 2 / (mu_e_k_iso * np.sqrt(mu_rho_ref))
                    # mu = np.sqrt(a / mu_rho_ref) / max(r_dec, r_other)
                    mu = t_cross
                    std_lo = 0
                    std_up = 0

        return mu, std_lo, std_up

    def get_dictionary(self, inc_grbs=[], exc_grbs=[], parameter_names=[], posterior=False, n_sample=-1,
                       return_names=False, clean=False):
        grb_names = [x for x in inc_grbs if x not in exc_grbs]
        grb_dict = {}
        for parameter_name in parameter_names:
            fancy_param_name = self.fancy_parameter_names(parameter_name, self.parameters[parameter_name]['log10'])
            grb_dict[fancy_param_name] = {}
            grb_dict[fancy_param_name]['mean'] = []
            grb_dict[fancy_param_name]['std_up'] = []
            grb_dict[fancy_param_name]['std_lo'] = []
        for grb_name in grb_names:
            for parameter_name in parameter_names:
                islog10 = self.parameters[parameter_name]['log10']
                fancy_param_name = self.fancy_parameter_names(parameter_name, islog10)

                mu_local, std_lo_local, std_up_local = self.physical_parameter_value(grb_name, parameter_name,
                                                                                     posterior, n_sample, clean=clean)
                if islog10:
                    if not hasattr(mu_local, '__len__'):
                        mu = np.log10(mu_local) if mu_local != 0 else np.nan
                        if mu_local != 0:
                            std_lo = np.abs(np.log10(np.e) * std_lo_local / mu_local) if std_lo_local != 0 else 0
                            std_up = np.abs(np.log10(np.e) * std_up_local / mu_local) if std_up_local != 0 else 0
                        else:
                            std_lo = 0
                            std_up = 0
                    else:
                        mu = np.log10(mu_local)
                        std_lo = np.zeros(len(mu))  # std_lo_local
                        std_up = np.zeros(len(mu))  # std_up_local
                else:
                    mu = mu_local
                    std_lo = np.abs(std_lo_local)
                    std_up = np.abs(std_up_local)
                grb_dict[fancy_param_name]['mean'].append(mu)
                grb_dict[fancy_param_name]['std_lo'].append(std_lo)
                grb_dict[fancy_param_name]['std_up'].append(std_up)
        if not return_names:
            return grb_dict
        else:
            return grb_dict, grb_names

    def get_mean_array(self, inc_grbs=[], exc_grbs=[], parameter_names=[], posterior=False, n_sample=100,
                       clean=False, return_names=False):
        temp = self.get_dictionary(inc_grbs=inc_grbs, exc_grbs=exc_grbs, parameter_names=parameter_names,
                                   posterior=posterior, n_sample=n_sample, return_names=return_names,
                                   clean=clean)
        if return_names:
            grb_dict, names = temp
        else:
            grb_dict = temp
        result = []
        for parameter_name in parameter_names:
            islog10 = self.parameters[parameter_name]['log10']
            fancy_param_name = self.fancy_parameter_names(parameter_name, islog10)
            result.append(np.array(grb_dict[fancy_param_name]['mean']))
        if return_names:
            return np.array(result), names
        else:
            return np.array(result)

    def get_std_lo_array(self, inc_grbs=[], exc_grbs=[], parameter_names=[], posterior=False, n_sample=100,
                         clean=False, return_names=False):
        temp = self.get_dictionary(inc_grbs=inc_grbs, exc_grbs=exc_grbs, parameter_names=parameter_names,
                                   posterior=posterior, n_sample=n_sample, return_names=return_names,
                                   clean=clean)
        if return_names:
            grb_dict, names = temp
        else:
            grb_dict = temp
        result = []
        for parameter_name in parameter_names:
            islog10 = self.parameters[parameter_name]['log10']
            fancy_param_name = self.fancy_parameter_names(parameter_name, islog10)
            result.append(np.array(grb_dict[fancy_param_name]['std_lo']))
        if return_names:
            return np.array(result), names
        else:
            return np.array(result)

    def get_std_up_array(self, inc_grbs=[], exc_grbs=[], parameter_names=[], posterior=False, n_sample=100,
                         clean=False, return_names=False):
        temp = self.get_dictionary(inc_grbs=inc_grbs, exc_grbs=exc_grbs, parameter_names=parameter_names,
                                   posterior=posterior, n_sample=n_sample, return_names=return_names,
                                   clean=clean)
        if return_names:
            grb_dict, names = temp
        else:
            grb_dict = temp
        result = []
        for parameter_name in parameter_names:
            islog10 = self.parameters[parameter_name]['log10']
            fancy_param_name = self.fancy_parameter_names(parameter_name, islog10)
            result.append(np.array(grb_dict[fancy_param_name]['std_up']))
        if return_names:
            return np.array(result), names
        else:
            return np.array(result)

    @staticmethod
    def make_table(grb_names, dictionary, filename='table.txt'):
        file = open(filename, 'w')
        for i, grb_name in enumerate(grb_names):
            file.write('{:s}'.format(grb_name))
            for param_name in dictionary:
                mu = dictionary[param_name]['mean'][i]
                std_lo = dictionary[param_name]['std_lo'][i]
                std_up = dictionary[param_name]['std_up'][i]
                file.write(r' & ${:.2f}'.format(mu) + '^{+' + '{:.2f}'.format(std_up) + '}_{-' + '{:.2f}'.format(
                    std_lo) + '}$')
            file.write('\n')
        file.close()

    @staticmethod
    def fancy_parameter_names(parameter_name: str, islog10: bool):
        if parameter_name == 'theta_0':
            fancy_name = r'$\theta_0$ (rad)'
            log10_fancy_name = r'$\log_{10}\theta_0$ [rad]'
        elif parameter_name == 'E_K_iso':
            fancy_name = r'$E_{K, \mathrm{iso}}$ (erg)'
            log10_fancy_name = r'$\log_{10}E_{K, \mathrm{iso}}$ [erg]'
        elif parameter_name == 'n_ref':
            fancy_name = r'$n_{\mathrm{ref}} \ (\mathrm{cm^{-3}})$'
            log10_fancy_name = r'$\log_{10}n_{\mathrm{ref}} \ [\mathrm{cm^{-3}}]$'
        elif parameter_name == 'theta_obs_frac':
            fancy_name = r'$\theta_{\mathrm{obs}} / \theta_0$'
            log10_fancy_name = r'$\log_{10}\theta_{\mathrm{obs}} / \theta_0$'
        elif parameter_name == 'p':
            fancy_name = r'$p$'
            log10_fancy_name = r'$\log_{10}p$'
        elif parameter_name == 'epsilon_B':
            fancy_name = r'$\epsilon_B$'
            log10_fancy_name = r'$\log_{10}\epsilon_B$'
        elif parameter_name == 'epsilon_e_bar':
            fancy_name = r'$\bar{\epsilon}_e$'
            log10_fancy_name = r'$\log_{10}\bar{\epsilon}_e$'
        elif parameter_name == 'z':
            fancy_name = r'$z + 1$'
            log10_fancy_name = r'$\log_{10}(z + 1)$'
        elif parameter_name == 'A_V':
            fancy_name = r'$A_V$'
            log10_fancy_name = r'$\log_{10}A_V$'
        elif parameter_name == 'E_K_true':
            fancy_name = r'$E_{K, \mathrm{true}}$ (erg)'
            log10_fancy_name = r'$\log_{10}E_{K, \mathrm{true}}$ [erg]'
        elif parameter_name == 'E_true':
            fancy_name = r'$E_{\gamma, \mathrm{true}} + E_{K, \mathrm{true}}$ (erg)'
            log10_fancy_name = r'$\log_{10}(E_{\gamma, \mathrm{true}} + E_{K, \mathrm{true}})$ [erg]'
        elif parameter_name == 'E_gamma_iso':
            fancy_name = r'$E_{\gamma, \mathrm{iso}}$ (erg)'
            log10_fancy_name = r'$\log_{10}E_{\gamma, \mathrm{iso}}$ [erg]'
        elif parameter_name == 'E_gamma_true':
            fancy_name = r'$E_{\gamma, \mathrm{true}}$ (erg)'
            log10_fancy_name = r'$\log_{10}E_{\gamma, \mathrm{true}}$ [erg]'
        elif parameter_name == 'epsilon_gamma':
            fancy_name = r'$\epsilon_\gamma$'
            log10_fancy_name = r'$\log_{10}\epsilon_\gamma$'
        elif parameter_name == 'M_BH':
            fancy_name = r'$M_{\mathrm{BH}}\ (M_\odot)$'
            log10_fancy_name = r'$\log_{10}M_{\mathrm{BH}}\ [M_\odot]$'
        elif parameter_name == 'epsilon_e':
            fancy_name = r'$\epsilon_e$'
            log10_fancy_name = r'$\log_{10}\epsilon_e$'
        elif parameter_name == 't_90':
            fancy_name = r'$t_{90}\ \mathrm{s}$'
            log10_fancy_name = r'$\log_{10}t_{90}\ [\mathrm{s}]$'
        elif parameter_name == 'E_K_iso_over_n_ref':
            fancy_name = r'$E_{K, \mathrm{iso}} / n_\mathrm{ref}$'
            log10_fancy_name = r'$\log_{10}(E_{K, \mathrm{iso}} / n_\mathrm{ref})$'
        elif parameter_name == 'v_*':
            fancy_name = r'$t_\times \ (\mathrm{days})$'
            log10_fancy_name = r'$\log_{10}t_\times \ (\mathrm{days})$'
            # fancy_name = r'$v_* \ (\mathrm{km / s})$'
            # log10_fancy_name = r'$\log_{10}v_* \ (\mathrm{km / s})$'
        else:
            fancy_name = parameter_name
            log10_fancy_name = 'log10_' + parameter_name
        return fancy_name if not islog10 else log10_fancy_name