<commit_before># -*- coding: utf-8 -*-
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import multiprocessing
import warnings
from collections.abc import Iterable

import astropy.units as u
import emcee
import numpy as np
from astropy import log

from .utils import sed_conversion, validate_data_table

__all__ = [
    "normal_prior",
    "uniform_prior",
    "log_uniform_prior",
    "get_sampler",
    "run_sampler",
]

# Define phsyical types used in plot and utils.validate_data_table
u.def_physical_type(u.erg / u.cm ** 2 / u.s, "flux")
u.def_physical_type(u.Unit("1/(s cm2 erg)"), "differential flux")
u.def_physical_type(u.Unit("1/(s erg)"), "differential power")
u.def_physical_type(u.Unit("1/TeV"), "differential energy")
u.def_physical_type(u.Unit("1/cm3"), "number density")
u.def_physical_type(u.Unit("1/(eV cm3)"), "differential number density")

# Prior functions


def uniform_prior(value, umin, umax):
    """Uniform prior distribution."""
    if umin <= value <= umax:
        return 0.0
    else:
        return -np.inf


def normal_prior(value, mean, sigma):
    """Normal prior distribution."""
    return -0.5 * (2 * np.pi * sigma) - (value - mean) ** 2 / (2.0 * sigma)


def log_uniform_prior(value, umin=0, umax=None):
    """Log-uniform prior distribution."""
    if value > 0 and value >= umin:
        if umax is not None:
            if value <= umax:
                return 1 / value
            else:
                return -np.inf
        else:
            return 1 / value
    else:
        return -np.inf


# Probability function


def lnprobmodel(model, data):

    # Check if conversion is required
    model_is_sed = model.unit.physical_type in ["power", "flux"]
    data_is_sed = data["flux"].unit.physical_type in ["power", "flux"]

    if model_is_sed != data_is_sed:
        unit, sed_factor = sed_conversion(
            data["energy"], model.unit, data_is_sed
        )
        model = (model * sed_factor).to(data["flux"].unit)

    ul = data["ul"]
    notul = ~ul

    difference = model[notul] - data["flux"][notul]

    # use different errors for model above or below data
    sign = difference > 0
    loerr, hierr = 1 * ~sign, 1 * sign
    logprob = -(difference ** 2) / (
        2.0
        * (
            loerr * data["flux_error_lo"][notul]
            + hierr * data["flux_error_hi"][notul]
        )
        ** 2
    )

    totallogprob = np.sum(logprob)

    if np.sum(ul) > 0:
        # deal with upper limits at CL set by data['cl']
        violated_uls = np.sum(model[ul] > data["flux"][ul])
        totallogprob += violated_uls * np.log(1.0 - data["cl"][violated_uls])

    return totallogprob


def lnprob(pars, data, modelfunc, priorfunc):

    if priorfunc is None:
        lnprob_priors = 0.0
    else:
        lnprob_priors = priorfunc(pars)

    modelout = modelfunc(pars, data)

    # Save blobs or save model if no blobs given
    if isinstance(modelout, Iterable) and not isinstance(modelout, np.ndarray):
        model = modelout[0]
        blob = modelout
    else:
        model = modelout
        # Avoid squeezing of the model output in order to retain units by
        # adding a second dummy object blob
        blob = (modelout, np.nan)

    if not np.isinf(lnprob_priors):
        lnprob_model = lnprobmodel(model, data)
        total_lnprob = lnprob_model + lnprob_priors
    else:
        total_lnprob = lnprob_priors

    return total_lnprob, *blob


# Sampler funcs


def _run_mcmc(sampler, pos, nrun):
    for i, state in enumerate(
        sampler.sample(pos, iterations=nrun, store=True)
    ):
        progress = 100.0 * float(i) / float(nrun)
        if progress % 5 < (5.0 / float(nrun)):
            print(
                "\nProgress of the run: {0:.0f} percent"
                " ({1} of {2} steps)".format(int(progress), i, nrun)
            )
            paravg, parstd = [], []
            npars = sampler.get_chain().shape[-1]
            for npar in range(npars):
                paravg.append(np.median(state.coords[:, npar]))
                parstd.append(np.std(state.coords[:, npar]))
            print(
                "                           "
                + (" ".join(["{%i:-^15}" % i for i in range(npars)])).format(
                    *sampler.labels
                )
            )
            print(
                "  Last ensemble median : "
                + (" ".join(["{%i:^15.3g}" % i for i in range(npars)])).format(<commit_msg>try to fix for Python 3.7<commit_after># -*- coding: utf-8 -*-
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import multiprocessing
import warnings
from collections.abc import Iterable

import astropy.units as u
import emcee
import numpy as np
from astropy import log

from .utils import sed_conversion, validate_data_table

__all__ = [
    "normal_prior",
    "uniform_prior",
    "log_uniform_prior",
    "get_sampler",
    "run_sampler",
]

# Define phsyical types used in plot and utils.validate_data_table
u.def_physical_type(u.erg / u.cm ** 2 / u.s, "flux")
u.def_physical_type(u.Unit("1/(s cm2 erg)"), "differential flux")
u.def_physical_type(u.Unit("1/(s erg)"), "differential power")
u.def_physical_type(u.Unit("1/TeV"), "differential energy")
u.def_physical_type(u.Unit("1/cm3"), "number density")
u.def_physical_type(u.Unit("1/(eV cm3)"), "differential number density")

# Prior functions


def uniform_prior(value, umin, umax):
    """Uniform prior distribution."""
    if umin <= value <= umax:
        return 0.0
    else:
        return -np.inf


def normal_prior(value, mean, sigma):
    """Normal prior distribution."""
    return -0.5 * (2 * np.pi * sigma) - (value - mean) ** 2 / (2.0 * sigma)


def log_uniform_prior(value, umin=0, umax=None):
    """Log-uniform prior distribution."""
    if value > 0 and value >= umin:
        if umax is not None:
            if value <= umax:
                return 1 / value
            else:
                return -np.inf
        else:
            return 1 / value
    else:
        return -np.inf


# Probability function


def lnprobmodel(model, data):

    # Check if conversion is required
    model_is_sed = model.unit.physical_type in ["power", "flux"]
    data_is_sed = data["flux"].unit.physical_type in ["power", "flux"]

    if model_is_sed != data_is_sed:
        unit, sed_factor = sed_conversion(
            data["energy"], model.unit, data_is_sed
        )
        model = (model * sed_factor).to(data["flux"].unit)

    ul = data["ul"]
    notul = ~ul

    difference = model[notul] - data["flux"][notul]

    # use different errors for model above or below data
    sign = difference > 0
    loerr, hierr = 1 * ~sign, 1 * sign
    logprob = -(difference ** 2) / (
        2.0
        * (
            loerr * data["flux_error_lo"][notul]
            + hierr * data["flux_error_hi"][notul]
        )
        ** 2
    )

    totallogprob = np.sum(logprob)

    if np.sum(ul) > 0:
        # deal with upper limits at CL set by data['cl']
        violated_uls = np.sum(model[ul] > data["flux"][ul])
        totallogprob += violated_uls * np.log(1.0 - data["cl"][violated_uls])

    return totallogprob


def lnprob(pars, data, modelfunc, priorfunc):

    if priorfunc is None:
        lnprob_priors = 0.0
    else:
        lnprob_priors = priorfunc(pars)

    modelout = modelfunc(pars, data)

    # Save blobs or save model if no blobs given
    if isinstance(modelout, Iterable) and not isinstance(modelout, np.ndarray):
        model = modelout[0]
        blob = modelout
    else:
        model = modelout
        # Avoid squeezing of the model output in order to retain units by
        # adding a second dummy object blob
        blob = (modelout, np.nan)

    if not np.isinf(lnprob_priors):
        lnprob_model = lnprobmodel(model, data)
        total_lnprob = lnprob_model + lnprob_priors
    else:
        total_lnprob = lnprob_priors

    return (total_lnprob, *blob)


# Sampler funcs


def _run_mcmc(sampler, pos, nrun):
    for i, state in enumerate(
        sampler.sample(pos, iterations=nrun, store=True)
    ):
        progress = 100.0 * float(i) / float(nrun)
        if progress % 5 < (5.0 / float(nrun)):
            print(
                "\nProgress of the run: {0:.0f} percent"
                " ({1} of {2} steps)".format(int(progress), i, nrun)
            )
            paravg, parstd = [], []
            npars = sampler.get_chain().shape[-1]
            for npar in range(npars):
                paravg.append(np.median(state.coords[:, npar]))
                parstd.append(np.std(state.coords[:, npar]))
            print(
                "                           "
                + (" ".join(["{%i:-^15}" % i for i in range(npars)])).format(
                    *sampler.labels
                )
            )
            print(
                "  Last ensemble median : "
                + (" ".join(["{%i:^15.3g}" % i for i in range(npars)])).format(
{'fraction_non_alphanumeric': 0.12690071556350627, 'fraction_numerical': 0.012075134168157423, 'mean_word_length': 2.198069360028602, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1953825', 'n_tokens_mistral': 3077, 'n_tokens_neox': 2860, 'n_words': 760}