From 445d7f0bfe04ea7babb102bec6f37c16203846f7 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Tue, 29 Apr 2025 18:28:40 +0200 Subject: [PATCH 1/9] Add preliminary version of 3d plot --- pdfplotter/pdf_set_nuclear.py | 308 +++++++++++++++++++++++++++++++++- 1 file changed, 307 insertions(+), 1 deletion(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index b042d26..32590e3 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -1,13 +1,23 @@ from __future__ import annotations from itertools import zip_longest -from typing import Sequence +from math import log +from typing import Sequence, Literal, Any +#from mpl_toolkits.mplot3d import Axes3D +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import matplotlib.ticker as mticker +import matplotlib.cm as cm import numpy as np import numpy.typing as npt import pandas as pd +import matplotlib.pyplot as plt +import sympy as sp from pdfplotter.pdf_set import PDFSet +from pdfplotter.util import update_kwargs +from pdfplotter.util import log_tick_formatter +from pdfplotter import util class NuclearPDFSet(PDFSet): @@ -153,3 +163,299 @@ def get( raise ValueError(f"Multiple PDFSets found for Z = {Z}") else: return pdf_set.iloc[0]["pdf_set"] + + def plot_A_dep_3d( + self, + ax: plt.Axes | npt.NDArray[plt.Axes], # pyright: ignore[reportInvalidTypeForm] + A: int | float | list[int | float], + observables: ( + sp.Basic + | npt.NDArray[sp.Basic] # pyright: ignore[reportInvalidTypeForm] + | list[sp.Basic] + ), + Q: float | None = None, + Q2: float | None = None, + colors: list = [], + logA: bool = True, + plot_uncertainty: bool = True, + plot_ratio: bool = False, + pdf_label: Literal["ylabel", "annotate"] = "annotate", + A_label: Literal["legend", "ticks"] = "ticks", + proj_type: Literal["ortho", "persp"] = "ortho", + view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), + kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_uncertainty_edges: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_title: dict[str, Any] = {}, + kwargs_notation: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_ylabel: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_xlabel: dict[str, Any] = {}, + kwargs_zlabel: dict[str, Any] = {}, + kwargs_legend: dict[str, Any] = {}, + ): + + my_sets = self.pdf_sets + x = self.get(A=my_sets["A"][0]).x_values + + + if Q is None and Q2 is None: + raise ValueError("Please pass either `Q` or `Q2`") + + elif Q is not None and Q2 is not None: + raise ValueError("Please pass either `Q` or `Q2`, not both") + + elif Q is not None: + if Q not in self.get(A=my_sets["A"][0]).Q_values and Q not in np.sqrt( + np.array(self.get(A=my_sets["A"][0]).Q2_values) + ): + raise ValueError( + f"Chosen Q value {Q} was not used for defining nuclear pdf set. \n Please choose Q that was used in initialization" + ) + else: + if ( + Q2 not in self.get(A=my_sets["A"][0]).Q2_values + and Q2 not in np.array(self.get(A=my_sets["A"][0]).Q_values) ** 2 + ): + raise ValueError( + f"Chosen Q2 value {Q2} was not used for defining nuclear pdf set. \n Please choose Q2 that was used in initialization" + ) + + if not isinstance(A, list): + A = [A] + + if isinstance(observables, np.ndarray): + observables = list(observables.flatten()) + + if not isinstance(observables, list): + observables = [observables] + + if not isinstance(ax, np.ndarray): + ax = np.array([ax]) + + if colors == []: + cmap = cm.get_cmap("tab10", lut=len(A)) + colors = [cmap(i) for i in range(len(A))] + + for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): + + ax_i.set_proj_type(proj_type) + ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) + for j, (A_j, col_j) in enumerate(zip(A, colors)): + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + kwargs_default = { + "color": col_j, + "label": f"A={A_j}", + "linewidth": 1.5, + } + if not isinstance(kwargs_theory, list): + kwargs = update_kwargs( + kwargs_default, + kwargs_theory, + ) + else: + kwargs = update_kwargs( + kwargs_default, + kwargs_theory, + i=j, + ) + if logA: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + if plot_uncertainty: + kwargs_uncertainty_default = { + "color": col_j, + "alpha": 0.3, + } + if not isinstance(kwargs_uncertainty, list): + kwargs = update_kwargs( + kwargs_uncertainty_default, + kwargs_uncertainty, + ) + else: + kwargs = update_kwargs( + kwargs_uncertainty_default, + kwargs_uncertainty, + i=j, + ) + + vertices = [] + z_lower = np.array(z_lower) + z_upper = np.array(z_upper) + if not logA: + for xi, ai, zl, zu in zip( + np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper + ): + vertices.append([xi, ai, zl]) + + for xi, ai, zl, zu in reversed( + list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) + ): + vertices.append([xi, ai, zu]) + else: + for xi, ai, zl, zu in zip( + np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper + ): + vertices.append([xi, ai, zl]) + + for xi, ai, zl, zu in reversed( + list(zip(np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper)) + ): + vertices.append([xi, ai, zu]) + poly = Poly3DCollection([vertices], **kwargs) + ax_i.add_collection3d(poly) + + kwargs_uncertainty_edges_default = { + "color": col_j, + "alpha": 1, + } + if not isinstance(kwargs_uncertainty_edges, list): + kwargs = update_kwargs( + kwargs_uncertainty_edges_default, + kwargs_uncertainty_edges, + ) + else: + kwargs = update_kwargs( + kwargs_uncertainty_edges_default, + kwargs_uncertainty_edges, + i=j, + ) + if not logA: + ax_i.plot(np.log10(x), len(x) * [A_j], z_upper, **kwargs) + ax_i.plot(np.log10(x), len(x) * [A_j], z_lower, **kwargs) + else: + ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) + ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) + if pdf_label == "annotate": + kwargs_notation_default = { + "fontsize": 12, + "xy": (0.97, 0.96), + "xycoords": "axes fraction", + "va": "top", + "ha": "right", + "bbox": dict( + facecolor=(1, 1, 1), + edgecolor=(0.8, 0.8, 0.8), + lw=0.9, + boxstyle="round, pad=0.2", + ), + } + if not isinstance(kwargs_notation, list): + kwargs_n = update_kwargs( + kwargs_notation_default, + kwargs_notation, + ) + else: + kwargs_n = update_kwargs( + kwargs_notation_default, + kwargs_notation, + i=i, + ) + ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) + + if pdf_label == "ylabel": + kwargs_ylabel_default = { + "fontsize": 14, + "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", + #"labelpad":-200 + } + if not isinstance(kwargs_ylabel, list): + kwargs = update_kwargs( + kwargs_ylabel_default, + kwargs_ylabel, + ) + else: + kwargs = update_kwargs( + kwargs_ylabel_default, + kwargs_ylabel, + i=i, + ) + ax_i.set_zlabel(**kwargs) + + else: + kwargs_notation_default = { + "fontsize": 12, + "xy": (0.47, 0.96), + "xycoords": "axes fraction", + "va": "top", + "ha": "right", + "bbox": dict( + facecolor=(1, 1, 1), + edgecolor=(0.8, 0.8, 0.8), + lw=0.9, + boxstyle="round, pad=0.2", + ), + } + kwargs_n = update_kwargs(kwargs_notation_default, kwargs_notation, i=i) + + ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) + ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) + if A_label == "ticks": + if logA: + ax_i.set_yticks(np.log10(A),A) + else: + ax_i.set_yticks(A,A) + kwargs_zlabel_default = { + "fontsize": 14, + "ylabel": f"$A$", + + } + kwargs = update_kwargs( + kwargs_zlabel_default, + kwargs_zlabel, + ) + + ax_i.set_ylabel(**kwargs) + else: + ax_i.set_yticks([]) + kwargs_legend_default = { + "fontsize": 12, + "bbox_to_anchor": (0.95, 0.95), + "frameon": False, + } + kwargs = update_kwargs( + kwargs_legend_default, + kwargs_legend, + ) + ax_i.legend() + kwargs_zlabel_default = { + "fontsize": 14, + "ylabel": f"$A$", + "labelpad":-10 + } + kwargs = update_kwargs( + kwargs_zlabel_default, + kwargs_zlabel, + ) + ax_i.set_ylabel(**kwargs) + ax_i.xaxis.pane.fill=False + ax_i.yaxis.pane.fill=False + ax_i.zaxis.pane.fill=False + ax_i.xaxis.pane.set_edgecolor("white") + ax_i.yaxis.pane.set_edgecolor("white") + ax_i.zaxis.pane.set_edgecolor("white") + + ax_i.zaxis._axinfo["juggled"]=(1,2,0) + + kwargs_xlabel_default = { + "fontsize": 14, + "xlabel": f"$x$", + + } + kwargs = update_kwargs( + kwargs_xlabel_default, + kwargs_xlabel, + ) + ax_i.set_xlabel(**kwargs) \ No newline at end of file From ebdc04d8c300c324a3a3afe3ca4306e39488ce71 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Wed, 30 Apr 2025 10:12:05 +0200 Subject: [PATCH 2/9] Add log_tick_formatter to util.py. Add possibility to plot ratios. --- pdfplotter/pdf_set_nuclear.py | 129 +++++++++++++++++++++++++++------- pdfplotter/util.py | 3 + 2 files changed, 107 insertions(+), 25 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 32590e3..7011a95 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -175,13 +175,14 @@ def plot_A_dep_3d( ), Q: float | None = None, Q2: float | None = None, + x_lines: float | list[float] | None = None, colors: list = [], logA: bool = True, plot_uncertainty: bool = True, plot_ratio: bool = False, pdf_label: Literal["ylabel", "annotate"] = "annotate", A_label: Literal["legend", "ticks"] = "ticks", - proj_type: Literal["ortho", "persp"] = "ortho", + proj_type: Literal["ortho", "persp"] = "persp", view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, @@ -192,6 +193,7 @@ def plot_A_dep_3d( kwargs_xlabel: dict[str, Any] = {}, kwargs_zlabel: dict[str, Any] = {}, kwargs_legend: dict[str, Any] = {}, + kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {}, ): my_sets = self.pdf_sets @@ -223,6 +225,11 @@ def plot_A_dep_3d( if not isinstance(A, list): A = [A] + if 1 not in A and plot_ratio: + raise ValueError( + "Please pass A=1 if you want to plot the ratio to Proton." + ) + if isinstance(observables, np.ndarray): observables = list(observables.flatten()) @@ -238,12 +245,17 @@ def plot_A_dep_3d( for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): - ax_i.set_proj_type(proj_type) - ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) for j, (A_j, col_j) in enumerate(zip(A, colors)): - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) + if not plot_ratio: + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + else: + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + z_lower = z_lower / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) + z_upper = z_upper / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) kwargs_default = { "color": col_j, "label": f"A={A_j}", @@ -261,19 +273,37 @@ def plot_A_dep_3d( i=j, ) if logA: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) + if plot_ratio: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) + / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) else: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) + if plot_ratio: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) + / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) if plot_uncertainty: kwargs_uncertainty_default = { "color": col_j, @@ -295,15 +325,16 @@ def plot_A_dep_3d( z_lower = np.array(z_lower) z_upper = np.array(z_upper) if not logA: + for xi, ai, zl, zu in zip( np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper ): vertices.append([xi, ai, zl]) - for xi, ai, zl, zu in reversed( list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) ): vertices.append([xi, ai, zu]) + else: for xi, ai, zl, zu in zip( np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper @@ -338,6 +369,50 @@ def plot_A_dep_3d( else: ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) + + centrals={} + if x_lines is not None: + if not isinstance(x_lines, list): + x_lines = [x_lines] + for k,x_line in enumerate(x_lines): + if x_line not in x: + raise ValueError( + f"Chosen x value {x_line} was not used for defining nuclear pdf set. \n Please choose x that was used in initialization" + ) + kwargs_xlines_default = { + "color": "black", + "linestyle": "--", + "linewidth": 1.5, + } + if not isinstance(kwargs_xlines, list): + kwargs = update_kwargs( + kwargs_xlines_default, + kwargs_xlines, + ) + else: + kwargs = update_kwargs( + kwargs_xlines_default, + kwargs_xlines, + i=k, + ) + for a in A: + if not plot_ratio: + if x_line not in centrals.keys(): + centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + else: + centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + else: + if x_line not in centrals.keys(): + centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + else: + centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + if logA: + ax_i.plot( + np.ones(len(A))*np.log10(x_line), np.log10(A), centrals[x_line],**kwargs) + else: + ax_i.plot( + np.ones(len(A))*np.log10(x_line), A, centrals[x_line],**kwargs) + if pdf_label == "annotate": kwargs_notation_default = { "fontsize": 12, @@ -369,7 +444,6 @@ def plot_A_dep_3d( kwargs_ylabel_default = { "fontsize": 14, "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", - #"labelpad":-200 } if not isinstance(kwargs_ylabel, list): kwargs = update_kwargs( @@ -433,7 +507,8 @@ def plot_A_dep_3d( kwargs_zlabel_default = { "fontsize": 14, "ylabel": f"$A$", - "labelpad":-10 + "labelpad":-10, + "linespacing": -4 } kwargs = update_kwargs( kwargs_zlabel_default, @@ -443,9 +518,9 @@ def plot_A_dep_3d( ax_i.xaxis.pane.fill=False ax_i.yaxis.pane.fill=False ax_i.zaxis.pane.fill=False - ax_i.xaxis.pane.set_edgecolor("white") - ax_i.yaxis.pane.set_edgecolor("white") - ax_i.zaxis.pane.set_edgecolor("white") + ax_i.xaxis.pane.set_edgecolor("w") + ax_i.yaxis.pane.set_edgecolor("w") + ax_i.zaxis.pane.set_edgecolor("w") ax_i.zaxis._axinfo["juggled"]=(1,2,0) @@ -458,4 +533,8 @@ def plot_A_dep_3d( kwargs_xlabel_default, kwargs_xlabel, ) - ax_i.set_xlabel(**kwargs) \ No newline at end of file + ax_i.set_xlabel(**kwargs) + ax_i.set_xlim(np.log10(x[0])-np.log10(x[0])/100, np.log10(x[-1])) + ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 + ax_i.set_proj_type(proj_type) + ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file diff --git a/pdfplotter/util.py b/pdfplotter/util.py index 16903fc..8903026 100644 --- a/pdfplotter/util.py +++ b/pdfplotter/util.py @@ -201,3 +201,6 @@ def update_kwargs( return kwargs else: raise ValueError("kwargs_user must be a dict or a list") + +def log_tick_formatter(val,pos=None): + return r"$10^{{{:.0f}}}$".format(val) \ No newline at end of file From 3ebe02b37eb9187404729541596b84bf6e3b63d4 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Tue, 6 May 2025 10:03:06 +0200 Subject: [PATCH 3/9] Remove y and x axis offset. Add vertical grid lines for case: A_label:yticks --- pdfplotter/pdf_set_nuclear.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 7011a95..4c7c021 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -492,18 +492,20 @@ def plot_A_dep_3d( ) ax_i.set_ylabel(**kwargs) + ax_i.set_xlim(np.log10(x[0])*0.98) else: ax_i.set_yticks([]) - kwargs_legend_default = { - "fontsize": 12, - "bbox_to_anchor": (0.95, 0.95), - "frameon": False, - } - kwargs = update_kwargs( - kwargs_legend_default, - kwargs_legend, - ) - ax_i.legend() + if i==len(observables)-1: + kwargs_legend_default = { + "fontsize": 12, + "bbox_to_anchor": (0.95, 0.95), + "frameon": False, + } + kwargs = update_kwargs( + kwargs_legend_default, + kwargs_legend, + ) + ax_i.legend() kwargs_zlabel_default = { "fontsize": 14, "ylabel": f"$A$", @@ -534,7 +536,8 @@ def plot_A_dep_3d( kwargs_xlabel, ) ax_i.set_xlabel(**kwargs) - ax_i.set_xlim(np.log10(x[0])-np.log10(x[0])/100, np.log10(x[-1])) - ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 + #, np.log10(x[-1])) + ax_i.set_zlim(ax_i.get_zlim()[1]*0.02) + #ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file From 75af81d19458071a0fc5b78cb91468eafa2c66a6 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Mon, 12 May 2025 11:50:26 +0200 Subject: [PATCH 4/9] Add description, format code --- pdfplotter/pdf_set_nuclear.py | 257 ++++++++++++++++++++++++---------- 1 file changed, 181 insertions(+), 76 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 4c7c021..36d51c9 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -4,7 +4,7 @@ from math import log from typing import Sequence, Literal, Any -#from mpl_toolkits.mplot3d import Axes3D +# from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.ticker as mticker import matplotlib.cm as cm @@ -181,25 +181,68 @@ def plot_A_dep_3d( plot_uncertainty: bool = True, plot_ratio: bool = False, pdf_label: Literal["ylabel", "annotate"] = "annotate", - A_label: Literal["legend", "ticks"] = "ticks", + A_label: Literal["legend", "ticks", "both"] = "ticks", proj_type: Literal["ortho", "persp"] = "persp", view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_uncertainty_edges: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_title: dict[str, Any] = {}, - kwargs_notation: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_ylabel: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_xlabel: dict[str, Any] = {}, - kwargs_zlabel: dict[str, Any] = {}, - kwargs_legend: dict[str, Any] = {}, - kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_annotation: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_xlabel: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_ylabel: dict[str, Any] = {}, + kwargs_zlabel: dict[str, Any] = {}, + kwargs_legend: dict[str, Any] = {}, + kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {}, ): + """Plot the A-dependence of this PDF set in a 3d projection axis. + + Parameters + ---------- + ax : matplotlib.axes.Axes | numpy.ndarray[matplotlib.axes.Axes] + The axes to plot on. + A : float | list[float] + The A values to plot for. + observables : sympy.Basic | numpy.ndarray[sympy.Basic] | list[sympy.Basic] + The observables to plot. + Q : float, optional + The scale at which to plot the PDFs + Q2 : float, optional + The Q^2 scale at which to plot the PDFs. Either Q or Q2 has to be passed. + x_lines: + x values, ah which scattered lines are plotted to point out the A dep at fixed x + colors : list, optional + The colors to use for the different x values, by default [], tab color palette is used if == []. + plot_ratio : bool, optional + If True, plot the ratio of the PDFs to the Proton PDF, by default False + pdf_label : str, optional + The label for the PDF, by default "annotate". If "ylabel", it is used in ax.set_title(). If "annotate", the label is set as an annotation. + A_label: + If "ticks", the values for A are chosen as z-ticks. If "legend", a legend is plottet. if "both" both is realised + kwargs_theory : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the plot function for the central PDF, by default {}. + kwargs_uncertainty : dict[str, Any] | list[dict[str, Any] | None], optional + Additional keyword arguments for the PDF uncertainty band that should be passed to `plt.Axes.fill_between`, by default {}. If a list of keyword arguments is given, the i-th element is used for the i-th A value. + kwargs_uncertainty_edges : dict[str, Any] | list[dict[str, Any] | None], optional + Additional keyword arguments for the edges of the PDF uncertainty band that should be passed to `plt.Axes.plot`, by default {}. If a list of keyword arguments is given, the i-th element is used for the i-th A value. + kwargs_legend : dict[str, Any], optional + The keyword arguments to pass to the legend function, by default {}. + kwargs_xlabel : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the xlabel function, by default {}. + kwargs_ylabel : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the zlabel function, the A-axis, by default {}. + kwargs_zlabel : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the ylabel function, the f(x,Q)-axis,, by default {}. + kwargs_title : dict[str, Any], optional + The keyword arguments to pass to the title function, by default {}. + kwargs_annotation : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the annotate function, by default {}. + kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {} + The keyword arguments to pass to the plot function for plotting the lines at fixed x + """ my_sets = self.pdf_sets x = self.get(A=my_sets["A"][0]).x_values - if Q is None and Q2 is None: raise ValueError("Please pass either `Q` or `Q2`") @@ -226,10 +269,8 @@ def plot_A_dep_3d( A = [A] if 1 not in A and plot_ratio: - raise ValueError( - "Please pass A=1 if you want to plot the ratio to Proton." - ) - + raise ValueError("Please pass A=1 if you want to plot the ratio to Proton.") + if isinstance(observables, np.ndarray): observables = list(observables.flatten()) @@ -250,12 +291,16 @@ def plot_A_dep_3d( z_lower, z_upper = self.get(A=A_j).get_uncertainties( observable=obs_i, x=x, Q=Q, Q2=Q2 ) - else: + else: z_lower, z_upper = self.get(A=A_j).get_uncertainties( observable=obs_i, x=x, Q=Q, Q2=Q2 ) - z_lower = z_lower / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) - z_upper = z_upper / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) + z_lower = z_lower / self.get(A=1).get_central( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + z_upper = z_upper / self.get(A=1).get_central( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) kwargs_default = { "color": col_j, "label": f"A={A_j}", @@ -277,15 +322,21 @@ def plot_A_dep_3d( ax_i.plot( np.log10(x), np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) - / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ) + / self.get(A=1).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ), **kwargs, ) else: ax_i.plot( np.log10(x), np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ), **kwargs, ) else: @@ -293,17 +344,23 @@ def plot_A_dep_3d( ax_i.plot( np.log10(x), len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) - / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ) + / self.get(A=1).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ), **kwargs, ) else: ax_i.plot( np.log10(x), len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i + ), **kwargs, - ) + ) if plot_uncertainty: kwargs_uncertainty_default = { "color": col_j, @@ -331,20 +388,34 @@ def plot_A_dep_3d( ): vertices.append([xi, ai, zl]) for xi, ai, zl, zu in reversed( - list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) + list( + zip( + np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper + ) + ) ): vertices.append([xi, ai, zu]) - + else: for xi, ai, zl, zu in zip( - np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper + np.log10(x), + np.ones(len(x)) * np.log10(A_j), + z_lower, + z_upper, ): vertices.append([xi, ai, zl]) for xi, ai, zl, zu in reversed( - list(zip(np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper)) + list( + zip( + np.log10(x), + np.ones(len(x)) * np.log10(A_j), + z_lower, + z_upper, + ) + ) ): - vertices.append([xi, ai, zu]) + vertices.append([xi, ai, zu]) poly = Poly3DCollection([vertices], **kwargs) ax_i.add_collection3d(poly) @@ -367,14 +438,18 @@ def plot_A_dep_3d( ax_i.plot(np.log10(x), len(x) * [A_j], z_upper, **kwargs) ax_i.plot(np.log10(x), len(x) * [A_j], z_lower, **kwargs) else: - ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) - ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) + ax_i.plot( + np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs + ) + ax_i.plot( + np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs + ) - centrals={} + centrals = {} if x_lines is not None: if not isinstance(x_lines, list): x_lines = [x_lines] - for k,x_line in enumerate(x_lines): + for k, x_line in enumerate(x_lines): if x_line not in x: raise ValueError( f"Chosen x value {x_line} was not used for defining nuclear pdf set. \n Please choose x that was used in initialization" @@ -398,23 +473,53 @@ def plot_A_dep_3d( for a in A: if not plot_ratio: if x_line not in centrals.keys(): - centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + centrals[x_line] = [ + self.get(A=a).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + ] else: - centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + centrals[x_line].append( + self.get(A=a).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + ) else: if x_line not in centrals.keys(): - centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + centrals[x_line] = [ + self.get(A=a).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + / self.get(A=1).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + ] else: - centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + centrals[x_line].append( + self.get(A=a).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + / self.get(A=1).get_central( + x=x_line, Q=Q, Q2=Q2, observable=obs_i + ) + ) if logA: ax_i.plot( - np.ones(len(A))*np.log10(x_line), np.log10(A), centrals[x_line],**kwargs) + np.ones(len(A)) * np.log10(x_line), + np.log10(A), + centrals[x_line], + **kwargs, + ) else: ax_i.plot( - np.ones(len(A))*np.log10(x_line), A, centrals[x_line],**kwargs) + np.ones(len(A)) * np.log10(x_line), + A, + centrals[x_line], + **kwargs, + ) if pdf_label == "annotate": - kwargs_notation_default = { + kwargs_annotation_default = { "fontsize": 12, "xy": (0.97, 0.96), "xycoords": "axes fraction", @@ -427,15 +532,15 @@ def plot_A_dep_3d( boxstyle="round, pad=0.2", ), } - if not isinstance(kwargs_notation, list): + if not isinstance(kwargs_annotation, list): kwargs_n = update_kwargs( - kwargs_notation_default, - kwargs_notation, + kwargs_annotation_default, + kwargs_annotation, ) else: kwargs_n = update_kwargs( - kwargs_notation_default, - kwargs_notation, + kwargs_annotation_default, + kwargs_annotation, i=i, ) ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) @@ -459,7 +564,7 @@ def plot_A_dep_3d( ax_i.set_zlabel(**kwargs) else: - kwargs_notation_default = { + kwargs_annotation_default = { "fontsize": 12, "xy": (0.47, 0.96), "xycoords": "axes fraction", @@ -472,19 +577,18 @@ def plot_A_dep_3d( boxstyle="round, pad=0.2", ), } - kwargs_n = update_kwargs(kwargs_notation_default, kwargs_notation, i=i) + kwargs_n = update_kwargs(kwargs_annotation_default, kwargs_annotation, i=i) ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) - ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) - if A_label == "ticks": + ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) + if A_label == "ticks" or A_label == "both": if logA: - ax_i.set_yticks(np.log10(A),A) + ax_i.set_yticks(np.log10(A), A) else: - ax_i.set_yticks(A,A) + ax_i.set_yticks(A, A) kwargs_zlabel_default = { "fontsize": 14, "ylabel": f"$A$", - } kwargs = update_kwargs( kwargs_zlabel_default, @@ -492,10 +596,22 @@ def plot_A_dep_3d( ) ax_i.set_ylabel(**kwargs) - ax_i.set_xlim(np.log10(x[0])*0.98) - else: - ax_i.set_yticks([]) - if i==len(observables)-1: + ax_i.set_xlim(np.log10(x[0]) * 0.98) + if A_label == "legend" or A_label == "both": + if A_label == "legend": + ax_i.set_yticks([]) + kwargs_zlabel_default = { + "fontsize": 14, + "ylabel": f"$A$", + "labelpad": -10, + "linespacing": -4, + } + kwargs = update_kwargs( + kwargs_zlabel_default, + kwargs_zlabel, + ) + ax_i.set_ylabel(**kwargs) + if i == len(observables) - 1: kwargs_legend_default = { "fontsize": 12, "bbox_to_anchor": (0.95, 0.95), @@ -506,38 +622,27 @@ def plot_A_dep_3d( kwargs_legend, ) ax_i.legend() - kwargs_zlabel_default = { - "fontsize": 14, - "ylabel": f"$A$", - "labelpad":-10, - "linespacing": -4 - } - kwargs = update_kwargs( - kwargs_zlabel_default, - kwargs_zlabel, - ) - ax_i.set_ylabel(**kwargs) - ax_i.xaxis.pane.fill=False - ax_i.yaxis.pane.fill=False - ax_i.zaxis.pane.fill=False + + ax_i.xaxis.pane.fill = False + ax_i.yaxis.pane.fill = False + ax_i.zaxis.pane.fill = False ax_i.xaxis.pane.set_edgecolor("w") ax_i.yaxis.pane.set_edgecolor("w") ax_i.zaxis.pane.set_edgecolor("w") - ax_i.zaxis._axinfo["juggled"]=(1,2,0) + ax_i.zaxis._axinfo["juggled"] = (1, 2, 0) kwargs_xlabel_default = { "fontsize": 14, "xlabel": f"$x$", - } kwargs = update_kwargs( kwargs_xlabel_default, kwargs_xlabel, ) ax_i.set_xlabel(**kwargs) - #, np.log10(x[-1])) - ax_i.set_zlim(ax_i.get_zlim()[1]*0.02) - #ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 + # , np.log10(x[-1])) + ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) + # ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) - ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file + ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) From 8bffa435dd36217d2a226f30f48a8a658e5cf6f1 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Thu, 22 May 2025 14:17:01 +0200 Subject: [PATCH 5/9] Add option to plot ratios, add option to choose lin or log as scale for A --- pdfplotter/pdf_set_nuclear.py | 253 ++++++++++++++-------------------- 1 file changed, 105 insertions(+), 148 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 36d51c9..13c1280 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -177,9 +177,9 @@ def plot_A_dep_3d( Q2: float | None = None, x_lines: float | list[float] | None = None, colors: list = [], - logA: bool = True, + A_scale: Literal["log", "linlog", "lin"] = "log", plot_uncertainty: bool = True, - plot_ratio: bool = False, + ratio_to: PDFSet | None = None, pdf_label: Literal["ylabel", "annotate"] = "annotate", A_label: Literal["legend", "ticks", "both"] = "ticks", proj_type: Literal["ortho", "persp"] = "persp", @@ -216,7 +216,7 @@ def plot_A_dep_3d( If True, plot the ratio of the PDFs to the Proton PDF, by default False pdf_label : str, optional The label for the PDF, by default "annotate". If "ylabel", it is used in ax.set_title(). If "annotate", the label is set as an annotation. - A_label: + A_label: If "ticks", the values for A are chosen as z-ticks. If "legend", a legend is plottet. if "both" both is realised kwargs_theory : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the plot function for the central PDF, by default {}. @@ -229,7 +229,7 @@ def plot_A_dep_3d( kwargs_xlabel : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the xlabel function, by default {}. kwargs_ylabel : dict[str, Any] | list[dict[str, Any] | None], optional - The keyword arguments to pass to the zlabel function, the A-axis, by default {}. + The keyword arguments to pass to the zlabel function, the A-axis, by default {}. kwargs_zlabel : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the ylabel function, the f(x,Q)-axis,, by default {}. kwargs_title : dict[str, Any], optional @@ -268,9 +268,6 @@ def plot_A_dep_3d( if not isinstance(A, list): A = [A] - if 1 not in A and plot_ratio: - raise ValueError("Please pass A=1 if you want to plot the ratio to Proton.") - if isinstance(observables, np.ndarray): observables = list(observables.flatten()) @@ -287,20 +284,13 @@ def plot_A_dep_3d( for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): for j, (A_j, col_j) in enumerate(zip(A, colors)): - if not plot_ratio: - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) - else: - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) - z_lower = z_lower / self.get(A=1).get_central( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) - z_upper = z_upper / self.get(A=1).get_central( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) + z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to + )[0] + z_lower = [k if k>0 else 0 for k in self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to + )[1]] + kwargs_default = { "color": col_j, "label": f"A={A_j}", @@ -317,50 +307,17 @@ def plot_A_dep_3d( kwargs_theory, i=j, ) - if logA: - if plot_ratio: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ) - / self.get(A=1).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ), - **kwargs, - ) - else: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ), - **kwargs, - ) - else: - if plot_ratio: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ) - / self.get(A=1).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ), - **kwargs, - ) - else: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central( - x=x, Q=Q, Q2=Q2, observable=obs_i - ), - **kwargs, - ) + if A_scale == "log": + Aj_arr = np.log10(len(x) * [A_j]) + if A_scale == "lin": + Aj_arr = len(x) * [A_j] + ax_i.plot( + np.log10(x), + Aj_arr, + [k if k>0 else 0 for k in self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i,ratio_to=ratio_to)], + **kwargs, + ) + if plot_uncertainty: kwargs_uncertainty_default = { "color": col_j, @@ -381,41 +338,26 @@ def plot_A_dep_3d( vertices = [] z_lower = np.array(z_lower) z_upper = np.array(z_upper) - if not logA: - - for xi, ai, zl, zu in zip( - np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper - ): - vertices.append([xi, ai, zl]) - for xi, ai, zl, zu in reversed( - list( - zip( - np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper - ) - ) - ): - vertices.append([xi, ai, zu]) - else: - for xi, ai, zl, zu in zip( - np.log10(x), - np.ones(len(x)) * np.log10(A_j), - z_lower, - z_upper, - ): - vertices.append([xi, ai, zl]) - - for xi, ai, zl, zu in reversed( - list( - zip( - np.log10(x), - np.ones(len(x)) * np.log10(A_j), - z_lower, - z_upper, - ) + for xi, ai, zl, zu in zip( + np.log10(x), + Aj_arr, + z_lower, + z_upper, + ): + vertices.append([xi, ai, zl]) + + for xi, ai, zl, zu in reversed( + list( + zip( + np.log10(x), + Aj_arr, + z_lower, + z_upper, ) - ): - vertices.append([xi, ai, zu]) + ) + ): + vertices.append([xi, ai, zu]) poly = Poly3DCollection([vertices], **kwargs) ax_i.add_collection3d(poly) @@ -434,16 +376,9 @@ def plot_A_dep_3d( kwargs_uncertainty_edges, i=j, ) - if not logA: - ax_i.plot(np.log10(x), len(x) * [A_j], z_upper, **kwargs) - ax_i.plot(np.log10(x), len(x) * [A_j], z_lower, **kwargs) - else: - ax_i.plot( - np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs - ) - ax_i.plot( - np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs - ) + + ax_i.plot(np.log10(x), Aj_arr, z_upper, **kwargs) + ax_i.plot(np.log10(x), Aj_arr, z_lower, **kwargs) centrals = {} if x_lines is not None: @@ -471,53 +406,63 @@ def plot_A_dep_3d( i=k, ) for a in A: - if not plot_ratio: - if x_line not in centrals.keys(): - centrals[x_line] = [ - self.get(A=a).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) - ] - else: - centrals[x_line].append( - self.get(A=a).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) + + if x_line not in centrals.keys(): + centrals[x_line] = [ + self.get(A=a).get_central( + x=x_line, + Q=Q, + Q2=Q2, + observable=obs_i, + ratio_to=ratio_to, ) + ] else: - if x_line not in centrals.keys(): - centrals[x_line] = [ - self.get(A=a).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) - / self.get(A=1).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) - ] - else: - centrals[x_line].append( - self.get(A=a).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) - / self.get(A=1).get_central( - x=x_line, Q=Q, Q2=Q2, observable=obs_i - ) + centrals[x_line].append( + self.get(A=a).get_central( + x=x_line, + Q=Q, + Q2=Q2, + observable=obs_i, + ratio_to=ratio_to, ) - if logA: + ) + if A_scale == "log": ax_i.plot( np.ones(len(A)) * np.log10(x_line), np.log10(A), centrals[x_line], **kwargs, ) - else: + ax_i.plot( + [np.log10(x_line),np.log10(x_line)], + [np.log10(A[0]),np.log10(A[0])], + [0,self.get(A=A[0]).get_central( + x=x_line, + Q=Q, + Q2=Q2, + observable=obs_i, + ratio_to=ratio_to, + )], + **kwargs,) + elif A_scale == "lin": ax_i.plot( np.ones(len(A)) * np.log10(x_line), A, centrals[x_line], **kwargs, ) - + ax_i.plot( + [np.log10(x_line),np.log10(x_line)], + [A[0],A[0]], + [0,self.get(A=A[0]).get_central( + x=x_line, + Q=Q, + Q2=Q2, + observable=obs_i, + ratio_to=ratio_to, + )], + **kwargs,) if pdf_label == "annotate": kwargs_annotation_default = { "fontsize": 12, @@ -546,10 +491,16 @@ def plot_A_dep_3d( ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) if pdf_label == "ylabel": - kwargs_ylabel_default = { - "fontsize": 14, - "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", - } + if ratio_to: + kwargs_ylabel_default = { + "fontsize": 14, + "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2,R=True)}$", + } + else: + kwargs_ylabel_default = { + "fontsize": 14, + "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2,R=False)}$", + } if not isinstance(kwargs_ylabel, list): kwargs = update_kwargs( kwargs_ylabel_default, @@ -577,14 +528,16 @@ def plot_A_dep_3d( boxstyle="round, pad=0.2", ), } - kwargs_n = update_kwargs(kwargs_annotation_default, kwargs_annotation, i=i) + kwargs_n = update_kwargs( + kwargs_annotation_default, kwargs_annotation, i=i + ) ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) if A_label == "ticks" or A_label == "both": - if logA: + if A_scale == "log": ax_i.set_yticks(np.log10(A), A) - else: + elif A_scale == "lin": ax_i.set_yticks(A, A) kwargs_zlabel_default = { "fontsize": 14, @@ -642,7 +595,11 @@ def plot_A_dep_3d( ) ax_i.set_xlabel(**kwargs) # , np.log10(x[-1])) - ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) + if not ratio_to: + ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) + else: + ax_i.set_zlim(ax_i.get_zlim()[0], 2*np.median(self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i,ratio_to=ratio_to))) + ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) # ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) From fc29adf1cce2d86841adda754e34ce21bf716365 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Thu, 12 Jun 2025 14:11:43 +0200 Subject: [PATCH 6/9] Add asym uncertainties --- pdfplotter/pdf_set_nuclear.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 13c1280..42eee7e 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -179,6 +179,7 @@ def plot_A_dep_3d( colors: list = [], A_scale: Literal["log", "linlog", "lin"] = "log", plot_uncertainty: bool = True, + uncertainty_convention: Literal["sym", "asym"] = "asym", ratio_to: PDFSet | None = None, pdf_label: Literal["ylabel", "annotate"] = "annotate", A_label: Literal["legend", "ticks", "both"] = "ticks", @@ -285,10 +286,10 @@ def plot_A_dep_3d( for j, (A_j, col_j) in enumerate(zip(A, colors)): z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to + observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to, convention=uncertainty_convention, )[0] z_lower = [k if k>0 else 0 for k in self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to + observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to, convention=uncertainty_convention, )[1]] kwargs_default = { From 2e79f07a98e21355295fd1bc65227c879b32d0d4 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Thu, 12 Jun 2025 14:27:14 +0200 Subject: [PATCH 7/9] Change x_values to x, Q_values to Q, similarly to PDFSet class --- pdfplotter/pdf_set_nuclear.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 42eee7e..41a597b 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -242,7 +242,7 @@ def plot_A_dep_3d( """ my_sets = self.pdf_sets - x = self.get(A=my_sets["A"][0]).x_values + x = self.get(A=my_sets["A"][0]).x if Q is None and Q2 is None: raise ValueError("Please pass either `Q` or `Q2`") @@ -251,16 +251,16 @@ def plot_A_dep_3d( raise ValueError("Please pass either `Q` or `Q2`, not both") elif Q is not None: - if Q not in self.get(A=my_sets["A"][0]).Q_values and Q not in np.sqrt( - np.array(self.get(A=my_sets["A"][0]).Q2_values) + if Q not in self.get(A=my_sets["A"][0]).Q and Q not in np.sqrt( + np.array(self.get(A=my_sets["A"][0]).Q2) ): raise ValueError( f"Chosen Q value {Q} was not used for defining nuclear pdf set. \n Please choose Q that was used in initialization" ) else: if ( - Q2 not in self.get(A=my_sets["A"][0]).Q2_values - and Q2 not in np.array(self.get(A=my_sets["A"][0]).Q_values) ** 2 + Q2 not in self.get(A=my_sets["A"][0]).Q2 + and Q2 not in np.array(self.get(A=my_sets["A"][0]).Q) ** 2 ): raise ValueError( f"Chosen Q2 value {Q2} was not used for defining nuclear pdf set. \n Please choose Q2 that was used in initialization" From 538507fe5d78c1d56cdc6a3ea2319d27b77f176a Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Thu, 12 Jun 2025 14:34:27 +0200 Subject: [PATCH 8/9] Add formatting --- pdfplotter/pdf_set_nuclear.py | 70 ++++++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 18 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 41a597b..6c46cf1 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -286,11 +286,24 @@ def plot_A_dep_3d( for j, (A_j, col_j) in enumerate(zip(A, colors)): z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to, convention=uncertainty_convention, + observable=obs_i, + x=x, + Q=Q, + Q2=Q2, + ratio_to=ratio_to, + convention=uncertainty_convention, )[0] - z_lower = [k if k>0 else 0 for k in self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2, ratio_to=ratio_to, convention=uncertainty_convention, - )[1]] + z_lower = [ + k if k > 0 else 0 + for k in self.get(A=A_j).get_uncertainties( + observable=obs_i, + x=x, + Q=Q, + Q2=Q2, + ratio_to=ratio_to, + convention=uncertainty_convention, + )[1] + ] kwargs_default = { "color": col_j, @@ -315,7 +328,12 @@ def plot_A_dep_3d( ax_i.plot( np.log10(x), Aj_arr, - [k if k>0 else 0 for k in self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i,ratio_to=ratio_to)], + [ + k if k > 0 else 0 + for k in self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i, ratio_to=ratio_to + ) + ], **kwargs, ) @@ -435,17 +453,21 @@ def plot_A_dep_3d( centrals[x_line], **kwargs, ) - ax_i.plot( - [np.log10(x_line),np.log10(x_line)], - [np.log10(A[0]),np.log10(A[0])], - [0,self.get(A=A[0]).get_central( + ax_i.plot( + [np.log10(x_line), np.log10(x_line)], + [np.log10(A[0]), np.log10(A[0])], + [ + 0, + self.get(A=A[0]).get_central( x=x_line, Q=Q, Q2=Q2, observable=obs_i, ratio_to=ratio_to, - )], - **kwargs,) + ), + ], + **kwargs, + ) elif A_scale == "lin": ax_i.plot( np.ones(len(A)) * np.log10(x_line), @@ -453,17 +475,21 @@ def plot_A_dep_3d( centrals[x_line], **kwargs, ) - ax_i.plot( - [np.log10(x_line),np.log10(x_line)], - [A[0],A[0]], - [0,self.get(A=A[0]).get_central( + ax_i.plot( + [np.log10(x_line), np.log10(x_line)], + [A[0], A[0]], + [ + 0, + self.get(A=A[0]).get_central( x=x_line, Q=Q, Q2=Q2, observable=obs_i, ratio_to=ratio_to, - )], - **kwargs,) + ), + ], + **kwargs, + ) if pdf_label == "annotate": kwargs_annotation_default = { "fontsize": 12, @@ -599,7 +625,15 @@ def plot_A_dep_3d( if not ratio_to: ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) else: - ax_i.set_zlim(ax_i.get_zlim()[0], 2*np.median(self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i,ratio_to=ratio_to))) + ax_i.set_zlim( + ax_i.get_zlim()[0], + 2 + * np.median( + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i, ratio_to=ratio_to + ) + ), + ) ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) # ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) From 947872438b2f7b72636b561dd4d75be02172c8b3 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Thu, 12 Jun 2025 16:06:50 +0200 Subject: [PATCH 9/9] Fix zlim offset --- pdfplotter/pdf_set_nuclear.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 6c46cf1..fb8f94b 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -283,7 +283,7 @@ def plot_A_dep_3d( colors = [cmap(i) for i in range(len(A))] for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): - + zlim_bottom = np.inf for j, (A_j, col_j) in enumerate(zip(A, colors)): z_upper = self.get(A=A_j).get_uncertainties( observable=obs_i, @@ -336,6 +336,14 @@ def plot_A_dep_3d( ], **kwargs, ) + zlim_bottom = min( + list( + self.get(A=A_j).get_central( + x=x, Q=Q, Q2=Q2, observable=obs_i, ratio_to=ratio_to + ) + ) + + [zlim_bottom] + ) if plot_uncertainty: kwargs_uncertainty_default = { @@ -623,7 +631,7 @@ def plot_A_dep_3d( ax_i.set_xlabel(**kwargs) # , np.log10(x[-1])) if not ratio_to: - ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) + ax_i.set_zlim(zlim_bottom) else: ax_i.set_zlim( ax_i.get_zlim()[0], @@ -634,7 +642,7 @@ def plot_A_dep_3d( ) ), ) - ax_i.set_zlim(ax_i.get_zlim()[1] * 0.02) + ax_i.set_zlim(zlim_bottom) # ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init)