Skip to content

Python: put scatplot bias and std dev on figure and remove pval #15705

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rmcdermo merged 1 commit into from
Dec 11, 2025
+44 −49
Merged

Python: put scatplot bias and std dev on figure and remove pval #15705

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Utilities/Python/FDS_validation_scatterplot_inputs.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,4 +39,4 @@ Plume Height,Empirical Plume Height (m),Predicted Plume Height (m),10,1000,0.03
Rate of Spread,Measured Rate of Spread (m/s),Predicted Rate of Spread (m/s),0.001,3,0.03 0.96,EastOutside,1.6,10,yes,loglog,FDS_Validation_Guide/SCRIPT_FIGURES/Scatterplots/FDS_Rate_of_Spread
Condensation,Measured Heat Flux (kW/m²),Predicted Heat Flux (kW/m²),0,12,0.03 0.96,EastOutside,1.5,12,yes,linear,FDS_Validation_Guide/SCRIPT_FIGURES/Scatterplots/FDS_Condensation
Aerosol Deposition,Measured Mass (mg),Predicted Mass (mg),0,0.15,0.03 0.96,EastOutside,1.5,12,yes,linear,FDS_Validation_Guide/SCRIPT_FIGURES/Scatterplots/FDS_Aerosol_Deposition
Evaporation Rate,Measured Rate (dD²/dt),Predicted Rate (dD²/dt),1.00E-09,1.00E-07,0.03 0.95,EastOutside,1.5,-1,yes,loglog,FDS_Validation_Guide/SCRIPT_FIGURES/Scatterplots/FDS_Evaporation_Rate
Evaporation Rate,Measured Rate (dD²/dt),Predicted Rate (dD²/dt),1.00E-09,1.00E-07,0.03 0.95,EastOutside,1.3,-1,yes,loglog,FDS_Validation_Guide/SCRIPT_FIGURES/Scatterplots/FDS_Evaporation_Rate
91 changes: 43 additions & 48 deletions Utilities/Python/fdsplotlib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1126,6 +1126,21 @@ def plot_to_fig(x_data,y_data,**kwargs):
if isinstance(data_label, str) and data_label.lower() == 'blank':
data_label = None

# ------------------------------------------------------------------
# PRE-PATCH FOR LOG AXES — PREVENTS log(0) WARNINGS
# ------------------------------------------------------------------
xmin = kwargs.get('x_min')
xmax = kwargs.get('x_max')
ymin = kwargs.get('y_min')
ymax = kwargs.get('y_max')

if plot_type in ('loglog', 'semilogx', 'semilogy'):
eps = 1e-12
if xmin is not None and xmin <= 0: xmin = eps
if ymin is not None and ymin <= 0: ymin = eps
if xmin is not None and xmax is not None: ax.set_xlim(xmin, xmax)
if ymin is not None and ymax is not None: ax.set_ylim(ymin, ymax)

# generate the main x,y plot
if plot_type=='linear':
ax.plot(x_data,y_data,
Expand Down Expand Up @@ -1225,13 +1240,11 @@ def plot_to_fig(x_data,y_data,**kwargs):
capthick=linewidth,
)


ticklabel_fontsize=kwargs.get('ticklabel_fontsize',default_ticklabel_fontsize)
plt.setp( ax.xaxis.get_majorticklabels(), rotation=0, fontsize=ticklabel_fontsize )
plt.setp( ax.yaxis.get_majorticklabels(), rotation=0, fontsize=ticklabel_fontsize )

axeslabel_fontsize=kwargs.get('axeslabel_fontsize',default_axeslabel_fontsize)

if not using_existing_figure:
plt.xlabel(kwargs.get('x_label'), fontsize=axeslabel_fontsize)
plt.ylabel(kwargs.get('y_label'), fontsize=axeslabel_fontsize)
Expand Down Expand Up @@ -2141,7 +2154,7 @@ def _label_at(k, labels):

