6543527074

Committed 17 Oct 2023 06:45AM UTC coverage: 39.83% (+4.2%) from 35.621%

Build # 6543527074

Build Type

push

github

Committed by

web-flow

Commit Message

Add CI Tests (#185)

Closes: #69
Closes: #66
Closes: #47
Closes: #90

Run Details

69 of 69 new or added lines in 7 files covered. (100.0%)

797 of 2001 relevant lines covered (39.83%)

0.4 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

15.06

/hepi/plot/plot.py

from collections.abc import Iterable

import matplotlib as mpl
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
import pyslha
from matplotlib import colors
from matplotlib.image import NonUniformImage
from matplotlib.ticker import NullFormatter
from scipy import integrate as auc
from scipy.interpolate import BSpline, make_interp_spline
from smpl import plot as splot

from ..input import Input, get_output_dir, replace_macros
from ..particles import get_name


def title(
    i: Input,
    axe=None,
    scenario=None,
    diff_L_R=None,
    extra="",
    cms_energy=True,
    pdf_info=True,
    id=False,
    **kwargs
):
    """Sets the title on axis `axe`."""
    if axe is None:
        axe = plt.gca()
    axe.set_title(
        "$pp\\to"
        + get_name(i.particle1)
        + get_name(i.particle2)
        + "$"
        + (" at $\\sqrt{S} = " + str(i.energy / 1000) + "$ TeV" if cms_energy else "")
        + ((" for " + i.slha.split(".")[0]) if scenario is None else (" " +scenario))
        + (" with " + i.pdf_nlo if pdf_info else "")
        + " "
        + extra
        + ((" [" + i.id + "]") if id else "")
    )


def energy_plot(
    dict_list,
    y,
    yscale=1.0,
    xaxis="E [GeV]",
    yaxis="$\\sigma$ [pb]",
    label=None,
    **kwargs
):
    """Plot energy on the x-axis."""
    plot(
        dict_list,
        "energy",
        y,
        label=label,
        xaxis=xaxis,
        yaxis=yaxis,
        logy=True,
        yscale=yscale,
        **kwargs
    )


def combined_mass_plot(dict_list, y, part, label=None, **kwargs):
    if label == "":
        kwargs["label"] = ""
    mass_plot(
        dict_list,
        y + "_NOERR",
        part,
        next_color=False,
        fmt=" ",
        fill=False,
        interpolate_label=label,
        **kwargs
    )
    if "label" in kwargs:
        del kwargs["label"]
    mass_plot(
        dict_list,
        y + "_COMBINED",
        part,
        next_color=True,
        fit_fmt=" ",
        fill=True,
        label="",
        interpolate_lower_uncertainty=False,
        **kwargs
    )


def combined_plot(dict_list, x, y, label=None, **kwargs):
    if label == "":
        kwargs["label"] = ""
    plot(
        dict_list,
        x,
        y + "_NOERR",
        next_color=False,
        fmt=" ",
        fill=False,
        interpolate_label=label,
        **kwargs
    )
    if "label" in kwargs:
        del kwargs["label"]
    plot(
        dict_list,
        x,
        y + "_COMBINED",
        next_color=True,
        fit_fmt=" ",
        fill=True,
        label="",
        interpolate_lower_uncertainty=False,
        **kwargs
    )


def mass_plot(
    dict_list,
    y,
    part,
    logy=True,
    yaxis="$\\sigma$ [pb]",
    yscale=1.0,
    label=None,
    xaxis=None,
    **kwargs
):
    if not "mass_" + str(part) in dict_list:
        dict_list["mass_" + str(part)] = get_mass(dict_list, abs(part))
    if xaxis is None:
        xaxis = "$m_{" + get_name(part) + "}$ [GeV]"
    plot(
        dict_list,
        "mass_" + str(part),
        y,
        label=label,
        yaxis=yaxis,
        xaxis=xaxis,
        logy=logy,
        yscale=yscale,
        **kwargs
    )


def mass_vplot(
    dict_list,
    y,
    part,
    logy=True,
    yaxis="$\\sigma$ [pb]",
    yscale=1.0,
    label=None,
    mask=None,
    **kwargs
):
    vplot(
        get_mass(dict_list, part)[mask],
        y[mask],
        label=label,
        xaxis="$m_{" + get_name(part) + "}$ [GeV]",
        yaxis=yaxis,
        logy=logy,
        yscale=yscale,
        mask=mask,
        **kwargs
    )


def get_mass(l: dict, iid: int):
    """
    Get the mass of particle with id `iid` out of the list in the "slha" element in the dict.

    Returns
        :obj:`list` of float : masses of particles in each element of the dict list.

    """
    ret = []
    for s in l["slha"]:
        d = pyslha.read(get_output_dir() + s)
        ret.append(d.blocks["MASS"][abs(iid)])

    return np.array(ret)


def plot(
    dict_list,
    x,
    y,
    label=None,
    xaxis="M [GeV]",
    yaxis="$\\sigma$ [pb]",
    ratio=False,
    K=False,
    K_plus_1=False,
    logy=True,
    yscale=1.0,
    mask=None,
    **kwargs
) -> None:
    """
    Creates a plot based on the entries `x`and `y` in `dict_list`.

    Examples

    .. plot::
        :include-source:

        >>> import urllib.request
        >>> import hepi
        >>> dl = hepi.load(urllib.request.urlopen(
        ... "https://raw.githubusercontent.com/fuenfundachtzig/xsec/master/json/pp13_hino_NLO%2BNLL.json"
        ... ))
        >>> hepi.plot(dl,"N1","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{\\\\chi}_1^0}$ [GeV]")
    """
    if isinstance(y, Iterable) and not isinstance(y, str):
        for yi in y:
            plot(
                dict_list,
                x,
                yi,
                label,
                xaxis,
                yaxis,
                ratio,
                K,
                K_plus_1,
                logy,
                yscale,
                mask,
                **kwargs
            )
        return
    # TODO use kwargs
    if label is None:
        label = y.replace("_PLUS_", "+").replace("_OVER_", "/")
    if label == "":
        label = None

    vx = dict_list[x].to_numpy()[mask]
    vy = dict_list[y].to_numpy()[mask]

    if K:
        yaxis = "$K$"
        yscale = 1.0
        vy = vy / splot.unv(dict_list["LO"].to_numpy()[mask])
        if K_plus_1:
            vy = vy + vy / vy
    if ratio:
        yaxis = "Ratio"
        yscale = 1.0
        vy = vy / splot.unv(vy)

    vplot(vx, vy, label, xaxis, yaxis, logy, yscale, mask=mask, **kwargs)


def index_open(var, idx):
    if len(idx) == 1:
        return var[idx]
    return index_open(var[idx[0]], idx[1:])


def slha_data(li, index_list):
    vx = []
    for l in li:
        b = pyslha.read(get_output_dir() + l.slha)
        vx.append(index_open(b.blocks, index_list))
    return np.array(vx)


def slha_plot(li, x, y, **kwargs):
    vx = slha_data(li, x)
    vy = slha_data(li, y)

    vplot(np.array(vx), np.array(vy), **kwargs)


def vplot(
    x,
    y,
    label=None,
    xaxis="E [GeV]",
    yaxis="$\\sigma$ [pb]",
    logy=True,
    yscale=1.0,
    interpolate=True,
    plot_data=True,
    data_color=None,
    mask=-1,
    fill=False,
    data_fmt=".",
    fmt="-",
    print_area=False,
    sort=True,
    **kwargs
):
    """
    Creates a plot based on the values in `x`and `y`.

    """
    color = data_color
    if label is None:
        # label = "??"
        pass
    if mask is None:
        x = x[0]
        y = y[0]
    if sort:
        permute = x.argsort(kind="stable")
        vx = x[permute]
        vy = y[permute]
    else:
        vx = x
        vy = y

    if data_color is None and "next_color" not in kwargs:
        if "axes" in kwargs and kwargs["axes"] is not None:
            (bl,) = kwargs["axes"].plot([], [])
        else:
            (bl,) = plt.gca().plot([], [])
        color = bl.get_color()
    iii = splot.data(
        vx,
        vy * yscale,
        label=label,
        xaxis=xaxis,
        yaxis=yaxis,
        logy=logy,
        data_color=color,
        fit_color=color,
        also_data=plot_data,
        interpolate=interpolate,
        fmt=data_fmt,
        sigmas=0 if not fill else 1,
        **kwargs
    )
    if iii is not None:
        # ii = iii[0]
        ix = iii[1]
        iy = iii[2]
    if print_area:
        print(
            "computed AUC using scipy.integrate.simpson: {}".format(auc.simpson(iy, ix))
        )
    if (np.any(np.less(vy, 0)) or (interpolate and np.any(np.less(iy, 0)))) and logy:
        splot.data(
            vx,
            -vy * yscale,
            label="-" + label,
            xaxis=xaxis,
            yaxis=yaxis,
            logy=logy,
            data_color=color,
            fit_color=color,
            also_data=plot_data,
            interpolate=interpolate,
            interpolate_fmt=(0, (3, 1, 3, 1, 1, 1)),
            fmt=data_fmt,
            sigmas=0 if not fill else 1,
            **kwargs
        )


def mass_mapplot(
    dict_list,
    part1,
    part2,
    z,
    logz=True,
    zaxis="$\\sigma$ [pb]",
    zscale=1.0,
    label=None,
):
    mapplot(
        dict_list,
        "mass_" + str(part1),
        "mass_" + str(part2),
        z,
        xaxis="$M_{" + get_name(part1) + "}$ [GeV]",
        yaxis="$M_{" + get_name(part2) + "}$ [GeV]",
        zaxis=zaxis,
        logz=logz,
        zscale=zscale,
    )


def mapplot(dict_list, x, y, z, xaxis=None, yaxis=None, zaxis=None, **kwargs):
    """

    Examples

    .. plot::
        :include-source:

        >>> import urllib.request
        >>> import hepi

        >>> dl = hepi.load(urllib.request.urlopen(
        ... "https://raw.githubusercontent.com/APN-Pucky/xsec/master/json/pp13_SGmodel_GGxsec_NLO%2BNLL.json"
        ... ),dimensions=2)
        >>> hepi.mapplot(dl,"gl","sq","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{g}}$ [GeV]",yaxis="$m_{\\\\tilde{q}}$ [GeV]" , zaxis="$\\\\sigma_{\\\\mathrm{NLO+NLL}}$ [pb]")
    """
    if xaxis is None:
        xaxis = x
    if yaxis is None:
        yaxis = y
    if zaxis is None:
        zaxis = replace_macros(z)
    vx = dict_list[x]
    vy = dict_list[y]
    vz = dict_list[z]
    splot.plot2d(
        vx, vy, vz, style="pcolormesh", xaxis=xaxis, yaxis=yaxis, zaxis=zaxis, **kwargs
    )


map_vplot = lambda *a, **da: splot.plot2d(*a, style="pcolormesh", **da)
scatter_vplot = lambda *a, **da: splot.plot2d(*a, style="scatter", **da)


def scatterplot(dict_list, x, y, z, xaxis=None, yaxis=None, zaxis=None, **kwargs):
    """
    Scatter map 2d.
    Central color is the central value, while the inner and outer ring are lower and upper bounds of the uncertainty interval.

    Examples

    .. plot::
        :include-source:

        >>> import urllib.request
        >>> import hepi
        >>> dl = hepi.load(urllib.request.urlopen(
        ... "https://raw.githubusercontent.com/APN-Pucky/xsec/master/json/pp13_hinosplit_N2N1_NLO%2BNLL.json"
        ... ),dimensions=2)
        >>> hepi.scatterplot(dl,"N1","N2","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{\\\\chi}_1^0}$ [GeV]",yaxis="$m_{\\\\tilde{\\\\chi}_2^0}$ [GeV]" , zaxis="$\\\\sigma_{\\\\mathrm{NLO+NLL}}$ [pb]")

    """
    if xaxis is None:
        xaxis = x
    if yaxis is None:
        yaxis = y
    if zaxis is None:
        zaxis = replace_macros(z)
    vx = dict_list[x]
    vy = dict_list[y]
    vz = dict_list[z]
    splot.plot2d(
        vx, vy, vz, style="scatter", xaxis=xaxis, yaxis=yaxis, zaxis=zaxis, **kwargs
    )


fig = None
axs = None
lines = []
labels = []


def err_plt(axes, x, y, label=None, error=False):
    v = label
    ind = np.argsort(splot.unv(x), kind="stable")
    if error:
        l, _, _ = axes.errorbar(
            x.to_numpy()[ind], splot.unv(y)[ind], yerr=splot.usd(y), capsize=5, label=v
        )
        return l
    else:
        l = axes.plot(x.to_numpy()[ind], splot.unv(y)[ind], label=v)
        return l[0]


def scale_plot(
    dict_list,
    vl,
    seven_point_band=False,
    cont=False,
    error=True,
    li=None,
    plehn_color=False,
    yscale=1.0,
    unit="pb",
    yaxis=None,
    **kwargs
):
    """Creates a scale variance plot with 5 panels (xline)."""
    global fig, axs, lines, labels
    cycle_safe = mpl.rcParams["axes.prop_cycle"]
    if plehn_color:
        mpl.rcParams["axes.prop_cycle"] = mpl.cycler(color=["b", "r", "g"])
    if not cont:
        fig, axs = plt.subplots(1, 5, figsize=(12, 3), sharey=True)
        # Remove horizontal space between axes
        fig.subplots_adjust(wspace=0)
        if li is not None:
            title(axs[2], li[0], **kwargs)

    mr = dict_list["mu_r"]
    mf = dict_list["mu_f"]

    if not cont:
        lines = []
        labels = []
        axs[0].plot([], [], " ", color="k", label="$\\mu_R=" + "\\mu_F$")
        axs[1].plot(
            [], [], " ", color="k", label="$\\mu_R=" + str(np.max(mr)) + "\\mu_0$"
        )
        axs[2].plot(
            [], [], " ", color="k", label="$\\mu_F=" + str(np.min(mf)) + "\\mu_0$"
        )
        axs[3].plot(
            [], [], " ", color="k", label="$\\mu_R=" + str(np.min(mr)) + "\\mu_0$"
        )
        axs[4].plot(
            [], [], " ", color="k", label="$\\mu_F=" + str(np.max(mf)) + "\\mu_0$"
        )

    for v in vl:
        mv = dict_list[v] * yscale
        if seven_point_band:
            # mask = (mf/mr < 4.) & (mf/mr > 1./4.) & (mf <= 2) & (mf >= 1./2.) & (mr <= 2) & (mr >= 1./2.)
            mask = (
                ((mf == 2.0) & (mr == 2.0))
                | ((mf == 2.0) & (mr == 1.0))
                | ((mf == 1.0) & (mr == 1.0))
                | ((mf == 1.0) & (mr == 2.0))
                | ((mf == 1 / 2.0) & (mr == 1 / 2.0))
                | ((mf == 1 / 2.0) & (mr == 1.0))
                | ((mf == 1.0) & (mr == 1 / 2.0))
            )
            mvmax = np.max(splot.unv(mv[mask]))
            mvmin = np.min(splot.unv(mv[mask]))

        mask = mf == mr
        l = err_plt(axs[0], mf[mask], mv[mask], error=error)
        # l, _, _ = axs[0].errorbar(mf[mask], splot.unv(mv[mask]),
        #                          yerr=splot.usd(mv[mask]), capsize=5, label=v)
        if seven_point_band:
            axs[0].fill_between(
                mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
            )

        mask = mr == np.max(mr)
        l = err_plt(axs[1], mf[mask], mv[mask], error=error)
        lines.append(l)
        labels.append(
            "$\\sigma_{\\mathrm{"
            + v.replace("_PLUS_", "+").replace(" ", "\\ ")
            + "} }$"
        )
        # l, _, _ = axs[1].errorbar(mf[mask], splot.unv(mv[mask]),
        #                          yerr=splot.usd(mv[mask]), capsize=5)
        if seven_point_band:
            axs[1].fill_between(
                mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
            )

        mask = mf == np.min(mf)
        l = err_plt(axs[2], mr[mask], mv[mask], error=error)
        # l, _, _ = axs[2].errorbar(mr[mask], splot.unv(mv[mask]),
        #                          yerr=splot.usd(mv[mask]), capsize=5)
        if seven_point_band:
            axs[2].fill_between(
                mr[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
            )

        mask = mr == np.min(mr)
        l = err_plt(axs[3], mf[mask], mv[mask], error=error)
        # l, _, _ = axs[3].errorbar(mf[mask], splot.unv(mv[mask]),
        #                          yerr=splot.usd(mv[mask]), capsize=5)
        if seven_point_band:
            axs[3].fill_between(
                mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
            )

        mask = mf == np.max(mf)
        l = err_plt(axs[4], mr[mask], mv[mask], error=error)
        # l, _, _ = axs[4].errorbar(mr[mask], splot.unv(mv[mask]),
        #                          yerr=splot.usd(mv[mask]), capsize=5)
        if seven_point_band:
            f = axs[4].fill_between(
                mr[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
            )
            lines.append(f)
            labels.append(
                "$\\Delta \\sigma_{\\mathrm{"
                + v.replace("NLO_PLUS_NLL", "NLO+NLL").replace(" ", "\\<space>")
                + "} }$"
            )

    axs[0].set_ylabel("$\\sigma$ [" + unit + "]" if yaxis is None else yaxis)

    axs[0].set_xscale("log")
    axs[0].set_xlim(np.min(mf), np.max(mf))
    axs[0].set_xlabel("$\\mu_{R,F}/\\mu_0$")

    # axs[1].plot(t, s2)
    axs[1].set_xscale("log")
    axs[1].set_xlim(np.max(mf), np.min(mf))
    axs[1].set_xticks([1.0])
    axs[1].xaxis.set_minor_formatter(NullFormatter())
    axs[1].set_xlabel("$\\mu_{F}/\\mu_0$")

    # axs[2].plot(t, s3)
    axs[2].set_xscale("log")
    axs[2].set_xlim(np.max(mf), np.min(mf))
    axs[2].set_xticks([1.0])
    axs[2].xaxis.set_minor_formatter(NullFormatter())
    axs[2].set_xlabel("$\\mu_{R}/\\mu_0$")

    # axs[3].plot(t, s3)
    axs[3].set_xscale("log")
    axs[3].set_xlim(np.min(mf), np.max(mf))
    axs[3].set_xticks([1.0])
    axs[3].xaxis.set_minor_formatter(NullFormatter())
    axs[3].set_xlabel("$\\mu_{F}/\\mu_0$")

    # axs[4].plot(t, s3)
    axs[4].set_xscale("log")
    axs[4].set_xlim(np.min(mf), np.max(mf))
    axs[4].set_xticks([1.0])
    axs[4].xaxis.set_minor_formatter(NullFormatter())
    axs[4].set_xlabel("$\\mu_{R}/\\mu_0$")

    axs[0].legend(handletextpad=-0.0, handlelength=0, fancybox=False)
    axs[1].legend(handletextpad=-0.0, handlelength=0, fancybox=False)
    axs[2].legend(handletextpad=-0.0, handlelength=0, fancybox=False)
    axs[3].legend(handletextpad=-0.0, handlelength=0, fancybox=False)
    axs[4].legend(handletextpad=-0.0, handlelength=0, fancybox=False)
    fig.legend()
    fig.legends = []
    fig.legend(handles=lines, labels=labels, loc="center right")
    plt.subplots_adjust(right=0.84)
    # plt.show()
    mpl.rcParams["axes.prop_cycle"] = cycle_safe


def central_scale_plot(
    dict_list, vl, cont=False, error=True, yscale=1.0, unit="pb", yaxis=None
):
    """Creates a scale variance plot with 3 panels (ystacked)."""
    global fig, axs
    if not cont:
        fig, axs = plt.subplots(3, 1, figsize=(12, 8), sharex=True)
        # Remove horizontal space between axes
        fig.subplots_adjust(hspace=0)

    mr = dict_list["mu_r"]
    mf = dict_list["mu_f"]

    for v in vl:
        mv = dict_list[v] * yscale

        mask = mf == mr
        l = err_plt(axs[0], mf[mask], mv[mask], label=v, error=error)

        mask = mf == 1.0
        l = err_plt(axs[1], mr[mask], mv[mask], error=error)

        mask = mr == 1.0
        l = err_plt(axs[2], mf[mask], mv[mask], error=error)

    if not cont:
        axs[0].plot([], [], " ", label="$\\mu_R=\\mu_F=\\mu$")
        axs[1].plot([], [], " ", label="$\\mu_F=\\mu_0$, $\\mu_R=\\mu$")
        axs[2].plot([], [], " ", label="$\\mu_R=\\mu_0$, $\\mu_F=\\mu$")

    axs[1].set_ylabel("$\\sigma$ [" + unit + "]" if yaxis is None else yaxis)

    axs[0].set_xscale("log")
    # axs[0].set_xlim(np.min(mf), np.max(mf))
    axs[2].set_xlabel("$\\mu/\\mu_0$")

    axs[0].legend()
    axs[1].legend()
    axs[2].legend()
    # plt.show()


def init_double_plot(
    figsize=(6, 8), sharex=True, sharey=False, gridspec_kw={"height_ratios": [3, 1]}
):
    """Initialze subplot for Ratio/K plots with another figure below."""
    fig, axs = plt.subplots(
        2, 1, figsize=figsize, sharex=sharex, sharey=sharey, gridspec_kw=gridspec_kw
    )
    # Remove horizontal space between axes
    fig.subplots_adjust(hspace=0)
    return fig, axs

1	from collections.abc import Iterable	1✔
2
3	import matplotlib as mpl	1✔
4	import matplotlib.cm as cm	1✔
5	import matplotlib.pyplot as plt	1✔
6	import numpy as np	1✔
7	import pyslha	1✔
8	from matplotlib import colors	1✔
9	from matplotlib.image import NonUniformImage	1✔
10	from matplotlib.ticker import NullFormatter	1✔
11	from scipy import integrate as auc	1✔
12	from scipy.interpolate import BSpline, make_interp_spline	1✔
13	from smpl import plot as splot	1✔
14
15	from ..input import Input, get_output_dir, replace_macros	1✔
16	from ..particles import get_name	1✔
17
18
19	def title(	1✔
20	i: Input,
21	axe=None,
22	scenario=None,
23	diff_L_R=None,
24	extra="",
25	cms_energy=True,
26	pdf_info=True,
27	id=False,
28	**kwargs
29	):
30	"""Sets the title on axis `axe`."""
31	if axe is None:	×
32	axe = plt.gca()	×
33	axe.set_title(	×
34	"$pp\\to"
35	+ get_name(i.particle1)
36	+ get_name(i.particle2)
37	+ "$"
38	+ (" at $\\sqrt{S} = " + str(i.energy / 1000) + "$ TeV" if cms_energy else "")
39	+ ((" for " + i.slha.split(".")[0]) if scenario is None else (" " +scenario))
40	+ (" with " + i.pdf_nlo if pdf_info else "")
41	+ " "
42	+ extra
43	+ ((" [" + i.id + "]") if id else "")
44	)
45
46
47	def energy_plot(	1✔
48	dict_list,
49	y,
50	yscale=1.0,
51	xaxis="E [GeV]",
52	yaxis="$\\sigma$ [pb]",
53	label=None,
54	**kwargs
55	):
56	"""Plot energy on the x-axis."""
57	plot(	×
58	dict_list,
59	"energy",
60	y,
61	label=label,
62	xaxis=xaxis,
63	yaxis=yaxis,
64	logy=True,
65	yscale=yscale,
66	**kwargs
67	)
68
69
70	def combined_mass_plot(dict_list, y, part, label=None, **kwargs):	1✔
71	if label == "":	×
72	kwargs["label"] = ""	×
73	mass_plot(	×
74	dict_list,
75	y + "_NOERR",
76	part,
77	next_color=False,
78	fmt=" ",
79	fill=False,
80	interpolate_label=label,
81	**kwargs
82	)
83	if "label" in kwargs:	×
84	del kwargs["label"]	×
85	mass_plot(	×
86	dict_list,
87	y + "_COMBINED",
88	part,
89	next_color=True,
90	fit_fmt=" ",
91	fill=True,
92	label="",
93	interpolate_lower_uncertainty=False,
94	**kwargs
95	)
96
97
98	def combined_plot(dict_list, x, y, label=None, **kwargs):	1✔
99	if label == "":	×
100	kwargs["label"] = ""	×
101	plot(	×
102	dict_list,
103	x,
104	y + "_NOERR",
105	next_color=False,
106	fmt=" ",
107	fill=False,
108	interpolate_label=label,
109	**kwargs
110	)
111	if "label" in kwargs:	×
112	del kwargs["label"]	×
113	plot(	×
114	dict_list,
115	x,
116	y + "_COMBINED",
117	next_color=True,
118	fit_fmt=" ",
119	fill=True,
120	label="",
121	interpolate_lower_uncertainty=False,
122	**kwargs
123	)
124
125
126	def mass_plot(	1✔
127	dict_list,
128	y,
129	part,
130	logy=True,
131	yaxis="$\\sigma$ [pb]",
132	yscale=1.0,
133	label=None,
134	xaxis=None,
135	**kwargs
136	):
137	if not "mass_" + str(part) in dict_list:	×
138	dict_list["mass_" + str(part)] = get_mass(dict_list, abs(part))	×
139	if xaxis is None:	×
140	xaxis = "$m_{" + get_name(part) + "}$ [GeV]"	×
141	plot(	×
142	dict_list,
143	"mass_" + str(part),
144	y,
145	label=label,
146	yaxis=yaxis,
147	xaxis=xaxis,
148	logy=logy,
149	yscale=yscale,
150	**kwargs
151	)
152
153
154	def mass_vplot(	1✔
155	dict_list,
156	y,
157	part,
158	logy=True,
159	yaxis="$\\sigma$ [pb]",
160	yscale=1.0,
161	label=None,
162	mask=None,
163	**kwargs
164	):
165	vplot(	×
166	get_mass(dict_list, part)[mask],
167	y[mask],
168	label=label,
169	xaxis="$m_{" + get_name(part) + "}$ [GeV]",
170	yaxis=yaxis,
171	logy=logy,
172	yscale=yscale,
173	mask=mask,
174	**kwargs
175	)
176
177
178	def get_mass(l: dict, iid: int):	1✔
179	"""
180	Get the mass of particle with id `iid` out of the list in the "slha" element in the dict.
181
182	Returns
183	:obj:`list` of float : masses of particles in each element of the dict list.
184
185	"""
186	ret = []	×
187	for s in l["slha"]:	×
188	d = pyslha.read(get_output_dir() + s)	×
189	ret.append(d.blocks["MASS"][abs(iid)])	×
190
191	return np.array(ret)	×
192
193
194	def plot(	1✔
195	dict_list,
196	x,
197	y,
198	label=None,
199	xaxis="M [GeV]",
200	yaxis="$\\sigma$ [pb]",
201	ratio=False,
202	K=False,
203	K_plus_1=False,
204	logy=True,
205	yscale=1.0,
206	mask=None,
207	**kwargs
208	) -> None:
209	"""
210	Creates a plot based on the entries `x`and `y` in `dict_list`.
211
212	Examples
213
214	.. plot::
215	:include-source:
216
217	>>> import urllib.request
218	>>> import hepi
219	>>> dl = hepi.load(urllib.request.urlopen(
220	... "https://raw.githubusercontent.com/fuenfundachtzig/xsec/master/json/pp13_hino_NLO%2BNLL.json"
221	... ))
222	>>> hepi.plot(dl,"N1","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{\\\\chi}_1^0}$ [GeV]")
223	"""
224	if isinstance(y, Iterable) and not isinstance(y, str):	×
225	for yi in y:	×
226	plot(	×
227	dict_list,
228	x,
229	yi,
230	label,
231	xaxis,
232	yaxis,
233	ratio,
234	K,
235	K_plus_1,
236	logy,
237	yscale,
238	mask,
239	**kwargs
240	)
241	return	×
242	# TODO use kwargs
243	if label is None:	×
244	label = y.replace("_PLUS_", "+").replace("_OVER_", "/")	×
245	if label == "":	×
246	label = None	×
247
248	vx = dict_list[x].to_numpy()[mask]	×
249	vy = dict_list[y].to_numpy()[mask]	×
250
251	if K:	×
252	yaxis = "$K$"	×
253	yscale = 1.0	×
254	vy = vy / splot.unv(dict_list["LO"].to_numpy()[mask])	×
255	if K_plus_1:	×
256	vy = vy + vy / vy	×
257	if ratio:	×
258	yaxis = "Ratio"	×
259	yscale = 1.0	×
260	vy = vy / splot.unv(vy)	×
261
262	vplot(vx, vy, label, xaxis, yaxis, logy, yscale, mask=mask, **kwargs)	×
263
264
265	def index_open(var, idx):	1✔
266	if len(idx) == 1:	×
267	return var[idx]	×
268	return index_open(var[idx[0]], idx[1:])	×
269
270
271	def slha_data(li, index_list):	1✔
272	vx = []	×
273	for l in li:	×
274	b = pyslha.read(get_output_dir() + l.slha)	×
275	vx.append(index_open(b.blocks, index_list))	×
276	return np.array(vx)	×
277
278
279	def slha_plot(li, x, y, **kwargs):	1✔
280	vx = slha_data(li, x)	×
281	vy = slha_data(li, y)	×
282
283	vplot(np.array(vx), np.array(vy), **kwargs)	×
284
285
286	def vplot(	1✔
287	x,
288	y,
289	label=None,
290	xaxis="E [GeV]",
291	yaxis="$\\sigma$ [pb]",
292	logy=True,
293	yscale=1.0,
294	interpolate=True,
295	plot_data=True,
296	data_color=None,
297	mask=-1,
298	fill=False,
299	data_fmt=".",
300	fmt="-",
301	print_area=False,
302	sort=True,
303	**kwargs
304	):
305	"""
306	Creates a plot based on the values in `x`and `y`.
307
308	"""
309	color = data_color	×
310	if label is None:	×
311	# label = "??"
312	pass	×
313	if mask is None:	×
314	x = x[0]	×
315	y = y[0]	×
316	if sort:	×
317	permute = x.argsort(kind="stable")	×
318	vx = x[permute]	×
319	vy = y[permute]	×
320	else:
321	vx = x	×
322	vy = y	×
323
324	if data_color is None and "next_color" not in kwargs:	×
325	if "axes" in kwargs and kwargs["axes"] is not None:	×
326	(bl,) = kwargs["axes"].plot([], [])	×
327	else:
328	(bl,) = plt.gca().plot([], [])	×
329	color = bl.get_color()	×
330	iii = splot.data(	×
331	vx,
332	vy * yscale,
333	label=label,
334	xaxis=xaxis,
335	yaxis=yaxis,
336	logy=logy,
337	data_color=color,
338	fit_color=color,
339	also_data=plot_data,
340	interpolate=interpolate,
341	fmt=data_fmt,
342	sigmas=0 if not fill else 1,
343	**kwargs
344	)
345	if iii is not None:	×
346	# ii = iii[0]
347	ix = iii[1]	×
348	iy = iii[2]	×
349	if print_area:	×
350	print(	×
351	"computed AUC using scipy.integrate.simpson: {}".format(auc.simpson(iy, ix))
352	)
353	if (np.any(np.less(vy, 0)) or (interpolate and np.any(np.less(iy, 0)))) and logy:	×
354	splot.data(	×
355	vx,
356	-vy * yscale,
357	label="-" + label,
358	xaxis=xaxis,
359	yaxis=yaxis,
360	logy=logy,
361	data_color=color,
362	fit_color=color,
363	also_data=plot_data,
364	interpolate=interpolate,
365	interpolate_fmt=(0, (3, 1, 3, 1, 1, 1)),
366	fmt=data_fmt,
367	sigmas=0 if not fill else 1,
368	**kwargs
369	)
370
371
372	def mass_mapplot(	1✔
373	dict_list,
374	part1,
375	part2,
376	z,
377	logz=True,
378	zaxis="$\\sigma$ [pb]",
379	zscale=1.0,
380	label=None,
381	):
382	mapplot(	×
383	dict_list,
384	"mass_" + str(part1),
385	"mass_" + str(part2),
386	z,
387	xaxis="$M_{" + get_name(part1) + "}$ [GeV]",
388	yaxis="$M_{" + get_name(part2) + "}$ [GeV]",
389	zaxis=zaxis,
390	logz=logz,
391	zscale=zscale,
392	)
393
394
395	def mapplot(dict_list, x, y, z, xaxis=None, yaxis=None, zaxis=None, **kwargs):	1✔
396	"""
397
398	Examples
399
400	.. plot::
401	:include-source:
402
403	>>> import urllib.request
404	>>> import hepi
405
406	>>> dl = hepi.load(urllib.request.urlopen(
407	... "https://raw.githubusercontent.com/APN-Pucky/xsec/master/json/pp13_SGmodel_GGxsec_NLO%2BNLL.json"
408	... ),dimensions=2)
409	>>> hepi.mapplot(dl,"gl","sq","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{g}}$ [GeV]",yaxis="$m_{\\\\tilde{q}}$ [GeV]" , zaxis="$\\\\sigma_{\\\\mathrm{NLO+NLL}}$ [pb]")
410	"""
411	if xaxis is None:	×
412	xaxis = x	×
413	if yaxis is None:	×
414	yaxis = y	×
415	if zaxis is None:	×
416	zaxis = replace_macros(z)	×
417	vx = dict_list[x]	×
418	vy = dict_list[y]	×
419	vz = dict_list[z]	×
420	splot.plot2d(	×
421	vx, vy, vz, style="pcolormesh", xaxis=xaxis, yaxis=yaxis, zaxis=zaxis, **kwargs
422	)
423
424
425	map_vplot = lambda a, da: splot.plot2d(a, style="pcolormesh", **da)	1✔
426	scatter_vplot = lambda a, da: splot.plot2d(a, style="scatter", **da)	1✔
427
428
429	def scatterplot(dict_list, x, y, z, xaxis=None, yaxis=None, zaxis=None, **kwargs):	1✔
430	"""
431	Scatter map 2d.
432	Central color is the central value, while the inner and outer ring are lower and upper bounds of the uncertainty interval.
433
434	Examples
435
436	.. plot::
437	:include-source:
438
439	>>> import urllib.request
440	>>> import hepi
441	>>> dl = hepi.load(urllib.request.urlopen(
442	... "https://raw.githubusercontent.com/APN-Pucky/xsec/master/json/pp13_hinosplit_N2N1_NLO%2BNLL.json"
443	... ),dimensions=2)
444	>>> hepi.scatterplot(dl,"N1","N2","NLO_PLUS_NLL_COMBINED",xaxis="$m_{\\\\tilde{\\\\chi}_1^0}$ [GeV]",yaxis="$m_{\\\\tilde{\\\\chi}_2^0}$ [GeV]" , zaxis="$\\\\sigma_{\\\\mathrm{NLO+NLL}}$ [pb]")
445
446	"""
447	if xaxis is None:	×
448	xaxis = x	×
449	if yaxis is None:	×
450	yaxis = y	×
451	if zaxis is None:	×
452	zaxis = replace_macros(z)	×
453	vx = dict_list[x]	×
454	vy = dict_list[y]	×
455	vz = dict_list[z]	×
456	splot.plot2d(	×
457	vx, vy, vz, style="scatter", xaxis=xaxis, yaxis=yaxis, zaxis=zaxis, **kwargs
458	)
459
460
461	fig = None	1✔
462	axs = None	1✔
463	lines = []	1✔
464	labels = []	1✔
465
466
467	def err_plt(axes, x, y, label=None, error=False):	1✔
468	v = label	×
469	ind = np.argsort(splot.unv(x), kind="stable")	×
470	if error:	×
471	l, _, _ = axes.errorbar(	×
472	x.to_numpy()[ind], splot.unv(y)[ind], yerr=splot.usd(y), capsize=5, label=v
473	)
474	return l	×
475	else:
476	l = axes.plot(x.to_numpy()[ind], splot.unv(y)[ind], label=v)	×
477	return l[0]	×
478
479
480	def scale_plot(	1✔
481	dict_list,
482	vl,
483	seven_point_band=False,
484	cont=False,
485	error=True,
486	li=None,
487	plehn_color=False,
488	yscale=1.0,
489	unit="pb",
490	yaxis=None,
491	**kwargs
492	):
493	"""Creates a scale variance plot with 5 panels (xline)."""
494	global fig, axs, lines, labels
495	cycle_safe = mpl.rcParams["axes.prop_cycle"]	×
496	if plehn_color:	×
497	mpl.rcParams["axes.prop_cycle"] = mpl.cycler(color=["b", "r", "g"])	×
498	if not cont:	×
499	fig, axs = plt.subplots(1, 5, figsize=(12, 3), sharey=True)	×
500	# Remove horizontal space between axes
501	fig.subplots_adjust(wspace=0)	×
502	if li is not None:	×
503	title(axs[2], li[0], **kwargs)	×
504
505	mr = dict_list["mu_r"]	×
506	mf = dict_list["mu_f"]	×
507
508	if not cont:	×
509	lines = []	×
510	labels = []	×
511	axs[0].plot([], [], " ", color="k", label="$\\mu_R=" + "\\mu_F$")	×
512	axs[1].plot(	×
513	[], [], " ", color="k", label="$\\mu_R=" + str(np.max(mr)) + "\\mu_0$"
514	)
515	axs[2].plot(	×
516	[], [], " ", color="k", label="$\\mu_F=" + str(np.min(mf)) + "\\mu_0$"
517	)
518	axs[3].plot(	×
519	[], [], " ", color="k", label="$\\mu_R=" + str(np.min(mr)) + "\\mu_0$"
520	)
521	axs[4].plot(	×
522	[], [], " ", color="k", label="$\\mu_F=" + str(np.max(mf)) + "\\mu_0$"
523	)
524
525	for v in vl:	×
526	mv = dict_list[v] * yscale	×
527	if seven_point_band:	×
528	# mask = (mf/mr < 4.) & (mf/mr > 1./4.) & (mf <= 2) & (mf >= 1./2.) & (mr <= 2) & (mr >= 1./2.)
529	mask = (	×
530	((mf == 2.0) & (mr == 2.0))
531	\| ((mf == 2.0) & (mr == 1.0))
532	\| ((mf == 1.0) & (mr == 1.0))
533	\| ((mf == 1.0) & (mr == 2.0))
534	\| ((mf == 1 / 2.0) & (mr == 1 / 2.0))
535	\| ((mf == 1 / 2.0) & (mr == 1.0))
536	\| ((mf == 1.0) & (mr == 1 / 2.0))
537	)
538	mvmax = np.max(splot.unv(mv[mask]))	×
539	mvmin = np.min(splot.unv(mv[mask]))	×
540
541	mask = mf == mr	×
542	l = err_plt(axs[0], mf[mask], mv[mask], error=error)	×
543	# l, _, _ = axs[0].errorbar(mf[mask], splot.unv(mv[mask]),
544	# yerr=splot.usd(mv[mask]), capsize=5, label=v)
545	if seven_point_band:	×
546	axs[0].fill_between(	×
547	mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
548	)
549
550	mask = mr == np.max(mr)	×
551	l = err_plt(axs[1], mf[mask], mv[mask], error=error)	×
552	lines.append(l)	×
553	labels.append(	×
554	"$\\sigma_{\\mathrm{"
555	+ v.replace("_PLUS_", "+").replace(" ", "\\ ")
556	+ "} }$"
557	)
558	# l, _, _ = axs[1].errorbar(mf[mask], splot.unv(mv[mask]),
559	# yerr=splot.usd(mv[mask]), capsize=5)
560	if seven_point_band:	×
561	axs[1].fill_between(	×
562	mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
563	)
564
565	mask = mf == np.min(mf)	×
566	l = err_plt(axs[2], mr[mask], mv[mask], error=error)	×
567	# l, _, _ = axs[2].errorbar(mr[mask], splot.unv(mv[mask]),
568	# yerr=splot.usd(mv[mask]), capsize=5)
569	if seven_point_band:	×
570	axs[2].fill_between(	×
571	mr[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
572	)
573
574	mask = mr == np.min(mr)	×
575	l = err_plt(axs[3], mf[mask], mv[mask], error=error)	×
576	# l, _, _ = axs[3].errorbar(mf[mask], splot.unv(mv[mask]),
577	# yerr=splot.usd(mv[mask]), capsize=5)
578	if seven_point_band:	×
579	axs[3].fill_between(	×
580	mf[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
581	)
582
583	mask = mf == np.max(mf)	×
584	l = err_plt(axs[4], mr[mask], mv[mask], error=error)	×
585	# l, _, _ = axs[4].errorbar(mr[mask], splot.unv(mv[mask]),
586	# yerr=splot.usd(mv[mask]), capsize=5)
587	if seven_point_band:	×
588	f = axs[4].fill_between(	×
589	mr[mask], mvmax, mvmin, facecolor=l.get_color(), alpha=0.3
590	)
591	lines.append(f)	×
592	labels.append(	×
593	"$\\Delta \\sigma_{\\mathrm{"
594	+ v.replace("NLO_PLUS_NLL", "NLO+NLL").replace(" ", "\\<space>")
595	+ "} }$"
596	)
597
598	axs[0].set_ylabel("$\\sigma$ [" + unit + "]" if yaxis is None else yaxis)	×
599
600	axs[0].set_xscale("log")	×
601	axs[0].set_xlim(np.min(mf), np.max(mf))	×
602	axs[0].set_xlabel("$\\mu_{R,F}/\\mu_0$")	×
603
604	# axs[1].plot(t, s2)
605	axs[1].set_xscale("log")	×
606	axs[1].set_xlim(np.max(mf), np.min(mf))	×
607	axs[1].set_xticks([1.0])	×
608	axs[1].xaxis.set_minor_formatter(NullFormatter())	×
609	axs[1].set_xlabel("$\\mu_{F}/\\mu_0$")	×
610
611	# axs[2].plot(t, s3)
612	axs[2].set_xscale("log")	×
613	axs[2].set_xlim(np.max(mf), np.min(mf))	×
614	axs[2].set_xticks([1.0])	×
615	axs[2].xaxis.set_minor_formatter(NullFormatter())	×
616	axs[2].set_xlabel("$\\mu_{R}/\\mu_0$")	×
617
618	# axs[3].plot(t, s3)
619	axs[3].set_xscale("log")	×
620	axs[3].set_xlim(np.min(mf), np.max(mf))	×
621	axs[3].set_xticks([1.0])	×
622	axs[3].xaxis.set_minor_formatter(NullFormatter())	×
623	axs[3].set_xlabel("$\\mu_{F}/\\mu_0$")	×
624
625	# axs[4].plot(t, s3)
626	axs[4].set_xscale("log")	×
627	axs[4].set_xlim(np.min(mf), np.max(mf))	×
628	axs[4].set_xticks([1.0])	×
629	axs[4].xaxis.set_minor_formatter(NullFormatter())	×
630	axs[4].set_xlabel("$\\mu_{R}/\\mu_0$")	×
631
632	axs[0].legend(handletextpad=-0.0, handlelength=0, fancybox=False)	×
633	axs[1].legend(handletextpad=-0.0, handlelength=0, fancybox=False)	×
634	axs[2].legend(handletextpad=-0.0, handlelength=0, fancybox=False)	×
635	axs[3].legend(handletextpad=-0.0, handlelength=0, fancybox=False)	×
636	axs[4].legend(handletextpad=-0.0, handlelength=0, fancybox=False)	×
637	fig.legend()	×
638	fig.legends = []	×
639	fig.legend(handles=lines, labels=labels, loc="center right")	×
640	plt.subplots_adjust(right=0.84)	×
641	# plt.show()
642	mpl.rcParams["axes.prop_cycle"] = cycle_safe	×
643
644
645	def central_scale_plot(	1✔
646	dict_list, vl, cont=False, error=True, yscale=1.0, unit="pb", yaxis=None
647	):
648	"""Creates a scale variance plot with 3 panels (ystacked)."""
649	global fig, axs
650	if not cont:	×
651	fig, axs = plt.subplots(3, 1, figsize=(12, 8), sharex=True)	×
652	# Remove horizontal space between axes
653	fig.subplots_adjust(hspace=0)	×
654
655	mr = dict_list["mu_r"]	×
656	mf = dict_list["mu_f"]	×
657
658	for v in vl:	×
659	mv = dict_list[v] * yscale	×
660
661	mask = mf == mr	×
662	l = err_plt(axs[0], mf[mask], mv[mask], label=v, error=error)	×
663
664	mask = mf == 1.0	×
665	l = err_plt(axs[1], mr[mask], mv[mask], error=error)	×
666
667	mask = mr == 1.0	×
668	l = err_plt(axs[2], mf[mask], mv[mask], error=error)	×
669
670	if not cont:	×
671	axs[0].plot([], [], " ", label="$\\mu_R=\\mu_F=\\mu$")	×
672	axs[1].plot([], [], " ", label="$\\mu_F=\\mu_0$, $\\mu_R=\\mu$")	×
673	axs[2].plot([], [], " ", label="$\\mu_R=\\mu_0$, $\\mu_F=\\mu$")	×
674
675	axs[1].set_ylabel("$\\sigma$ [" + unit + "]" if yaxis is None else yaxis)	×
676
677	axs[0].set_xscale("log")	×
678	# axs[0].set_xlim(np.min(mf), np.max(mf))
679	axs[2].set_xlabel("$\\mu/\\mu_0$")	×
680
681	axs[0].legend()	×
682	axs[1].legend()	×
683	axs[2].legend()	×
684	# plt.show()
685
686
687	def init_double_plot(	1✔
688	figsize=(6, 8), sharex=True, sharey=False, gridspec_kw={"height_ratios": [3, 1]}
689	):
690	"""Initialze subplot for Ratio/K plots with another figure below."""
691	fig, axs = plt.subplots(	×
692	2, 1, figsize=figsize, sharex=sharex, sharey=sharey, gridspec_kw=gridspec_kw
693	)
694	# Remove horizontal space between axes
695	fig.subplots_adjust(hspace=0)	×
696	return fig, axs	×

APN-Pucky / HEPi / 6543527074

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous