Preserving isosurfaces with safeguards¶
In this example, we compute the isosurface from a 3D dataset of wind u and pressure anomaly p during the 2003 hurricane Isabel. We compare how three different lossy compressors (ZFP, SZ3, and SPERR) affect the isosurface. Finally, we apply safeguards to guarantee that the isosurface of interest is preserved. We also showcase how the safeguards can be used to preserve an arbitrary number of isosurfaces.
QPET supports preserving isosurfaces. However, at the time of writing, this is not yet implemented in QPET-SPERR. Therefore, we do not compare with QPET in this example.
Copied!
import ssl
ssl._create_default_https_context = ssl._create_stdlib_context
import ssl
ssl._create_default_https_context = ssl._create_stdlib_context
Copied!
from pathlib import Path
import humanize
import numpy as np
import pandas as pd
from matplotlib import gridspec
from matplotlib import pyplot as plt
from pathlib import Path
import humanize
import numpy as np
import pandas as pd
from matplotlib import gridspec
from matplotlib import pyplot as plt
Copied!
# Retrieve the data,
# which is stored as a big endian float32 block of shape z*y*x = 100*500*500
pf48 = (
np.fromfile(
Path() / "data" / "isabel" / "Pf48.bin",
dtype=">f4",
count=500 * 500 * 100,
sep="",
)
.reshape(100, 500, 500)
.astype(np.float32)
)
pf48[pf48 == 1.0e35] = np.nan
uf48 = (
np.fromfile(
Path() / "data" / "isabel" / "Uf48.bin",
dtype=">f4",
count=500 * 500 * 100,
sep="",
)
.reshape(100, 500, 500)
.astype(np.float32)
)
uf48[uf48 == 1.0e35] = np.nan
# Retrieve the data,
# which is stored as a big endian float32 block of shape z*y*x = 100*500*500
pf48 = (
np.fromfile(
Path() / "data" / "isabel" / "Pf48.bin",
dtype=">f4",
count=500 * 500 * 100,
sep="",
)
.reshape(100, 500, 500)
.astype(np.float32)
)
pf48[pf48 == 1.0e35] = np.nan
uf48 = (
np.fromfile(
Path() / "data" / "isabel" / "Uf48.bin",
dtype=">f4",
count=500 * 500 * 100,
sep="",
)
.reshape(100, 500, 500)
.astype(np.float32)
)
uf48[uf48 == 1.0e35] = np.nan
Copied!
def compute_corrections_percentage(my_f: np.ndarray, f: np.ndarray) -> float:
return np.mean(~((my_f == f) | (np.isnan(my_f) & np.isnan(f))))
def compute_corrections_percentage(my_f: np.ndarray, f: np.ndarray) -> float:
return np.mean(~((my_f == f) | (np.isnan(my_f) & np.isnan(f))))
Copied!
from numpy.lib.stride_tricks import sliding_window_view
def compute_failures(x, xnew, level):
cells_x = sliding_window_view(x < level, (2, 2, 2))[::2, ::2, ::2].reshape(
tuple(np.array(x.shape) // 2) + (2 * 2 * 2,)
)
cells_xnew = sliding_window_view(xnew < level, (2, 2, 2))[::2, ::2, ::2].reshape(
tuple(np.array(x.shape) // 2) + (2 * 2 * 2,)
)
all_x = np.logical_and.reduce(cells_x, axis=-1)
any_x = np.logical_or.reduce(cells_x, axis=-1)
any_xnew = np.logical_or.reduce(cells_xnew, axis=-1)
iso_x = any_x & ~all_x
neq = np.logical_or.reduce(cells_x != cells_xnew, axis=-1)
fn = any_x & ~any_xnew
fp = ~any_x & any_xnew
fs = any_x & any_xnew & neq
return iso_x, fn, fp, fs
from numpy.lib.stride_tricks import sliding_window_view
def compute_failures(x, xnew, level):
cells_x = sliding_window_view(x < level, (2, 2, 2))[::2, ::2, ::2].reshape(
tuple(np.array(x.shape) // 2) + (2 * 2 * 2,)
)
cells_xnew = sliding_window_view(xnew < level, (2, 2, 2))[::2, ::2, ::2].reshape(
tuple(np.array(x.shape) // 2) + (2 * 2 * 2,)
)
all_x = np.logical_and.reduce(cells_x, axis=-1)
any_x = np.logical_or.reduce(cells_x, axis=-1)
any_xnew = np.logical_or.reduce(cells_xnew, axis=-1)
iso_x = any_x & ~all_x
neq = np.logical_or.reduce(cells_x != cells_xnew, axis=-1)
fn = any_x & ~any_xnew
fp = ~any_x & any_xnew
fs = any_x & any_xnew & neq
return iso_x, fn, fp, fs
Copied!
from skimage.measure import marching_cubes
def plot_isosurface(
my_f,
f,
cr,
ax,
title,
var,
label,
level,
my_span,
span,
my_n=22,
n=15,
error=False,
corr=None,
):
ax.set_xlim(0, 500)
ax.set_ylim(0, 500)
ax.set_zlim(0, 100)
ax.set_xticks(np.arange(0, 501, 500 / 3))
ax.set_xticks(np.arange(0, 501, 500 / 6), minor=True)
ax.set_xticklabels(
[rf"${x}\degree$W" for x in np.linspace(83, 62, 4)], color="#777"
)
ax.set_yticks(np.arange(0, 501, 500 / 3)[:-1])
ax.set_yticks(np.arange(0, 501, 500 / 6), minor=True)
ax.set_yticklabels(
[rf"${y}\degree$N" for y in np.linspace(23.7, 41.7, 4)[:-1]],
color="#777",
va="bottom",
ha="left",
)
ax.set_zticks(np.arange(0, 101, 100 / 3))
ax.set_zticks(np.arange(0, 101, 100 / 6), minor=True)
ax.set_zticklabels(
[f"${z}$km" for z in np.round(np.linspace(0.035, 19.835, 4), 3)],
color="#777",
ha="left",
)
levels = np.linspace(-span, span + 1, n)
err_levels = np.linspace(-my_span, my_span, my_n)
with np.errstate(invalid="ignore"):
cb_f = (
np.where(
np.isnan(my_f) != np.isnan(f),
np.nan_to_num(my_f) - np.nan_to_num(f),
my_f - f,
)
if error
else f
)
ct = my_span if error else span
cb_b = my_n if error else n
extend_left = np.nanmin(cb_f) < -ct
extend_right = np.nanmax(cb_f) > ct
extend = {
(False, False): "neither",
(True, False): "min",
(False, True): "max",
(True, True): "both",
}[(extend_left, extend_right)]
x = np.broadcast_to(np.arange(500).reshape(1, 500), (500, 500))
y = np.broadcast_to(np.arange(500)[::-1].reshape(500, 1), (500, 500))
z = np.broadcast_to(0, (500, 500))
ax.fill_between(
[0, 490],
[10, 10],
[0, 0],
[0, 490],
[500, 500],
[0, 0],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
cm = ax.contourf(
x,
y,
cb_f[0].T,
zdir="z",
offset=0,
cmap="coolwarm",
levels=err_levels,
zorder=-503,
extend=extend,
)
cm.set_visible(False)
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[0].T)),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
cm = ax.contourf(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
x = np.broadcast_to(0, (100, 500))
y = np.broadcast_to(np.arange(500)[::-1].reshape(1, 500), (100, 500))
z = np.broadcast_to(np.arange(100).reshape(100, 1), (100, 500))
ax.fill_between(
[0, 0],
[10, 500],
[0, 0],
[0, 0],
[10, 500],
[98, 98],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[:, 0, :])),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
ax.contourf(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
x = np.broadcast_to(np.arange(500).reshape(1, 500), (100, 500))
y = np.broadcast_to(499, (100, 500))
z = np.broadcast_to(np.arange(100).reshape(100, 1), (100, 500))
ax.fill_between(
[0, 490],
[500, 500],
[0, 0],
[0, 490],
[500, 500],
[98, 98],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[:, :, 0])),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
ax.contourf(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
verts, faces, _, _ = marching_cubes(
my_f.transpose(1, 2, 0)[:, ::-1, :], level=level, step_size=2
)
ax.plot_trisurf(
verts[:, 0],
verts[:, 1],
faces,
verts[:, 2],
color=plt.cm.PuOr_r(np.linspace(0, 1, 14))[5],
edgecolor="none",
# alpha=0.5,
zorder=-499,
rasterized=True,
)
ax.set_rasterization_zorder(-501)
ax.set_title(title, color="#333333")
if error:
err_v = np.mean(
((my_f < level) != (f < level))
| ((my_f > level) != (f > level))
| (np.isnan(my_f) != np.isnan(f))
)
err_v = (
0
if err_v == 0
else np.format_float_positional(100 * err_v, precision=1, min_digits=1)
+ "%"
)
if err_v == "0.0%":
err_v = "<0.05%"
t = ax.text2D(
0.95,
-0.075,
f"V={err_v}",
ha="right",
va="bottom",
transform=ax.transAxes,
color="#333333",
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax.text2D(
0.95,
0.925,
rf"$\times$ {np.round(cr, 2)}"
if error
else humanize.naturalsize(f.nbytes, binary=True),
ha="right",
va="top",
transform=ax.transAxes,
color="#333333",
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax.grid(False)
for xy in np.arange(0, 501, 500 / 6):
ax.plot([xy, xy], [0, 500], [0, 0], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([xy, xy], [500, 500], [0, 100], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 500], [xy, xy], [0, 0], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 0], [xy, xy], [0, 100], color="#AAAAAA", lw=0.8, zorder=-500)
for z in np.arange(0, 101, 100 / 6):
ax.plot([0, 500], [500, 500], [z, z], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 0], [0, 500], [z, z], color="#AAAAAA", lw=0.8, zorder=-500)
cb = ax.figure.colorbar(
cm, ax=ax, location="bottom", pad=0.125, ticks=np.linspace(-ct, ct, 5)
)
cb.outline.set_color("#AAAAAA")
cb.ax.tick_params(pad=0, width=0, labelcolor="#333333")
cb.ax.set_xlabel(label, color="#333333")
counts, bins = np.histogram(cb_f, range=(-ct, ct), bins=cb_b - 1)
midpoints = bins[:-1] + np.diff(bins) / 2
extend_width = (bins[-1] - bins[-2]) / (bins[-1] - bins[0])
cax = cb.ax.inset_axes(
[
0.0 - extend_width * extend_left,
1.25,
1.0 + extend_width * (0 + extend_left + extend_right),
1.0,
]
)
cax.bar(
midpoints,
height=counts,
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints)),
)
if extend_left:
cax.bar(
bins[0] - (bins[1] - bins[0]) / 2,
height=np.sum(cb_f < -ct),
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints[0])),
)
if extend_right:
cax.bar(
bins[-1] + (bins[-1] - bins[-2]) / 2,
height=np.sum(cb_f > ct),
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints[-1])),
)
q1, q2, q3 = np.nanquantile(cb_f, [0.25, 0.5, 0.75])
cax.axvline(np.nanmean(cb_f), ls=":", ymin=0.1, ymax=0.9, c="w", lw=2)
cax.axvline(q1, ymin=0.25, ymax=0.75, c="w", lw=2)
cax.axvline(q2, ymin=0.1, ymax=0.9, c="w", lw=2)
cax.axvline(q3, ymin=0.25, ymax=0.75, c="w", lw=2)
cax.axvline(np.nanmean(cb_f), ymin=0.1, ymax=0.9, ls=":", c="k", lw=1)
cax.axvline(q1, ymin=0.25, ymax=0.75, c="k", lw=1)
cax.axvline(q2, ymin=0.1, ymax=0.9, c="k", lw=1)
cax.axvline(q3, ymin=0.25, ymax=0.75, c="k", lw=1)
cax.set_xlim(
-ct - (bins[-1] - bins[-2]) * extend_left,
ct + (bins[-1] - bins[-2]) * extend_right,
)
cax.set_xticks([])
cax.set_yticks([])
cax.spines[:].set_visible(False)
from skimage.measure import marching_cubes
def plot_isosurface(
my_f,
f,
cr,
ax,
title,
var,
label,
level,
my_span,
span,
my_n=22,
n=15,
error=False,
corr=None,
):
ax.set_xlim(0, 500)
ax.set_ylim(0, 500)
ax.set_zlim(0, 100)
ax.set_xticks(np.arange(0, 501, 500 / 3))
ax.set_xticks(np.arange(0, 501, 500 / 6), minor=True)
ax.set_xticklabels(
[rf"${x}\degree$W" for x in np.linspace(83, 62, 4)], color="#777"
)
ax.set_yticks(np.arange(0, 501, 500 / 3)[:-1])
ax.set_yticks(np.arange(0, 501, 500 / 6), minor=True)
ax.set_yticklabels(
[rf"${y}\degree$N" for y in np.linspace(23.7, 41.7, 4)[:-1]],
color="#777",
va="bottom",
ha="left",
)
ax.set_zticks(np.arange(0, 101, 100 / 3))
ax.set_zticks(np.arange(0, 101, 100 / 6), minor=True)
ax.set_zticklabels(
[f"${z}$km" for z in np.round(np.linspace(0.035, 19.835, 4), 3)],
color="#777",
ha="left",
)
levels = np.linspace(-span, span + 1, n)
err_levels = np.linspace(-my_span, my_span, my_n)
with np.errstate(invalid="ignore"):
cb_f = (
np.where(
np.isnan(my_f) != np.isnan(f),
np.nan_to_num(my_f) - np.nan_to_num(f),
my_f - f,
)
if error
else f
)
ct = my_span if error else span
cb_b = my_n if error else n
extend_left = np.nanmin(cb_f) < -ct
extend_right = np.nanmax(cb_f) > ct
extend = {
(False, False): "neither",
(True, False): "min",
(False, True): "max",
(True, True): "both",
}[(extend_left, extend_right)]
x = np.broadcast_to(np.arange(500).reshape(1, 500), (500, 500))
y = np.broadcast_to(np.arange(500)[::-1].reshape(500, 1), (500, 500))
z = np.broadcast_to(0, (500, 500))
ax.fill_between(
[0, 490],
[10, 10],
[0, 0],
[0, 490],
[500, 500],
[0, 0],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
cm = ax.contourf(
x,
y,
cb_f[0].T,
zdir="z",
offset=0,
cmap="coolwarm",
levels=err_levels,
zorder=-503,
extend=extend,
)
cm.set_visible(False)
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[0].T)),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
cm = ax.contourf(
x,
y,
my_f[0].T,
zdir="z",
offset=0,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
x = np.broadcast_to(0, (100, 500))
y = np.broadcast_to(np.arange(500)[::-1].reshape(1, 500), (100, 500))
z = np.broadcast_to(np.arange(100).reshape(100, 1), (100, 500))
ax.fill_between(
[0, 0],
[10, 500],
[0, 0],
[0, 0],
[10, 500],
[98, 98],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[:, 0, :])),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
ax.contourf(
my_f[:, 0, :],
y,
z,
zdir="x",
offset=0,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
x = np.broadcast_to(np.arange(500).reshape(1, 500), (100, 500))
y = np.broadcast_to(499, (100, 500))
z = np.broadcast_to(np.arange(100).reshape(100, 1), (100, 500))
ax.fill_between(
[0, 490],
[500, 500],
[0, 0],
[0, 490],
[500, 500],
[98, 98],
mode="polygon",
hatch="XX",
edgecolor="magenta",
facecolor="lavenderblush",
zorder=-504,
)
if error:
ax.plot_surface(
x,
y,
z,
rstride=1,
cstride=1,
facecolors=cm.cmap(cm.norm(cb_f[:, :, 0])),
shade=False,
zorder=-503,
)
with np.errstate(invalid="ignore"):
ax.contour(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
colors="black",
linewidths=2,
levels=levels,
zorder=-502,
extend=extend,
)
ax.contour(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
cmap="PuOr_r",
linewidths=1,
levels=levels,
zorder=-501,
extend=extend,
)
else:
ax.contourf(
x,
my_f[:, :, 0],
z,
zdir="y",
offset=499,
cmap="PuOr_r",
levels=levels,
zorder=-503,
extend=extend,
)
verts, faces, _, _ = marching_cubes(
my_f.transpose(1, 2, 0)[:, ::-1, :], level=level, step_size=2
)
ax.plot_trisurf(
verts[:, 0],
verts[:, 1],
faces,
verts[:, 2],
color=plt.cm.PuOr_r(np.linspace(0, 1, 14))[5],
edgecolor="none",
# alpha=0.5,
zorder=-499,
rasterized=True,
)
ax.set_rasterization_zorder(-501)
ax.set_title(title, color="#333333")
if error:
err_v = np.mean(
((my_f < level) != (f < level))
| ((my_f > level) != (f > level))
| (np.isnan(my_f) != np.isnan(f))
)
err_v = (
0
if err_v == 0
else np.format_float_positional(100 * err_v, precision=1, min_digits=1)
+ "%"
)
if err_v == "0.0%":
err_v = "<0.05%"
t = ax.text2D(
0.95,
-0.075,
f"V={err_v}",
ha="right",
va="bottom",
transform=ax.transAxes,
color="#333333",
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax.text2D(
0.95,
0.925,
rf"$\times$ {np.round(cr, 2)}"
if error
else humanize.naturalsize(f.nbytes, binary=True),
ha="right",
va="top",
transform=ax.transAxes,
color="#333333",
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax.grid(False)
for xy in np.arange(0, 501, 500 / 6):
ax.plot([xy, xy], [0, 500], [0, 0], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([xy, xy], [500, 500], [0, 100], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 500], [xy, xy], [0, 0], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 0], [xy, xy], [0, 100], color="#AAAAAA", lw=0.8, zorder=-500)
for z in np.arange(0, 101, 100 / 6):
ax.plot([0, 500], [500, 500], [z, z], color="#AAAAAA", lw=0.8, zorder=-500)
ax.plot([0, 0], [0, 500], [z, z], color="#AAAAAA", lw=0.8, zorder=-500)
cb = ax.figure.colorbar(
cm, ax=ax, location="bottom", pad=0.125, ticks=np.linspace(-ct, ct, 5)
)
cb.outline.set_color("#AAAAAA")
cb.ax.tick_params(pad=0, width=0, labelcolor="#333333")
cb.ax.set_xlabel(label, color="#333333")
counts, bins = np.histogram(cb_f, range=(-ct, ct), bins=cb_b - 1)
midpoints = bins[:-1] + np.diff(bins) / 2
extend_width = (bins[-1] - bins[-2]) / (bins[-1] - bins[0])
cax = cb.ax.inset_axes(
[
0.0 - extend_width * extend_left,
1.25,
1.0 + extend_width * (0 + extend_left + extend_right),
1.0,
]
)
cax.bar(
midpoints,
height=counts,
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints)),
)
if extend_left:
cax.bar(
bins[0] - (bins[1] - bins[0]) / 2,
height=np.sum(cb_f < -ct),
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints[0])),
)
if extend_right:
cax.bar(
bins[-1] + (bins[-1] - bins[-2]) / 2,
height=np.sum(cb_f > ct),
width=(bins[-1] - bins[0]) / len(counts),
color=cb.cmap(cb.norm(midpoints[-1])),
)
q1, q2, q3 = np.nanquantile(cb_f, [0.25, 0.5, 0.75])
cax.axvline(np.nanmean(cb_f), ls=":", ymin=0.1, ymax=0.9, c="w", lw=2)
cax.axvline(q1, ymin=0.25, ymax=0.75, c="w", lw=2)
cax.axvline(q2, ymin=0.1, ymax=0.9, c="w", lw=2)
cax.axvline(q3, ymin=0.25, ymax=0.75, c="w", lw=2)
cax.axvline(np.nanmean(cb_f), ymin=0.1, ymax=0.9, ls=":", c="k", lw=1)
cax.axvline(q1, ymin=0.25, ymax=0.75, c="k", lw=1)
cax.axvline(q2, ymin=0.1, ymax=0.9, c="k", lw=1)
cax.axvline(q3, ymin=0.25, ymax=0.75, c="k", lw=1)
cax.set_xlim(
-ct - (bins[-1] - bins[-2]) * extend_left,
ct + (bins[-1] - bins[-2]) * extend_right,
)
cax.set_xticks([])
cax.set_yticks([])
cax.spines[:].set_visible(False)
Copied!
def table_isosurface(
my_f,
f,
cr,
title,
var,
level,
corr,
) -> pd.DataFrame:
with np.errstate(over="ignore", invalid="ignore"):
err_diff = my_f - f
err_diff[np.isnan(my_f) & np.isnan(f)] = 0
err_inf = np.amax(np.abs(err_diff))
err_inf_fin = np.nanmax(np.abs(err_diff))
err_2 = np.sqrt(np.mean(np.square(err_diff)))
err_2_fin = np.sqrt(np.nanmean(np.square(err_diff)))
isof, fnv, fpv, fsv = compute_failures(f, my_f, level)
ison = int(np.sum(isof))
fn = int(np.sum(fnv)) / ison
fn = (
0
if fn == 0
else np.format_float_positional(100 * fn, precision=1, min_digits=1) + "%"
)
if fn == "0.0%":
fn = "<0.05%"
fp = int(np.sum(fpv)) / ison
fp = (
0
if fp == 0
else np.format_float_positional(100 * fp, precision=1, min_digits=1) + "%"
)
if fp == "0.0%":
fp = "<0.05%"
fs = int(np.sum(fsv)) / ison
fs = (
0
if fs == 0
else np.format_float_positional(100 * fs, precision=1, min_digits=1) + "%"
)
if fs == "0.0%":
fs = "<0.05%"
err_v = np.mean(
((my_f < level) != (f < level))
| ((my_f > level) != (f > level))
| (np.isnan(my_f) != np.isnan(f))
)
err_v = (
0
if err_v == 0
else np.format_float_positional(100 * err_v, precision=1, min_digits=1) + "%"
)
if err_v == "0.0%":
err_v = "<0.05%"
corr = None if corr is None else compute_corrections_percentage(my_f, corr)
corr = (
""
if corr is None
else (
0
if corr == 0
else np.format_float_positional(100 * corr, precision=1, min_digits=1) + "%"
)
)
if corr == "0.0%":
corr = "<0.05%"
return pd.DataFrame(
{
"Compressor": [title[0]],
"Safeguarded": [title[1]],
"Corrections": [title[2]],
rf"$L_{{\infty}}(\hat{{{var}}})$": [
f"{err_inf:.04}".replace("nan", "NaN")
+ ("" if np.isfinite(err_inf) else f" [{err_inf_fin:.04}]")
],
rf"$L_{{2}}(\hat{{{var}}})$": [
f"{err_2:.04}".replace("nan", "NaN")
+ ("" if np.isfinite(err_2) else f" [{err_2_fin:.04}]")
],
"FN": [fn],
"FP": [fp],
"FS": [fs],
"V": [err_v],
"C": [corr],
"CR": [
rf"$\times$ {np.round(cr, 2)}",
],
}
)
def table_isosurface(
my_f,
f,
cr,
title,
var,
level,
corr,
) -> pd.DataFrame:
with np.errstate(over="ignore", invalid="ignore"):
err_diff = my_f - f
err_diff[np.isnan(my_f) & np.isnan(f)] = 0
err_inf = np.amax(np.abs(err_diff))
err_inf_fin = np.nanmax(np.abs(err_diff))
err_2 = np.sqrt(np.mean(np.square(err_diff)))
err_2_fin = np.sqrt(np.nanmean(np.square(err_diff)))
isof, fnv, fpv, fsv = compute_failures(f, my_f, level)
ison = int(np.sum(isof))
fn = int(np.sum(fnv)) / ison
fn = (
0
if fn == 0
else np.format_float_positional(100 * fn, precision=1, min_digits=1) + "%"
)
if fn == "0.0%":
fn = "<0.05%"
fp = int(np.sum(fpv)) / ison
fp = (
0
if fp == 0
else np.format_float_positional(100 * fp, precision=1, min_digits=1) + "%"
)
if fp == "0.0%":
fp = "<0.05%"
fs = int(np.sum(fsv)) / ison
fs = (
0
if fs == 0
else np.format_float_positional(100 * fs, precision=1, min_digits=1) + "%"
)
if fs == "0.0%":
fs = "<0.05%"
err_v = np.mean(
((my_f < level) != (f < level))
| ((my_f > level) != (f > level))
| (np.isnan(my_f) != np.isnan(f))
)
err_v = (
0
if err_v == 0
else np.format_float_positional(100 * err_v, precision=1, min_digits=1) + "%"
)
if err_v == "0.0%":
err_v = "<0.05%"
corr = None if corr is None else compute_corrections_percentage(my_f, corr)
corr = (
""
if corr is None
else (
0
if corr == 0
else np.format_float_positional(100 * corr, precision=1, min_digits=1) + "%"
)
)
if corr == "0.0%":
corr = "<0.05%"
return pd.DataFrame(
{
"Compressor": [title[0]],
"Safeguarded": [title[1]],
"Corrections": [title[2]],
rf"$L_{{\infty}}(\hat{{{var}}})$": [
f"{err_inf:.04}".replace("nan", "NaN")
+ ("" if np.isfinite(err_inf) else f" [{err_inf_fin:.04}]")
],
rf"$L_{{2}}(\hat{{{var}}})$": [
f"{err_2:.04}".replace("nan", "NaN")
+ ("" if np.isfinite(err_2) else f" [{err_2_fin:.04}]")
],
"FN": [fn],
"FP": [fp],
"FS": [fs],
"V": [err_v],
"C": [corr],
"CR": [
rf"$\times$ {np.round(cr, 2)}",
],
}
)
Copied!
import observe
observations = []
import observe
observations = []
Example 1: Pressure anomaly¶
Lossless compression¶
We first compress the data losslessly with ZStandard at level 22, which gives maximum compression, to provide a baseline.
Copied!
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
with observe.observe(zstd, observations):
pf48_zstd_enc = zstd.encode(pf48)
pf48_zstd = zstd.decode(pf48_zstd_enc)
pf48_zstd_cr = pf48.nbytes / pf48_zstd_enc.nbytes
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
with observe.observe(zstd, observations):
pf48_zstd_enc = zstd.encode(pf48)
pf48_zstd = zstd.decode(pf48_zstd_enc)
pf48_zstd_cr = pf48.nbytes / pf48_zstd_enc.nbytes
Compressing p with lossy compressors¶
We configure each compressor with an absolute error bound of 50 Pa and aim to preserve the isosurface at 0 Pa anomaly. The error bound is chosen to be quite high so that compression artefacts are visually distinguishable.
Since SPERR does not support NaN values, we first replace NaN values with the data mean before applying the SPERR compressor.
Copied!
eb_abs_p = 50
level_p = 0
eb_abs_p = 50
level_p = 0
Copied!
from numcodecs_wasm_zfp import Zfp
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs_p, non_finite="allow-unsafe")
with observe.observe(zfp, observations):
pf48_zfp_enc = zfp.encode(pf48)
pf48_zfp = zfp.decode(pf48_zfp_enc)
pf48_zfp_cr = pf48.nbytes / pf48_zfp_enc.nbytes
from numcodecs_wasm_zfp import Zfp
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs_p, non_finite="allow-unsafe")
with observe.observe(zfp, observations):
pf48_zfp_enc = zfp.encode(pf48)
pf48_zfp = zfp.decode(pf48_zfp_enc)
pf48_zfp_cr = pf48.nbytes / pf48_zfp_enc.nbytes
Copied!
from numcodecs_wasm_sz3 import Sz3
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs_p)
with observe.observe(sz3, observations):
pf48_sz3_enc = sz3.encode(pf48)
pf48_sz3 = sz3.decode(pf48_sz3_enc)
pf48_sz3_cr = pf48.nbytes / pf48_sz3_enc.nbytes
from numcodecs_wasm_sz3 import Sz3
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs_p)
with observe.observe(sz3, observations):
pf48_sz3_enc = sz3.encode(pf48)
pf48_sz3 = sz3.decode(pf48_sz3_enc)
pf48_sz3_cr = pf48.nbytes / pf48_sz3_enc.nbytes
Copied!
from numcodecs_combinators.stack import CodecStack
from numcodecs_replace import ReplaceFilterCodec
from numcodecs_wasm_sperr import Sperr
# inspired by H5Z-SPERR's treatment of NaN values:
# https://github.com/NCAR/H5Z-SPERR/blob/72ebcb00e382886c229c5ef5a7e237fe451d5fb8/src/h5z-sperr.c#L464-L473
# https://github.com/NCAR/H5Z-SPERR/blob/72ebcb00e382886c229c5ef5a7e237fe451d5fb8/src/h5zsperr_helper.cpp#L179-L212
sperr = CodecStack(
ReplaceFilterCodec(replacements={np.nan: "nan_mean"}),
Sperr(mode="pwe", pwe=eb_abs_p),
)
with observe.observe(sperr, observations):
pf48_sperr_enc = sperr.encode(pf48)
pf48_sperr = sperr.decode(pf48_sperr_enc)
pf48_sperr_cr = pf48.nbytes / pf48_sperr_enc.nbytes
from numcodecs_combinators.stack import CodecStack
from numcodecs_replace import ReplaceFilterCodec
from numcodecs_wasm_sperr import Sperr
# inspired by H5Z-SPERR's treatment of NaN values:
# https://github.com/NCAR/H5Z-SPERR/blob/72ebcb00e382886c229c5ef5a7e237fe451d5fb8/src/h5z-sperr.c#L464-L473
# https://github.com/NCAR/H5Z-SPERR/blob/72ebcb00e382886c229c5ef5a7e237fe451d5fb8/src/h5zsperr_helper.cpp#L179-L212
sperr = CodecStack(
ReplaceFilterCodec(replacements={np.nan: "nan_mean"}),
Sperr(mode="pwe", pwe=eb_abs_p),
)
with observe.observe(sperr, observations):
pf48_sperr_enc = sperr.encode(pf48)
pf48_sperr = sperr.decode(pf48_sperr_enc)
pf48_sperr_cr = pf48.nbytes / pf48_sperr_enc.nbytes
Copied!
from numcodecs_zero import ZeroCodec
zero = ZeroCodec()
with observe.observe(zero, observations):
pf48_zero_enc = zero.encode(pf48)
pf48_zero = zero.decode(pf48_zero_enc)
from numcodecs_zero import ZeroCodec
zero = ZeroCodec()
with observe.observe(zero, observations):
pf48_zero_enc = zero.encode(pf48)
pf48_zero = zero.decode(pf48_zero_enc)
Compressing p using the safeguarded lossy compressors¶
We configure the safeguards to bound the pointwise absolute error and preserve the 0 Pa anomaly isosurface by using a sign safeguard that is offset by the value of the isosurface.
Copied!
from numcodecs_safeguards import SafeguardedCodec
pf48_sg = dict()
pf48_sg_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_p, equal_nan=True),
dict(kind="sign", offset=level_p),
],
)
with observe.observe(sg, observations):
pf48_sg_enc = sg.encode(pf48)
pf48_sg[codec_id] = sg.decode(pf48_sg_enc)
pf48_sg_cr[codec_id] = pf48.nbytes / np.asarray(pf48_sg_enc).nbytes
from numcodecs_safeguards import SafeguardedCodec
pf48_sg = dict()
pf48_sg_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_p, equal_nan=True),
dict(kind="sign", offset=level_p),
],
)
with observe.observe(sg, observations):
pf48_sg_enc = sg.encode(pf48)
pf48_sg[codec_id] = sg.decode(pf48_sg_enc)
pf48_sg_cr[codec_id] = pf48.nbytes / np.asarray(pf48_sg_enc).nbytes
Copied!
from numcodecs_safeguards import SafeguardedCodec
pf48_sg_lossless = dict()
pf48_sg_lossless_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_p, equal_nan=True),
dict(kind="sign", offset=level_p),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg, observations):
pf48_sg_lossless_enc = sg.encode(pf48)
pf48_sg_lossless[codec_id] = sg.decode(pf48_sg_lossless_enc)
pf48_sg_lossless_cr[codec_id] = (
pf48.nbytes / np.asarray(pf48_sg_lossless_enc).nbytes
)
from numcodecs_safeguards import SafeguardedCodec
pf48_sg_lossless = dict()
pf48_sg_lossless_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_p, equal_nan=True),
dict(kind="sign", offset=level_p),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg, observations):
pf48_sg_lossless_enc = sg.encode(pf48)
pf48_sg_lossless[codec_id] = sg.decode(pf48_sg_lossless_enc)
pf48_sg_lossless_cr[codec_id] = (
pf48.nbytes / np.asarray(pf48_sg_lossless_enc).nbytes
)
We also compare the outcomes of only safeguarding the isosurface (without any error bound) and of preserving all isosurfaces (that we plot) (again without any error bound).
Copied!
sg_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=level_p),
],
)
with observe.observe(sg_iso, observations):
pf48_sg_iso_enc = sg_iso.encode(pf48)
pf48_sg_iso = sg_iso.decode(pf48_sg_iso_enc)
pf48_sg_iso_cr = pf48.nbytes / np.asarray(pf48_sg_iso_enc).nbytes
sg_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=level_p),
],
)
with observe.observe(sg_iso, observations):
pf48_sg_iso_enc = sg_iso.encode(pf48)
pf48_sg_iso = sg_iso.decode(pf48_sg_iso_enc)
pf48_sg_iso_cr = pf48.nbytes / np.asarray(pf48_sg_iso_enc).nbytes
Copied!
sg_iso_iso = SafeguardedCodec(
codec=zero,
safeguards=[dict(kind="sign", offset=lvl) for lvl in np.linspace(-3500, 3500, 15)],
)
with observe.observe(sg_iso_iso, observations):
pf48_sg_iso_iso_enc = sg_iso_iso.encode(pf48)
pf48_sg_iso_iso = sg_iso_iso.decode(pf48_sg_iso_iso_enc)
pf48_sg_iso_iso_cr = pf48.nbytes / np.asarray(pf48_sg_iso_iso_enc).nbytes
sg_iso_iso = SafeguardedCodec(
codec=zero,
safeguards=[dict(kind="sign", offset=lvl) for lvl in np.linspace(-3500, 3500, 15)],
)
with observe.observe(sg_iso_iso, observations):
pf48_sg_iso_iso_enc = sg_iso_iso.encode(pf48)
pf48_sg_iso_iso = sg_iso_iso.decode(pf48_sg_iso_iso_enc)
pf48_sg_iso_iso_cr = pf48.nbytes / np.asarray(pf48_sg_iso_iso_enc).nbytes
Compressing p with OptZConfig¶
We configure OptZConfig with a custom safety violations metric, implemented in Python, that computes the percentage of violations $V$. We then maximise the score
$$ \textrm{score} = \begin{cases} -\textrm{V} \quad &\text{if } \textrm{V} > 0 \\ \textrm{CR} \quad &\text{otherwise} \end{cases} $$
using the FRAZ search algorithm with 25 iterations, where CR is the achieved compression ratio. Since FRAZ seems to struggle with finding sufficient absolute error bounds spread across several orders of magnitude, we search for bounds in logarithmic space by wrapping each codec in an Exponential<CODEC> meta-codec.
Copied!
import numcodecs
class SafetyViolationsMetric(numcodecs.abc.Codec):
codec_id = "safety-violations-metric"
def __init__(self, level: float):
self._data = None
self._level = level
def encode(self, buf):
# store the original data for later
self._data = np.array(buf, copy=True)
# return no metric
return np.empty(0, dtype=np.float64)
def decode(self, buf, out=None):
# compute the violations
violations = np.mean(
((buf < self._level) != (self._data < self._level))
| ((buf > self._level) != (self._data > self._level))
| (np.isnan(buf) != np.isnan(self._data))
)
self._data = None
# return the violations score metric
return numcodecs.compat.ndarray_copy(np.float64(violations), out)
def get_config(self):
return dict(id=type(self).codec_id, level=self._level)
numcodecs.registry.register_codec(SafetyViolationsMetric)
import numcodecs
class SafetyViolationsMetric(numcodecs.abc.Codec):
codec_id = "safety-violations-metric"
def __init__(self, level: float):
self._data = None
self._level = level
def encode(self, buf):
# store the original data for later
self._data = np.array(buf, copy=True)
# return no metric
return np.empty(0, dtype=np.float64)
def decode(self, buf, out=None):
# compute the violations
violations = np.mean(
((buf < self._level) != (self._data < self._level))
| ((buf > self._level) != (self._data > self._level))
| (np.isnan(buf) != np.isnan(self._data))
)
self._data = None
# return the violations score metric
return numcodecs.compat.ndarray_copy(np.float64(violations), out)
def get_config(self):
return dict(id=type(self).codec_id, level=self._level)
numcodecs.registry.register_codec(SafetyViolationsMetric)
Copied!
class ExponentialZfp(Zfp):
codec_id = "e-zfp.rs"
def __new__(cls, tolerance: float, **kwargs):
codec = super().__new__(cls, tolerance=np.exp(tolerance), **kwargs)
codec._tolerance = tolerance
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"tolerance": self._tolerance,
}
class ExponentialSz3(Sz3):
codec_id = "e-sz3.rs"
def __new__(cls, eb_abs: float, **kwargs):
codec = super().__new__(cls, eb_abs=np.exp(eb_abs), **kwargs)
codec._eb_abs = eb_abs
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"eb_abs": self._eb_abs,
}
# libpressio does not support the configuration format used by the CodecStack
# and ReplaceFilterCodec, so we inline their functionality here
class ExponentialMaskedSperr(Sperr):
codec_id = "e-sperr.rs"
def __new__(cls, pwe: float, **kwargs):
codec = super().__new__(cls, pwe=np.exp(pwe), **kwargs)
codec._pwe = pwe
codec._filter = ReplaceFilterCodec(replacements={np.nan: "nan_mean"})
return codec
def encode(self, buf):
return super().encode(self._filter.encode(buf))
def decode(self, buf, out=None):
return self._filter.decode(super().decode(buf, out=out), out=out)
def get_config(self):
return {**super().get_config(), "id": type(self).codec_id, "pwe": self._pwe}
numcodecs.registry.register_codec(ExponentialZfp)
numcodecs.registry.register_codec(ExponentialSz3)
numcodecs.registry.register_codec(ExponentialMaskedSperr)
class ExponentialZfp(Zfp):
codec_id = "e-zfp.rs"
def __new__(cls, tolerance: float, **kwargs):
codec = super().__new__(cls, tolerance=np.exp(tolerance), **kwargs)
codec._tolerance = tolerance
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"tolerance": self._tolerance,
}
class ExponentialSz3(Sz3):
codec_id = "e-sz3.rs"
def __new__(cls, eb_abs: float, **kwargs):
codec = super().__new__(cls, eb_abs=np.exp(eb_abs), **kwargs)
codec._eb_abs = eb_abs
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"eb_abs": self._eb_abs,
}
# libpressio does not support the configuration format used by the CodecStack
# and ReplaceFilterCodec, so we inline their functionality here
class ExponentialMaskedSperr(Sperr):
codec_id = "e-sperr.rs"
def __new__(cls, pwe: float, **kwargs):
codec = super().__new__(cls, pwe=np.exp(pwe), **kwargs)
codec._pwe = pwe
codec._filter = ReplaceFilterCodec(replacements={np.nan: "nan_mean"})
return codec
def encode(self, buf):
return super().encode(self._filter.encode(buf))
def decode(self, buf, out=None):
return self._filter.decode(super().decode(buf, out=out), out=out)
def get_config(self):
return {**super().get_config(), "id": type(self).codec_id, "pwe": self._pwe}
numcodecs.registry.register_codec(ExponentialZfp)
numcodecs.registry.register_codec(ExponentialSz3)
numcodecs.registry.register_codec(ExponentialMaskedSperr)
Copied!
from numcodecs_wasm_pressio import Pressio
pf48_optzconfig = dict()
pf48_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-12), # tiny bound
(sz3, "eb_abs", 1e-5), # decent guess
(sperr[-1], "pwe", 1e-12), # tiny bound
]:
optzconfig = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs_p),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
"numcodecs.rs-metric:level": level_p,
},
)
with observe.observe(optzconfig, observations):
pf48_optzconfig_enc = optzconfig.encode(pf48)
pf48_optzconfig[codec.codec_id] = optzconfig.decode(pf48_optzconfig_enc)
pf48_optzconfig_cr[codec.codec_id] = pf48.nbytes / pf48_optzconfig_enc.nbytes
from numcodecs_wasm_pressio import Pressio
pf48_optzconfig = dict()
pf48_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-12), # tiny bound
(sz3, "eb_abs", 1e-5), # decent guess
(sperr[-1], "pwe", 1e-12), # tiny bound
]:
optzconfig = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs_p),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
"numcodecs.rs-metric:level": level_p,
},
)
with observe.observe(optzconfig, observations):
pf48_optzconfig_enc = optzconfig.encode(pf48)
pf48_optzconfig[codec.codec_id] = optzconfig.decode(pf48_optzconfig_enc)
pf48_optzconfig_cr[codec.codec_id] = pf48.nbytes / pf48_optzconfig_enc.nbytes
rank={0,1,} iter={0} input={-11.8595,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={1} input={-20.2729,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={2} input={-3.57471,} output={-0.0041716,} objective={-0.0041716}
rank={0,1,} iter={3} input={-20.9567,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={4} input={1.75846,} output={-0.00429304,} objective={-0.00429304}
rank={0,1,} iter={5} input={-25.5032,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={6} input={0.770514,} output={-0.00423556,} objective={-0.00423556}
rank={0,1,} iter={7} input={1.39129,} output={-0.00429304,} objective={-0.00429304}
rank={0,1,} iter={8} input={-12.1628,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={9} input={-22.62,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={10} input={-2.0315,} output={-0.00417636,} objective={-0.00417636}
rank={0,1,} iter={11} input={-11.9693,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={12} input={-22.9839,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={13} input={-1.14688,} output={-0.00418116,} objective={-0.00418116}
rank={0,1,} iter={14} input={-24.9913,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={15} input={-20.4875,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={16} input={-16.1997,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={17} input={3.48382,} output={-0.0049368,} objective={-0.0049368}
rank={0,1,} iter={18} input={-11.3317,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={19} input={1.20116,} output={-0.00423556,} objective={-0.00423556}
rank={0,1,} iter={20} input={-27.376,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={21} input={-22.5664,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={22} input={-1.93138,} output={-0.00417636,} objective={-0.00417636}
rank={0,1,} iter={23} input={-25.0826,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={24} input={2.97377,} output={-0.00458564,} objective={-0.00458564}
final_iter={25} inputs={-11.8595,} output={-0.00417116,}
rank={0,1,} iter={0} input={-3.80045,} output={-8.52e-06,} objective={-8.52e-06}
rank={0,1,} iter={1} input={-7.91471,} output={-1.6e-07,} objective={-1.6e-07}
rank={0,1,} iter={2} input={0.250916,} output={-0.0003614,} objective={-0.0003614}
rank={0,1,} iter={3} input={-8.24911,} output={-1.6e-07,} objective={-1.6e-07}
rank={0,1,} iter={4} input={2.8589,} output={-0.00150972,} objective={-0.00150972}
rank={0,1,} iter={5} input={-10.4724,} output={1.74198,} objective={1.74198}
rank={0,1,} iter={6} input={-11.5129,} output={1.64105,} objective={1.64105}
rank={0,1,} iter={7} input={-10.9261,} output={1.76158,} objective={1.76158}
rank={0,1,} iter={8} input={-11.1404,} output={1.63799,} objective={1.63799}
rank={0,1,} iter={9} input={-10.7225,} output={1.75837,} objective={1.75837}
rank={0,1,} iter={10} input={-10.6087,} output={1.77314,} objective={1.77314}
rank={0,1,} iter={11} input={-10.8275,} output={1.78306,} objective={1.78306}
rank={0,1,} iter={12} input={-10.8638,} output={1.78995,} objective={1.78995}
rank={0,1,} iter={13} input={-10.7905,} output={1.7698,} objective={1.7698}
rank={0,1,} iter={14} input={-10.5615,} output={1.76306,} objective={1.76306}
rank={0,1,} iter={15} input={-10.6568,} output={1.78142,} objective={1.78142}
rank={0,1,} iter={16} input={-10.877,} output={1.75176,} objective={1.75176}
rank={0,1,} iter={17} input={-10.8472,} output={1.78696,} objective={1.78696}
rank={0,1,} iter={18} input={-10.8741,} output={1.75097,} objective={1.75097}
rank={0,1,} iter={19} input={-10.5601,} output={1.76248,} objective={1.76248}
rank={0,1,} iter={20} input={-10.8741,} output={1.75095,} objective={1.75095}
rank={0,1,} iter={21} input={-10.712,} output={1.75456,} objective={1.75456}
rank={0,1,} iter={22} input={-10.8726,} output={1.75057,} objective={1.75057}
rank={0,1,} iter={23} input={-10.7879,} output={1.76901,} objective={1.76901}
rank={0,1,} iter={24} input={-10.8773,} output={1.75175,} objective={1.75175}
final_iter={25} inputs={-10.8638,} output={1.78995,}
rank={0,1,} iter={0} input={-11.8595,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={1} input={-20.2729,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={2} input={-3.57471,} output={-0.00417668,} objective={-0.00417668}
rank={0,1,} iter={3} input={-20.9567,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={4} input={1.75846,} output={-0.00453244,} objective={-0.00453244}
rank={0,1,} iter={5} input={-25.5032,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={6} input={0.770514,} output={-0.00435444,} objective={-0.00435444}
rank={0,1,} iter={7} input={1.39129,} output={-0.00446008,} objective={-0.00446008}
rank={0,1,} iter={8} input={-12.1628,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={9} input={-22.62,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={10} input={-2.0315,} output={-0.0041916,} objective={-0.0041916}
rank={0,1,} iter={11} input={-11.9693,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={12} input={-22.9839,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={13} input={-1.14688,} output={-0.00421444,} objective={-0.00421444}
rank={0,1,} iter={14} input={-24.9913,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={15} input={-20.4875,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={16} input={-16.1997,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={17} input={3.48382,} output={-0.0050806,} objective={-0.0050806}
rank={0,1,} iter={18} input={-11.3317,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={19} input={1.20116,} output={-0.00442624,} objective={-0.00442624}
rank={0,1,} iter={20} input={-27.376,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={21} input={-22.5664,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={22} input={-1.93138,} output={-0.00419328,} objective={-0.00419328}
rank={0,1,} iter={23} input={-25.0826,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={24} input={2.97377,} output={-0.00490332,} objective={-0.00490332}
final_iter={25} inputs={-11.8595,} output={-0.00417116,}
Visual comparison of the pressure anomaly isosurfaces¶
Copied!
plt.rcParams["xtick.color"] = "#AAAAAA"
plt.rcParams["ytick.color"] = "#AAAAAA"
plt.rcParams["axes.edgecolor"] = "#AAAAAA"
plt.rcParams["xtick.color"] = "#AAAAAA"
plt.rcParams["ytick.color"] = "#AAAAAA"
plt.rcParams["axes.edgecolor"] = "#AAAAAA"
Copied!
fig = plt.figure(figsize=(25, 18))
gs = gridspec.GridSpec(
3, 5, left=0.035, right=0.925, top=0.915, bottom=0.025, wspace=0.3, hspace=0.2
)
plot_isosurface(
pf48,
pf48,
1.0,
fig.add_subplot(gs[0, 0], projection="3d", computed_zorder=False),
"Original",
"p",
"Pressure anomaly (Pa)",
level=level_p,
my_span=0,
span=3500,
)
plot_isosurface(
pf48_zfp,
pf48,
pf48_zfp_cr,
fig.add_subplot(gs[0, 1], projection="3d", computed_zorder=False),
r"ZFP($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=5,
span=3500,
error=True,
)
plot_isosurface(
pf48_sz3,
pf48,
pf48_sz3_cr,
fig.add_subplot(gs[0, 2], projection="3d", computed_zorder=False),
r"SZ3($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
error=True,
)
plot_isosurface(
pf48_sperr,
pf48,
pf48_sperr_cr,
fig.add_subplot(gs[0, 3], projection="3d", computed_zorder=False),
r"SPERR($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
error=True,
)
plot_isosurface(
pf48_sg_iso,
pf48,
pf48_sg_iso_cr,
fig.add_subplot(gs[0, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_sg["zero"],
pf48,
pf48_sg_cr["zero"],
fig.add_subplot(gs[1, 0], projection="3d", computed_zorder=False),
r"Safeguarded(0, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_sg["zfp.rs"],
pf48,
pf48_sg_cr["zfp.rs"],
fig.add_subplot(gs[1, 1], projection="3d", computed_zorder=False),
r"Safeguarded(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_zfp,
error=True,
)
plot_isosurface(
pf48_sg["sz3.rs"],
pf48,
pf48_sg_cr["sz3.rs"],
fig.add_subplot(gs[1, 2], projection="3d", computed_zorder=False),
r"Safeguarded(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_sz3,
error=True,
)
plot_isosurface(
pf48_sg["sperr.rs"],
pf48,
pf48_sg_cr["sperr.rs"],
fig.add_subplot(gs[1, 3], projection="3d", computed_zorder=False),
r"Safeguarded(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_sperr,
error=True,
)
plot_isosurface(
pf48_sg_iso_iso,
pf48,
pf48_sg_iso_iso_cr,
fig.add_subplot(gs[1, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso*)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_optzconfig["zfp.rs"],
pf48,
pf48_optzconfig_cr["zfp.rs"],
fig.add_subplot(gs[2, 1], projection="3d", computed_zorder=False),
r"OptZConfig(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
plot_isosurface(
pf48_optzconfig["sz3.rs"],
pf48,
pf48_optzconfig_cr["sz3.rs"],
fig.add_subplot(gs[2, 2], projection="3d", computed_zorder=False),
r"OptZConfig(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
plot_isosurface(
pf48_optzconfig["sperr.rs"],
pf48,
pf48_optzconfig_cr["sperr.rs"],
fig.add_subplot(gs[2, 3], projection="3d", computed_zorder=False),
r"OptZConfig(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
# plt.tight_layout()
plt.savefig(Path("plots") / "isosurface-p.pdf", dpi=300)
plt.show()
fig = plt.figure(figsize=(25, 18))
gs = gridspec.GridSpec(
3, 5, left=0.035, right=0.925, top=0.915, bottom=0.025, wspace=0.3, hspace=0.2
)
plot_isosurface(
pf48,
pf48,
1.0,
fig.add_subplot(gs[0, 0], projection="3d", computed_zorder=False),
"Original",
"p",
"Pressure anomaly (Pa)",
level=level_p,
my_span=0,
span=3500,
)
plot_isosurface(
pf48_zfp,
pf48,
pf48_zfp_cr,
fig.add_subplot(gs[0, 1], projection="3d", computed_zorder=False),
r"ZFP($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=5,
span=3500,
error=True,
)
plot_isosurface(
pf48_sz3,
pf48,
pf48_sz3_cr,
fig.add_subplot(gs[0, 2], projection="3d", computed_zorder=False),
r"SZ3($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
error=True,
)
plot_isosurface(
pf48_sperr,
pf48,
pf48_sperr_cr,
fig.add_subplot(gs[0, 3], projection="3d", computed_zorder=False),
r"SPERR($\epsilon_{abs}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
error=True,
)
plot_isosurface(
pf48_sg_iso,
pf48,
pf48_sg_iso_cr,
fig.add_subplot(gs[0, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_sg["zero"],
pf48,
pf48_sg_cr["zero"],
fig.add_subplot(gs[1, 0], projection="3d", computed_zorder=False),
r"Safeguarded(0, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_sg["zfp.rs"],
pf48,
pf48_sg_cr["zfp.rs"],
fig.add_subplot(gs[1, 1], projection="3d", computed_zorder=False),
r"Safeguarded(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_zfp,
error=True,
)
plot_isosurface(
pf48_sg["sz3.rs"],
pf48,
pf48_sg_cr["sz3.rs"],
fig.add_subplot(gs[1, 2], projection="3d", computed_zorder=False),
r"Safeguarded(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_sz3,
error=True,
)
plot_isosurface(
pf48_sg["sperr.rs"],
pf48,
pf48_sg_cr["sperr.rs"],
fig.add_subplot(gs[1, 3], projection="3d", computed_zorder=False),
r"Safeguarded(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
corr=pf48_sperr,
error=True,
)
plot_isosurface(
pf48_sg_iso_iso,
pf48,
pf48_sg_iso_iso_cr,
fig.add_subplot(gs[1, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso*)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=50,
span=3500,
corr=pf48_zero,
error=True,
)
plot_isosurface(
pf48_optzconfig["zfp.rs"],
pf48,
pf48_optzconfig_cr["zfp.rs"],
fig.add_subplot(gs[2, 1], projection="3d", computed_zorder=False),
r"OptZConfig(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
plot_isosurface(
pf48_optzconfig["sz3.rs"],
pf48,
pf48_optzconfig_cr["sz3.rs"],
fig.add_subplot(gs[2, 2], projection="3d", computed_zorder=False),
r"OptZConfig(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
plot_isosurface(
pf48_optzconfig["sperr.rs"],
pf48,
pf48_optzconfig_cr["sperr.rs"],
fig.add_subplot(gs[2, 3], projection="3d", computed_zorder=False),
r"OptZConfig(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"p",
"Absolute error over pressure anomaly (Pa)",
level=level_p,
my_span=eb_abs_p,
span=3500,
error=True,
)
# plt.tight_layout()
plt.savefig(Path("plots") / "isosurface-p.pdf", dpi=300)
plt.show()
Copied!
iso_p_table = pd.concat(
[
table_isosurface(
pf48_sg_lossless["zero"],
pf48,
pf48_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg["zero"],
pf48,
pf48_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg_iso,
pf48,
pf48_sg_iso_cr,
["0", r"$\text{iso}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg_iso_iso,
pf48,
pf48_sg_iso_iso_cr,
["0", r"$\text{iso*}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_zfp,
pf48,
pf48_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["zfp.rs"],
pf48,
pf48_sg_lossless_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_zfp,
),
table_isosurface(
pf48_sg["zfp.rs"],
pf48,
pf48_sg_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_zfp,
),
table_isosurface(
pf48_optzconfig["zfp.rs"],
pf48,
pf48_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sz3,
pf48,
pf48_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["sz3.rs"],
pf48,
pf48_sg_lossless_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_sz3,
),
table_isosurface(
pf48_sg["sz3.rs"],
pf48,
pf48_sg_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_sz3,
),
table_isosurface(
pf48_optzconfig["sz3.rs"],
pf48,
pf48_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sperr,
pf48,
pf48_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["sperr.rs"],
pf48,
pf48_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"lossless",
],
"p",
level_p,
pf48_sperr,
),
table_isosurface(
pf48_sg["sperr.rs"],
pf48,
pf48_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"one-shot",
],
"p",
level_p,
pf48_sperr,
),
table_isosurface(
pf48_optzconfig["sperr.rs"],
pf48,
pf48_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_zstd,
pf48,
pf48_zstd_cr,
["ZSTD(20)", "-", ""],
"p",
level_p,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections"])
Path("tables").joinpath("isosurface-p.tex").write_text(
iso_p_table.to_latex(escape=False)
.replace("%", r"\%")
.replace("\\cline{1-11} \\cline{2-11}\n\\bottomrule", "\\bottomrule")
)
iso_p_table
iso_p_table = pd.concat(
[
table_isosurface(
pf48_sg_lossless["zero"],
pf48,
pf48_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg["zero"],
pf48,
pf48_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg_iso,
pf48,
pf48_sg_iso_cr,
["0", r"$\text{iso}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_sg_iso_iso,
pf48,
pf48_sg_iso_iso_cr,
["0", r"$\text{iso*}$", "one-shot"],
"p",
level_p,
pf48_zero,
),
table_isosurface(
pf48_zfp,
pf48,
pf48_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["zfp.rs"],
pf48,
pf48_sg_lossless_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_zfp,
),
table_isosurface(
pf48_sg["zfp.rs"],
pf48,
pf48_sg_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_zfp,
),
table_isosurface(
pf48_optzconfig["zfp.rs"],
pf48,
pf48_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sz3,
pf48,
pf48_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["sz3.rs"],
pf48,
pf48_sg_lossless_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"p",
level_p,
pf48_sz3,
),
table_isosurface(
pf48_sg["sz3.rs"],
pf48,
pf48_sg_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"p",
level_p,
pf48_sz3,
),
table_isosurface(
pf48_optzconfig["sz3.rs"],
pf48,
pf48_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sperr,
pf48,
pf48_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_sg_lossless["sperr.rs"],
pf48,
pf48_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"lossless",
],
"p",
level_p,
pf48_sperr,
),
table_isosurface(
pf48_sg["sperr.rs"],
pf48,
pf48_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"one-shot",
],
"p",
level_p,
pf48_sperr,
),
table_isosurface(
pf48_optzconfig["sperr.rs"],
pf48,
pf48_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"p",
level_p,
None,
),
table_isosurface(
pf48_zstd,
pf48,
pf48_zstd_cr,
["ZSTD(20)", "-", ""],
"p",
level_p,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections"])
Path("tables").joinpath("isosurface-p.tex").write_text(
iso_p_table.to_latex(escape=False)
.replace("%", r"\%")
.replace("\\cline{1-11} \\cline{2-11}\n\\bottomrule", "\\bottomrule")
)
iso_p_table
| $L_{\infty}(\hat{p})$ | $L_{2}(\hat{p})$ | FN | FP | FS | V | C | CR | |||
|---|---|---|---|---|---|---|---|---|---|---|
| Compressor | Safeguarded | Corrections | ||||||||
| 0 | $\epsilon_{abs} \cup \text{iso}$ | lossless | 0.0 | 0.0 | 0 | 0 | 0 | 0 | 100.0% | $\times$ 1.59 |
| one-shot | 50.0 | 27.38 | 0 | 0 | 0 | 0 | 100.0% | $\times$ 74.32 | ||
| $\text{iso}$ | one-shot | 3.412e+03 | 629.2 | 0 | 0 | 0 | 0 | 100.0% | $\times$ 3549.5 | |
| $\text{iso*}$ | one-shot | 500.0 | 226.2 | 0 | 0 | 0 | 0 | 100.0% | $\times$ 456.07 | |
| ZFP($\epsilon_{abs}$) | - | NaN [12.41] | NaN [1.125] | 11.2% | 103.7% | 43.1% | 0.5% | $\times$ 25.77 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 12.41 | 1.122 | 0 | 0 | 0 | 0 | 0.5% | $\times$ 24.29 | |
| one-shot | 47.9 | 1.127 | 0 | 0 | 0 | 0 | 0.5% | $\times$ 25.09 | ||
| OptZConfig(ZFP) | $\epsilon_{abs} \cup \text{iso}$ | NaN [0.0001221] | NaN [1.195e-06] | 0 | 73.6% | 13.1% | 0.4% | $\times$ 1.6 | ||
| SZ3($\epsilon_{abs}$) | - | 50.0 | 11.01 | 21.2% | 29.5% | 68.2% | 0.2% | $\times$ 478.31 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 50.0 | 10.99 | 0 | 0 | 0 | 0 | 0.2% | $\times$ 196.09 | |
| one-shot | 50.0 | 11.0 | 0 | 0 | 0 | 0 | 0.2% | $\times$ 191.61 | ||
| OptZConfig(SZ3) | $\epsilon_{abs} \cup \text{iso}$ | 1.913e-05 | 7.603e-06 | 0 | 0 | 0 | 0 | $\times$ 1.79 | ||
| SPERR($\epsilon_{abs}$) | - | NaN [50.0] | NaN [3.929] | 12.4% | 11.4% | 57.1% | 0.5% | $\times$ 438.15 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 50.0 | 3.913 | 0 | 0 | 0 | 0 | 0.5% | $\times$ 137.77 | |
| one-shot | 50.0 | 3.922 | 0 | 0 | 0 | 0 | 0.5% | $\times$ 248.77 | ||
| OptZConfig(SPERR) | $\epsilon_{abs} \cup \text{iso}$ | NaN [7.629e-06] | NaN [1.839e-06] | 0 | 0 | 0 | 0.4% | $\times$ 2.25 | ||
| ZSTD(20) | - | 0.0 | 0.0 | 0 | 0 | 0 | 0 | $\times$ 1.17 |
Example 2: u wind¶
Lossless compression¶
We first compress the data losslessly with ZStandard at level 22, which gives maximum compression, to provide a baseline.
Copied!
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
with observe.observe(zstd, observations):
uf48_zstd_enc = zstd.encode(uf48)
uf48_zstd = zstd.decode(uf48_zstd_enc)
uf48_zstd_cr = uf48.nbytes / uf48_zstd_enc.nbytes
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
with observe.observe(zstd, observations):
uf48_zstd_enc = zstd.encode(uf48)
uf48_zstd = zstd.decode(uf48_zstd_enc)
uf48_zstd_cr = uf48.nbytes / uf48_zstd_enc.nbytes
Compressing u with lossy compressors¶
We configure each compressor with an absolute error bound of 5 m/s and aim to preserve the mean value isosurface. The error bound is chosen to be quite high so that compression artefacts are visually distinguishable.
Copied!
eb_abs_u = 5
level_u = float(np.nanmean(uf48))
eb_abs_u = 5
level_u = float(np.nanmean(uf48))
Copied!
from numcodecs_wasm_zfp import Zfp
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs_u, non_finite="allow-unsafe")
with observe.observe(zfp, observations):
uf48_zfp_enc = zfp.encode(uf48)
uf48_zfp = zfp.decode(uf48_zfp_enc)
uf48_zfp_cr = uf48.nbytes / uf48_zfp_enc.nbytes
from numcodecs_wasm_zfp import Zfp
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs_u, non_finite="allow-unsafe")
with observe.observe(zfp, observations):
uf48_zfp_enc = zfp.encode(uf48)
uf48_zfp = zfp.decode(uf48_zfp_enc)
uf48_zfp_cr = uf48.nbytes / uf48_zfp_enc.nbytes
Copied!
from numcodecs_wasm_sz3 import Sz3
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs_u)
with observe.observe(sz3, observations):
uf48_sz3_enc = sz3.encode(uf48)
uf48_sz3 = sz3.decode(uf48_sz3_enc)
uf48_sz3_cr = uf48.nbytes / uf48_sz3_enc.nbytes
from numcodecs_wasm_sz3 import Sz3
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs_u)
with observe.observe(sz3, observations):
uf48_sz3_enc = sz3.encode(uf48)
uf48_sz3 = sz3.decode(uf48_sz3_enc)
uf48_sz3_cr = uf48.nbytes / uf48_sz3_enc.nbytes
Copied!
from numcodecs_combinators.stack import CodecStack
from numcodecs_replace import ReplaceFilterCodec
from numcodecs_wasm_sperr import Sperr
sperr = CodecStack(
ReplaceFilterCodec(replacements={np.nan: "nan_mean"}),
Sperr(mode="pwe", pwe=eb_abs_u),
)
with observe.observe(sperr, observations):
uf48_sperr_enc = sperr.encode(uf48)
uf48_sperr = sperr.decode(uf48_sperr_enc)
uf48_sperr_cr = uf48.nbytes / uf48_sperr_enc.nbytes
from numcodecs_combinators.stack import CodecStack
from numcodecs_replace import ReplaceFilterCodec
from numcodecs_wasm_sperr import Sperr
sperr = CodecStack(
ReplaceFilterCodec(replacements={np.nan: "nan_mean"}),
Sperr(mode="pwe", pwe=eb_abs_u),
)
with observe.observe(sperr, observations):
uf48_sperr_enc = sperr.encode(uf48)
uf48_sperr = sperr.decode(uf48_sperr_enc)
uf48_sperr_cr = uf48.nbytes / uf48_sperr_enc.nbytes
Copied!
from numcodecs_zero import ZeroCodec
zero = ZeroCodec()
with observe.observe(zero, observations):
uf48_zero_enc = zero.encode(uf48)
uf48_zero = zero.decode(uf48_zero_enc)
from numcodecs_zero import ZeroCodec
zero = ZeroCodec()
with observe.observe(zero, observations):
uf48_zero_enc = zero.encode(uf48)
uf48_zero = zero.decode(uf48_zero_enc)
Compressing u using the safeguarded lossy compressors¶
We configure the safeguards to bound the pointwise absolute error and preserve the mean u wind isosurface by using a sign safeguard that is offset by the value of the isosurface.
Copied!
from numcodecs_safeguards import SafeguardedCodec
uf48_sg = dict()
uf48_sg_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_u, equal_nan=True),
dict(kind="sign", offset=level_u),
],
)
with observe.observe(sg, observations):
uf48_sg_enc = sg.encode(uf48)
uf48_sg[codec_id] = sg.decode(uf48_sg_enc)
uf48_sg_cr[codec_id] = uf48.nbytes / np.asarray(uf48_sg_enc).nbytes
from numcodecs_safeguards import SafeguardedCodec
uf48_sg = dict()
uf48_sg_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_u, equal_nan=True),
dict(kind="sign", offset=level_u),
],
)
with observe.observe(sg, observations):
uf48_sg_enc = sg.encode(uf48)
uf48_sg[codec_id] = sg.decode(uf48_sg_enc)
uf48_sg_cr[codec_id] = uf48.nbytes / np.asarray(uf48_sg_enc).nbytes
Copied!
from numcodecs_safeguards import SafeguardedCodec
uf48_sg_lossless = dict()
uf48_sg_lossless_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_u, equal_nan=True),
dict(kind="sign", offset=level_u),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg, observations):
uf48_sg_lossless_enc = sg.encode(uf48)
uf48_sg_lossless[codec_id] = sg.decode(uf48_sg_lossless_enc)
uf48_sg_lossless_cr[codec_id] = (
uf48.nbytes / np.asarray(uf48_sg_lossless_enc).nbytes
)
from numcodecs_safeguards import SafeguardedCodec
uf48_sg_lossless = dict()
uf48_sg_lossless_cr = dict()
for codec_id, codec in {
zero.codec_id: zero,
zfp.codec_id: zfp,
sz3.codec_id: sz3,
"sperr.rs": sperr,
}.items():
sg = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs_u, equal_nan=True),
dict(kind="sign", offset=level_u),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg, observations):
uf48_sg_lossless_enc = sg.encode(uf48)
uf48_sg_lossless[codec_id] = sg.decode(uf48_sg_lossless_enc)
uf48_sg_lossless_cr[codec_id] = (
uf48.nbytes / np.asarray(uf48_sg_lossless_enc).nbytes
)
Copied!
sg_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=level_u),
],
)
with observe.observe(sg_iso, observations):
uf48_sg_iso_enc = sg_iso.encode(uf48)
uf48_sg_iso = sg_iso.decode(uf48_sg_iso_enc)
uf48_sg_iso_cr = uf48.nbytes / np.asarray(uf48_sg_iso_enc).nbytes
sg_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=level_u),
],
)
with observe.observe(sg_iso, observations):
uf48_sg_iso_enc = sg_iso.encode(uf48)
uf48_sg_iso = sg_iso.decode(uf48_sg_iso_enc)
uf48_sg_iso_cr = uf48.nbytes / np.asarray(uf48_sg_iso_enc).nbytes
Copied!
sg_iso_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=float(uf48.dtype.type(lvl)))
for lvl in np.linspace(-np.abs(level_u) * 15, np.abs(level_u) * 15, 16)
],
)
with observe.observe(sg_iso_iso, observations):
uf48_sg_iso_iso_enc = sg_iso_iso.encode(uf48)
uf48_sg_iso_iso = sg_iso_iso.decode(uf48_sg_iso_iso_enc)
uf48_sg_iso_iso_cr = uf48.nbytes / np.asarray(uf48_sg_iso_iso_enc).nbytes
sg_iso_iso = SafeguardedCodec(
codec=zero,
safeguards=[
dict(kind="sign", offset=float(uf48.dtype.type(lvl)))
for lvl in np.linspace(-np.abs(level_u) * 15, np.abs(level_u) * 15, 16)
],
)
with observe.observe(sg_iso_iso, observations):
uf48_sg_iso_iso_enc = sg_iso_iso.encode(uf48)
uf48_sg_iso_iso = sg_iso_iso.decode(uf48_sg_iso_iso_enc)
uf48_sg_iso_iso_cr = uf48.nbytes / np.asarray(uf48_sg_iso_iso_enc).nbytes
Compressing u with OptZConfig¶
Copied!
from numcodecs_wasm_pressio import Pressio
uf48_optzconfig = dict()
uf48_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-12), # tiny bound
(sz3, "eb_abs", 1e-8), # decent guess
(sperr[-1], "pwe", 1e-12), # tiny bound
]:
optzconfig = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs_u),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
"numcodecs.rs-metric:level": level_u,
},
)
with observe.observe(optzconfig, observations):
uf48_optzconfig_enc = optzconfig.encode(uf48)
uf48_optzconfig[codec.codec_id] = optzconfig.decode(uf48_optzconfig_enc)
uf48_optzconfig_cr[codec.codec_id] = uf48.nbytes / uf48_optzconfig_enc.nbytes
from numcodecs_wasm_pressio import Pressio
uf48_optzconfig = dict()
uf48_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-12), # tiny bound
(sz3, "eb_abs", 1e-8), # decent guess
(sperr[-1], "pwe", 1e-12), # tiny bound
]:
optzconfig = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs_u),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
"numcodecs.rs-metric:level": level_u,
},
)
with observe.observe(optzconfig, observations):
uf48_optzconfig_enc = optzconfig.encode(uf48)
uf48_optzconfig[codec.codec_id] = optzconfig.decode(uf48_optzconfig_enc)
uf48_optzconfig_cr[codec.codec_id] = uf48.nbytes / uf48_optzconfig_enc.nbytes
rank={0,1,} iter={0} input={-13.0108,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={1} input={-20.81,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={2} input={-5.33078,} output={-0.00418524,} objective={-0.00418524}
rank={0,1,} iter={3} input={-21.4439,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={4} input={-0.386922,} output={-0.00540696,} objective={-0.00540696}
rank={0,1,} iter={5} input={-25.6586,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={6} input={-1.30275,} output={-0.004854,} objective={-0.004854}
rank={0,1,} iter={7} input={-0.727284,} output={-0.004854,} objective={-0.004854}
rank={0,1,} iter={8} input={-13.2919,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={9} input={-22.9858,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={10} input={-3.90022,} output={-0.00422184,} objective={-0.00422184}
rank={0,1,} iter={11} input={-13.1126,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={12} input={-23.3232,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={13} input={-3.08018,} output={-0.00427656,} objective={-0.00427656}
rank={0,1,} iter={14} input={-25.184,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={15} input={-21.0089,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={16} input={-17.0342,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={17} input={1.2125,} output={-0.00784632,} objective={-0.00784632}
rank={0,1,} iter={18} input={-12.5215,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={19} input={-0.903541,} output={-0.004854,} objective={-0.004854}
rank={0,1,} iter={20} input={-27.3946,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={21} input={-22.9361,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={22} input={-3.80741,} output={-0.00422184,} objective={-0.00422184}
rank={0,1,} iter={23} input={-25.2687,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={24} input={0.739675,} output={-0.00784632,} objective={-0.00784632}
final_iter={25} inputs={-13.0108,} output={-0.00417116,}
rank={0,1,} iter={0} input={-8.40562,} output={-7.36e-06,} objective={-7.36e-06}
rank={0,1,} iter={1} input={-13.7482,} output={1.66407,} objective={1.66407}
rank={0,1,} iter={2} input={-3.1447,} output={-0.00117916,} objective={-0.00117916}
rank={0,1,} iter={3} input={-18.4207,} output={1.38155,} objective={1.38155}
rank={0,1,} iter={4} input={-15.6296,} output={1.44239,} objective={1.44239}
rank={0,1,} iter={5} input={-8.71154,} output={-4.84e-06,} objective={-4.84e-06}
rank={0,1,} iter={6} input={-16.929,} output={1.38817,} objective={1.38817}
rank={0,1,} iter={7} input={-11.2299,} output={-6e-07,} objective={-6e-07}
rank={0,1,} iter={8} input={1.60908,} output={-0.0473198,} objective={-0.0473198}
rank={0,1,} iter={9} input={-14.356,} output={1.48442,} objective={1.48442}
rank={0,1,} iter={10} input={-14.9655,} output={1.50573,} objective={1.50573}
rank={0,1,} iter={11} input={-13.1186,} output={1.84762,} objective={1.84762}
rank={0,1,} iter={12} input={-5.77793,} output={-8.916e-05,} objective={-8.916e-05}
rank={0,1,} iter={13} input={-12.489,} output={-1.6e-07,} objective={-1.6e-07}
rank={0,1,} iter={14} input={-0.768595,} output={-0.0077354,} objective={-0.0077354}
rank={0,1,} iter={15} input={-13.3765,} output={1.72969,} objective={1.72969}
rank={0,1,} iter={16} input={-4.45948,} output={-0.00033684,} objective={-0.00033684}
rank={0,1,} iter={17} input={-12.9612,} output={-1.6e-07,} objective={-1.6e-07}
rank={0,1,} iter={18} input={-7.09037,} output={-2.476e-05,} objective={-2.476e-05}
rank={0,1,} iter={19} input={-13.1973,} output={1.67891,} objective={1.67891}
rank={0,1,} iter={20} input={-9.9718,} output={-1.84e-06,} objective={-1.84e-06}
rank={0,1,} iter={21} input={-13.1397,} output={1.84537,} objective={1.84537}
rank={0,1,} iter={22} input={0.420034,} output={-0.0275595,} objective={-0.0275595}
rank={0,1,} iter={23} input={-1.95687,} output={-0.00302584,} objective={-0.00302584}
rank={0,1,} iter={24} input={-17.6735,} output={1.38102,} objective={1.38102}
final_iter={25} inputs={-13.1186,} output={1.84762,}
rank={0,1,} iter={0} input={-13.0108,} output={-0.00417124,} objective={-0.00417124}
rank={0,1,} iter={1} input={-20.81,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={2} input={-5.33078,} output={-0.0042568,} objective={-0.0042568}
rank={0,1,} iter={3} input={-21.4439,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={4} input={-0.386922,} output={-0.00855468,} objective={-0.00855468}
rank={0,1,} iter={5} input={-25.6586,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={6} input={-1.30275,} output={-0.00655612,} objective={-0.00655612}
rank={0,1,} iter={7} input={-0.727284,} output={-0.00772588,} objective={-0.00772588}
rank={0,1,} iter={8} input={-13.2919,} output={-0.00417124,} objective={-0.00417124}
rank={0,1,} iter={9} input={-22.9858,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={10} input={-3.90022,} output={-0.00447516,} objective={-0.00447516}
rank={0,1,} iter={11} input={-13.1126,} output={-0.00417124,} objective={-0.00417124}
rank={0,1,} iter={12} input={-23.3232,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={13} input={-3.08018,} output={-0.00477972,} objective={-0.00477972}
rank={0,1,} iter={14} input={-25.184,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={15} input={-21.0089,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={16} input={-17.0342,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={17} input={1.2125,} output={-0.0152096,} objective={-0.0152096}
rank={0,1,} iter={18} input={-12.5215,} output={-0.0041712,} objective={-0.0041712}
rank={0,1,} iter={19} input={-0.903541,} output={-0.00732604,} objective={-0.00732604}
rank={0,1,} iter={20} input={-27.3946,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={21} input={-22.9361,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={22} input={-3.80741,} output={-0.00449856,} objective={-0.00449856}
rank={0,1,} iter={23} input={-25.2687,} output={-0.00417116,} objective={-0.00417116}
rank={0,1,} iter={24} input={0.739675,} output={-0.012838,} objective={-0.012838}
final_iter={25} inputs={-13.0108,} output={-0.00417124,}
Visual comparison of the u wind isosurfaces¶
Copied!
fig = plt.figure(figsize=(25, 18))
gs = gridspec.GridSpec(
3, 5, left=0.035, right=0.925, top=0.915, bottom=0.025, wspace=0.3, hspace=0.2
)
plot_isosurface(
uf48,
uf48,
1.0,
fig.add_subplot(gs[0, 0], projection="3d", computed_zorder=False),
"Original",
"u",
r"U component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=0,
span=np.abs(level_u) * 15,
n=16,
)
plot_isosurface(
uf48_zfp,
uf48,
uf48_zfp_cr,
fig.add_subplot(gs[0, 1], projection="3d", computed_zorder=False),
r"ZFP($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=0.5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sz3,
uf48,
uf48_sz3_cr,
fig.add_subplot(gs[0, 2], projection="3d", computed_zorder=False),
r"SZ3($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sperr,
uf48,
uf48_sperr_cr,
fig.add_subplot(gs[0, 3], projection="3d", computed_zorder=False),
r"SPERR($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sg_iso,
uf48,
uf48_sg_iso_cr,
fig.add_subplot(gs[0, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_sg["zero"],
uf48,
uf48_sg_cr["zero"],
fig.add_subplot(gs[1, 0], projection="3d", computed_zorder=False),
r"Safeguarded(0, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_sg["zfp.rs"],
uf48,
uf48_sg_cr["zfp.rs"],
fig.add_subplot(gs[1, 1], projection="3d", computed_zorder=False),
r"Safeguarded(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zfp,
error=True,
)
plot_isosurface(
uf48_sg["sz3.rs"],
uf48,
uf48_sg_cr["sz3.rs"],
fig.add_subplot(gs[1, 2], projection="3d", computed_zorder=False),
r"Safeguarded(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_sz3,
error=True,
)
plot_isosurface(
uf48_sg["sperr.rs"],
uf48,
uf48_sg_cr["sperr.rs"],
fig.add_subplot(gs[1, 3], projection="3d", computed_zorder=False),
r"Safeguarded(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_sperr,
error=True,
)
plot_isosurface(
uf48_sg_iso_iso,
uf48,
uf48_sg_iso_iso_cr,
fig.add_subplot(gs[1, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso*)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_optzconfig["zfp.rs"],
uf48,
uf48_optzconfig_cr["zfp.rs"],
fig.add_subplot(gs[2, 1], projection="3d", computed_zorder=False),
r"OptZConfig(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_optzconfig["sz3.rs"],
uf48,
uf48_optzconfig_cr["sz3.rs"],
fig.add_subplot(gs[2, 2], projection="3d", computed_zorder=False),
r"OptZConfig(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_optzconfig["sperr.rs"],
uf48,
uf48_optzconfig_cr["sperr.rs"],
fig.add_subplot(gs[2, 3], projection="3d", computed_zorder=False),
r"OptZConfig(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
# plt.tight_layout()
plt.savefig(Path("plots") / "isosurface-u.pdf", dpi=300)
plt.show()
fig = plt.figure(figsize=(25, 18))
gs = gridspec.GridSpec(
3, 5, left=0.035, right=0.925, top=0.915, bottom=0.025, wspace=0.3, hspace=0.2
)
plot_isosurface(
uf48,
uf48,
1.0,
fig.add_subplot(gs[0, 0], projection="3d", computed_zorder=False),
"Original",
"u",
r"U component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=0,
span=np.abs(level_u) * 15,
n=16,
)
plot_isosurface(
uf48_zfp,
uf48,
uf48_zfp_cr,
fig.add_subplot(gs[0, 1], projection="3d", computed_zorder=False),
r"ZFP($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=0.5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sz3,
uf48,
uf48_sz3_cr,
fig.add_subplot(gs[0, 2], projection="3d", computed_zorder=False),
r"SZ3($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sperr,
uf48,
uf48_sperr_cr,
fig.add_subplot(gs[0, 3], projection="3d", computed_zorder=False),
r"SPERR($\epsilon_{abs}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_sg_iso,
uf48,
uf48_sg_iso_cr,
fig.add_subplot(gs[0, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_sg["zero"],
uf48,
uf48_sg_cr["zero"],
fig.add_subplot(gs[1, 0], projection="3d", computed_zorder=False),
r"Safeguarded(0, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_sg["zfp.rs"],
uf48,
uf48_sg_cr["zfp.rs"],
fig.add_subplot(gs[1, 1], projection="3d", computed_zorder=False),
r"Safeguarded(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zfp,
error=True,
)
plot_isosurface(
uf48_sg["sz3.rs"],
uf48,
uf48_sg_cr["sz3.rs"],
fig.add_subplot(gs[1, 2], projection="3d", computed_zorder=False),
r"Safeguarded(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_sz3,
error=True,
)
plot_isosurface(
uf48_sg["sperr.rs"],
uf48,
uf48_sg_cr["sperr.rs"],
fig.add_subplot(gs[1, 3], projection="3d", computed_zorder=False),
r"Safeguarded(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_sperr,
error=True,
)
plot_isosurface(
uf48_sg_iso_iso,
uf48,
uf48_sg_iso_iso_cr,
fig.add_subplot(gs[1, 4], projection="3d", computed_zorder=False),
r"Safeguarded(0, iso*)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=5,
span=np.abs(level_u) * 15,
n=16,
corr=uf48_zero,
error=True,
)
plot_isosurface(
uf48_optzconfig["zfp.rs"],
uf48,
uf48_optzconfig_cr["zfp.rs"],
fig.add_subplot(gs[2, 1], projection="3d", computed_zorder=False),
r"OptZConfig(ZFP, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_optzconfig["sz3.rs"],
uf48,
uf48_optzconfig_cr["sz3.rs"],
fig.add_subplot(gs[2, 2], projection="3d", computed_zorder=False),
r"OptZConfig(SZ3, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
plot_isosurface(
uf48_optzconfig["sperr.rs"],
uf48,
uf48_optzconfig_cr["sperr.rs"],
fig.add_subplot(gs[2, 3], projection="3d", computed_zorder=False),
r"OptZConfig(SPERR, $\epsilon_{abs} \cup \text{iso}$)",
"u",
r"Absolute error over u component of wind ($m \cdot s^{-1}$)",
level=level_u,
my_span=eb_abs_u,
span=np.abs(level_u) * 15,
n=16,
error=True,
)
# plt.tight_layout()
plt.savefig(Path("plots") / "isosurface-u.pdf", dpi=300)
plt.show()
Copied!
iso_u_table = pd.concat(
[
table_isosurface(
uf48_sg_lossless["zero"],
uf48,
uf48_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg["zero"],
uf48,
uf48_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg_iso,
uf48,
uf48_sg_iso_cr,
["0", r"$\text{iso}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg_iso_iso,
uf48,
uf48_sg_iso_iso_cr,
["0", r"$\text{iso*}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_zfp,
uf48,
uf48_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["zfp.rs"],
uf48,
uf48_sg_lossless_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_zfp,
),
table_isosurface(
uf48_sg["zfp.rs"],
uf48,
uf48_sg_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_zfp,
),
table_isosurface(
uf48_optzconfig["zfp.rs"],
uf48,
uf48_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sz3,
uf48,
uf48_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["sz3.rs"],
uf48,
uf48_sg_lossless_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_sz3,
),
table_isosurface(
uf48_sg["sz3.rs"],
uf48,
uf48_sg_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_sz3,
),
table_isosurface(
uf48_optzconfig["sz3.rs"],
uf48,
uf48_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sperr,
uf48,
uf48_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["sperr.rs"],
uf48,
uf48_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"lossless",
],
"u",
level_u,
uf48_sperr,
),
table_isosurface(
uf48_sg["sperr.rs"],
uf48,
uf48_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"one-shot",
],
"u",
level_u,
uf48_sperr,
),
table_isosurface(
uf48_optzconfig["sperr.rs"],
uf48,
uf48_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_zstd,
uf48,
uf48_zstd_cr,
["ZSTD(20)", "-", ""],
"u",
level_u,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections"])
Path("tables").joinpath("isosurface-u.tex").write_text(
iso_u_table.to_latex(escape=False)
.replace("%", r"\%")
.replace("\\cline{1-11} \\cline{2-11}\n\\bottomrule", "\\bottomrule")
)
iso_u_table
iso_u_table = pd.concat(
[
table_isosurface(
uf48_sg_lossless["zero"],
uf48,
uf48_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg["zero"],
uf48,
uf48_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg_iso,
uf48,
uf48_sg_iso_cr,
["0", r"$\text{iso}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_sg_iso_iso,
uf48,
uf48_sg_iso_iso_cr,
["0", r"$\text{iso*}$", "one-shot"],
"u",
level_u,
uf48_zero,
),
table_isosurface(
uf48_zfp,
uf48,
uf48_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["zfp.rs"],
uf48,
uf48_sg_lossless_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_zfp,
),
table_isosurface(
uf48_sg["zfp.rs"],
uf48,
uf48_sg_cr["zfp.rs"],
[r"ZFP($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_zfp,
),
table_isosurface(
uf48_optzconfig["zfp.rs"],
uf48,
uf48_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sz3,
uf48,
uf48_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["sz3.rs"],
uf48,
uf48_sg_lossless_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "lossless"],
"u",
level_u,
uf48_sz3,
),
table_isosurface(
uf48_sg["sz3.rs"],
uf48,
uf48_sg_cr["sz3.rs"],
[r"SZ3($\epsilon_{abs}$)", r"$\epsilon_{abs} \cup \text{iso}$", "one-shot"],
"u",
level_u,
uf48_sz3,
),
table_isosurface(
uf48_optzconfig["sz3.rs"],
uf48,
uf48_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sperr,
uf48,
uf48_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_sg_lossless["sperr.rs"],
uf48,
uf48_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"lossless",
],
"u",
level_u,
uf48_sperr,
),
table_isosurface(
uf48_sg["sperr.rs"],
uf48,
uf48_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{iso}$",
"one-shot",
],
"u",
level_u,
uf48_sperr,
),
table_isosurface(
uf48_optzconfig["sperr.rs"],
uf48,
uf48_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{iso}$", ""],
"u",
level_u,
None,
),
table_isosurface(
uf48_zstd,
uf48,
uf48_zstd_cr,
["ZSTD(20)", "-", ""],
"u",
level_u,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections"])
Path("tables").joinpath("isosurface-u.tex").write_text(
iso_u_table.to_latex(escape=False)
.replace("%", r"\%")
.replace("\\cline{1-11} \\cline{2-11}\n\\bottomrule", "\\bottomrule")
)
iso_u_table
| $L_{\infty}(\hat{u})$ | $L_{2}(\hat{u})$ | FN | FP | FS | V | C | CR | |||
|---|---|---|---|---|---|---|---|---|---|---|
| Compressor | Safeguarded | Corrections | ||||||||
| 0 | $\epsilon_{abs} \cup \text{iso}$ | lossless | 5.0 | 1.263 | 0 | 0 | 0 | 0 | 70.4% | $\times$ 2.06 |
| one-shot | 5.0 | 2.622 | 0 | 0 | 0 | 0 | 70.4% | $\times$ 105.86 | ||
| $\text{iso}$ | one-shot | 3.689e+19 | inf [inf] | 0 | 0 | 0 | 0 | 50.8% | $\times$ 338.1 | |
| $\text{iso*}$ | one-shot | 3.689e+19 | inf [inf] | 0 | 0 | 0 | 0 | 79.7% | $\times$ 92.39 | |
| ZFP($\epsilon_{abs}$) | - | NaN [1.761] | NaN [0.1323] | 6.8% | 6.8% | 55.6% | 1.0% | $\times$ 44.13 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 1.761 | 0.131 | 0 | 0 | 0 | 0 | 1.0% | $\times$ 32.17 | |
| one-shot | 2.232 | 0.1833 | 0 | 0 | 0 | 0 | 1.0% | $\times$ 38.78 | ||
| OptZConfig(ZFP) | $\epsilon_{abs} \cup \text{iso}$ | NaN [1.907e-06] | NaN [1.386e-07] | 0 | 0 | 0 | 0.4% | $\times$ 1.7 | ||
| SZ3($\epsilon_{abs}$) | - | 4.999 | 1.042 | 66.9% | 46.1% | 78.3% | 4.7% | $\times$ 952.92 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 4.999 | 0.9953 | 0 | 0 | 0 | 0 | 4.7% | $\times$ 20.7 | |
| one-shot | 5.0 | 1.036 | 0 | 0 | 0 | 0 | 4.7% | $\times$ 150.68 | ||
| OptZConfig(SZ3) | $\epsilon_{abs} \cup \text{iso}$ | 2.007e-06 | 1.094e-06 | 0 | 0 | 0 | 0 | $\times$ 1.85 | ||
| SPERR($\epsilon_{abs}$) | - | NaN [4.992] | NaN [0.3102] | 17.8% | 21.4% | 75.0% | 1.8% | $\times$ 1446.38 | ||
| $\epsilon_{abs} \cup \text{iso}$ | lossless | 4.992 | 0.3041 | 0 | 0 | 0 | 0 | 1.8% | $\times$ 50.22 | |
| one-shot | 4.992 | 0.3516 | 0 | 0 | 0 | 0 | 1.8% | $\times$ 115.57 | ||
| OptZConfig(SPERR) | $\epsilon_{abs} \cup \text{iso}$ | NaN [3.815e-06] | NaN [9.336e-07] | 0 | 10.4% | 2.9% | 0.4% | $\times$ 2.34 | ||
| ZSTD(20) | - | 0.0 | 0.0 | 0 | 0 | 0 | 0 | $\times$ 1.12 |
Copied!
import json
with Path("observations").joinpath("isosurface.json").open("w") as f:
json.dump(observations, f)
import json
with Path("observations").joinpath("isosurface.json").open("w") as f:
json.dump(observations, f)
Preserving multiple isosurfaces¶
As shown above, the safeguards can also preserve an arbitrary number of isosurfaces by using multiple sign safeguards with offsets that match the isosurface values:
SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="sign", offset=-10),
dict(kind="sign", offset=0),
dict(kind="sign", offset=10),
...
],
)
These offsets can also be late-bound, e.g. offset="offset", to use a pointwise-varying isosurface values.
Copied!