# --- Call histogram BEFORE filtering to match MATLAB (includes zeros/Infs) ---
try:
hist_file, pval = statistics_histogram(
hist_file = statistics_histogram(
Measured_Values, Predicted_Values,
Plot_Filename, Manuals_Dir, Scatter_Plot_Title
)
Expand All @@ -2164,16 +2177,23 @@ def _label_at(k, labels):
print(f"[scatplot] Skipping {Scatter_Plot_Title} (no valid data)")
continue

# --- Compute statistics ---
# --- Compute statistics (MATLAB logic restored) ---
weight = np.ones_like(Measured_Values)
log_E_bar = np.sum(np.log(Measured_Values) * weight) / np.sum(weight)
log_M_bar = np.sum(np.log(Predicted_Values) * weight) / np.sum(weight)
u2 = np.sum((((np.log(Predicted_Values) - np.log(Measured_Values))
- (log_M_bar - log_E_bar)) ** 2) * weight) / (np.sum(weight) - 1)
u = np.sqrt(u2)
Sigma_E = min(u / np.sqrt(2), Sigma_E_input / 100.0)
Sigma_M = np.sqrt(max(0.0, u ** 2 - Sigma_E ** 2))
delta = np.exp(log_M_bar - log_E_bar + 0.5 * Sigma_M ** 2 - 0.5 * Sigma_E ** 2)

# Restore MATLAB logic:
# If no Sigma_E is supplied, experimental sigma = u/sqrt(2)
if Sigma_E_input > 0:
Sigma_E = Sigma_E_input / 100.0
else:
Sigma_E = u / np.sqrt(2)

Sigma_M = np.sqrt(max(0.0, u**2 - Sigma_E**2))
delta = np.exp(log_M_bar - log_E_bar + 0.5 * Sigma_M**2 - 0.5 * Sigma_E**2)

# --- Scatter Plot ---
fig = fdsplotlib.plot_to_fig(x_data=[-1], y_data=[-1],
Expand All @@ -2188,12 +2208,21 @@ def _label_at(k, labels):
legend_expand=row["Paper_Width_Factor"])

fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=[Plot_Min, Plot_Max], line_style="k-", figure_handle=fig)
if Sigma_E > 0:
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * (1 + 2 * Sigma_E), line_style="k--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) / (1 + 2 * Sigma_E), line_style="k--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta, line_style="r-", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta * (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta / (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * (1 + 2 * Sigma_E), line_style="k--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) / (1 + 2 * Sigma_E), line_style="k--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta, line_style="r-", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta * (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig)
fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta / (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig)

# --- Add statistics annotation ---
ax = fig.gca()
xnorm = 0.03 # 3% from left
ynorm = 0.97 # 3% from top
dy = 0.05 # vertical spacing in normalized units
ax.text(xnorm, ynorm,f"{Scatter_Plot_Title}",fontsize=plot_style["Scat_Title_Font_Size"],ha="left", va="top",transform=ax.transAxes)
ax.text(xnorm, ynorm - dy,f"Exp. Rel. Std. Dev.: {Sigma_E:.2f}",fontsize=plot_style["Scat_Label_Font_Size"],ha="left", va="top",transform=ax.transAxes)
ax.text(xnorm, ynorm - 2*dy,f"Model Rel. Std. Dev.: {Sigma_M:.2f}",fontsize=plot_style["Scat_Label_Font_Size"],ha="left", va="top",transform=ax.transAxes)
ax.text(xnorm, ynorm - 3*dy,f"Model Bias Factor: {delta:.2f}",fontsize=plot_style["Scat_Label_Font_Size"],ha="left", va="top",transform=ax.transAxes)

# --- Plot each dataset with its own marker and color ---
seen_labels = set()
Expand Down Expand Up @@ -2528,30 +2557,6 @@ def histogram_output(Histogram_Tex_Output, Output_Histograms):
print(f"[histogram_output] Wrote LaTeX histogram file: {Histogram_Tex_Output}")


def spiegel_test(x):
"""Exact translation of MATLAB spiegel_test.m with Inf/NaN handling identical to MATLAB."""
import numpy as np
from math import sqrt, erf

x = np.asarray(x, dtype=float)

# MATLAB arithmetic allows Inf and NaN to propagate but doesn't remove them
xm = np.nanmean(x)
xs = np.nanstd(x, ddof=0)
xz = (x - xm) / xs
xz2 = xz ** 2

# MATLAB behavior: Inf * 0 -> NaN, which is ignored by sum()
with np.errstate(invalid='ignore', divide='ignore'):
term = xz2 * np.log(xz2)

N = np.nansum(term)
n = np.sum(np.isfinite(xz2)) # count finite entries only
ts = (N - 0.73 * n) / (0.8969 * sqrt(n))
pval = 1 - abs(erf(ts / sqrt(2))) # 2-sided
return pval


def statistics_histogram(Measured_Values, Predicted_Values,
Plot_Filename, Manuals_Dir, Scatter_Plot_Title,
Figure_Visibility='off', Paper_Width=5.0,
Expand All @@ -2566,13 +2571,6 @@ def statistics_histogram(Measured_Values, Predicted_Values,
with np.errstate(divide='ignore', invalid='ignore'):
ln_M_E = np.log(Predicted_Values) - np.log(Measured_Values)

# Normality test before filtering — this is the key
if len(ln_M_E) >= 4:
pval = spiegel_test(ln_M_E)
else:
print(f"[statistics_histogram] Not enough data for {Scatter_Plot_Title}")
return None, None

# MATLAB's hist() ignores NaN/Inf implicitly
valid = np.isfinite(ln_M_E)
ln_M_E = ln_M_E[valid]
Expand Down Expand Up @@ -2600,8 +2598,6 @@ def statistics_histogram(Measured_Values, Predicted_Values,
ax.set_xticks(xcenters)
ax.set_xticklabels([str(i) for i in range(1, len(xcenters) + 1)])
ax.text(0.03, 0.90, Scatter_Plot_Title, transform=ax.transAxes)
ax.text(0.03, 0.82, "Normality Test", transform=ax.transAxes)
ax.text(0.03, 0.74, f"p-value = {pval:.2f}", transform=ax.transAxes)

outpath = os.path.join(Manuals_Dir, f"{Plot_Filename}_Histogram.pdf")
os.makedirs(os.path.dirname(outpath), exist_ok=True)
Expand All @@ -2610,8 +2606,7 @@ def statistics_histogram(Measured_Values, Predicted_Values,
plt.close(fig)
plt.figure().clear()

print(f"[statistics_histogram] {Scatter_Plot_Title}: p = {pval:.2f}")
return f"{os.path.basename(Plot_Filename)}_Histogram", pval
return f"{os.path.basename(Plot_Filename)}_Histogram"


def set_ticks_like_matlab(fig):
Expand Down