9. Uncertainty analysis

Here we will demonstrate a full uncertainty analysis of the inversion. We use a stochastic approach, where we 1) choose the important input parameters to the inversion, 2) define each of there uncertainty distributions, 3) run a series of inversions which sample these inputs from their uncertainty distributions, and 4) use the ensemble of inverted topography models to define the mean result and the uncertainty.

9.1. Import packages

%load_ext autoreload
%autoreload 2
import copy
import logging
import pathlib
import pickle
import string

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import verde as vd
import xarray as xr
from polartoolkit import maps, profiles
from polartoolkit import utils as polar_utils

from invert4geom import (
    plotting,
    regional,
    uncertainty,
    utils,
)

# set up logging to see what's going on
logging.basicConfig(level=logging.INFO)

/home/mdtanker/miniforge3/envs/invert4geom/lib/python3.12/site-packages/UQpy/__init__.py:6: UserWarning:

pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.

10. Load results from notebook Combining it all

# load pickle files
params = []
with pathlib.Path("../tmp/combining_it_all_results.pickle").open("rb") as file:
    while True:
        try:
            params.append(pickle.load(file))
        except EOFError:
            break
(
    grav_df,
    constraint_points,
    starting_topography_kwargs,
    grids,
    inversion_results,
    kwargs,
    regional_grav_kwargs,
    _,
) = params

grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()
grav_df.describe()

	northing	easting	upward	upper_surface_grav	lower_surface_grav	observed_grav	gravity_anomaly	uncert	starting_gravity	true_res	misfit	reg	res
count	4941.000000	4941.000000	4941.0	4941.000000	4941.000000	4941.000000	4941.000000	4.941000e+03	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000
mean	15000.000000	20000.000000	1001.0	-0.229678	0.004286	-0.225392	-0.225332	1.000000e-01	1.025253	-1.254871	-1.250585	-0.826227	-0.424358
std	8804.299419	11691.635129	0.0	7.103507	0.720692	7.247722	7.245591	1.064189e-10	6.860052	4.621594	4.852437	1.372783	4.483083
min	0.000000	0.000000	1001.0	-17.227928	-1.034071	-16.940400	-17.040892	1.000000e-01	-8.040797	-19.592049	-19.972713	-3.252807	-18.401653
25%	7500.000000	10000.000000	1001.0	-5.587227	-0.610336	-5.906763	-5.910496	1.000000e-01	-4.472542	-3.008015	-3.123793	-2.048533	-2.142099
50%	15000.000000	20000.000000	1001.0	-1.239580	-0.155380	-0.822890	-0.813296	1.000000e-01	-0.749905	-1.483188	-1.843795	-0.733369	-0.346187
75%	22500.000000	30000.000000	1001.0	3.490030	0.542849	3.738912	3.726649	1.000000e-01	5.319529	0.577570	0.952694	0.418303	1.567960
max	30000.000000	40000.000000	1001.0	19.754954	1.602413	19.932831	19.959187	1.000000e-01	23.360744	13.969991	14.301613	1.238055	14.384864

constraint_points.describe()

	easting	northing	true_upward	upward	uncert	weight	inverted_topography
count	15.000000	15.000000	15.000000	15.000000	15.000000	15.000000	15.000000
mean	18971.649226	16518.655913	489.096872	487.025168	9.781937	0.011132	490.371247
std	10839.669022	8540.046257	79.152932	81.118110	1.583059	0.003039	85.481722
min	3505.581182	1697.178371	406.379413	405.305719	8.127588	0.005624	396.314242
25%	11575.107009	10553.975772	429.926990	424.663596	8.598540	0.009240	429.575606
50%	19995.535031	16357.717739	454.924262	458.214230	9.098485	0.012080	455.379898
75%	27651.208305	24233.535276	520.183444	514.207570	10.403669	0.013525	527.969934
max	38163.601456	27603.626012	666.725645	668.402196	13.334513	0.015138	673.102522

true_topography = grids.true
lower_topography = grids.lower
starting_topography = grids.starting
final_topography = grids.inverted

topo_results, grav_results, parameters, elapsed_time = inversion_results

kwargs

{'max_iterations': 200,
 'l2_norm_tolerance': 0.3,
 'delta_l2_norm_tolerance': 1.008}

regional_grav_kwargs["constraints_df"] = constraint_points
regional_grav_kwargs

{'method': 'constraints',
 'grid_method': 'eq_sources',
 'constraints_weights_column': 'weight',
 'cv': True,
 'cv_kwargs': {'n_trials': 50,
  'damping_limits': (1e-10, 1),
  'progressbar': False,
  'fname': '../tmp/regional_sep'},
 'depth': 'default',
 'block_size': None,
 'constraints_df':          easting      northing  true_upward      upward     uncert    weight  \
 0    3899.714996  26468.618776   666.725645  668.402196  13.334513  0.005624
 1   30636.914105   6974.789900   406.379413  405.305719   8.127588  0.015138
 2   17659.550795  13659.470931   480.393921  486.547024   9.607878  0.010833
 3   28491.676758  27073.768551   517.183455  518.268507  10.343669  0.009347
 4   38163.601456   1697.178371   435.459979  430.794727   8.709200  0.013184
 5   21462.843076  17815.369689   454.924262  458.214230   9.098485  0.012080
 6   20042.577619  27603.626012   411.155008  421.877931   8.223100  0.014789
 7    3737.943068   7448.480613   582.902996  593.944123  11.658060  0.007358
 8   11200.681244  16357.717739   511.624141  504.423182  10.232483  0.009551
 9   19995.535031  26455.594497   450.505648  439.104058   9.010113  0.012318
 10  26810.739853   4728.744481   430.317216  424.953101   8.606344  0.013501
 11  31542.083372  15655.552259   424.023626  424.374090   8.480473  0.013905
 12  15475.763060  22011.476055   429.536764  409.563040   8.590735  0.013550
 13   3505.581182  19732.370745   612.137575  609.458963  12.242752  0.006672
 14  11949.532774  14097.080073   523.183432  510.146633  10.463669  0.009133

     inverted_topography
 0            673.102522
 1            396.314242
 2            484.695808
 3            511.140548
 4            439.034625
 5            455.379898
 6            399.686372
 7            602.485537
 8            522.819184
 9            443.471207
 10           433.165380
 11           415.195606
 12           425.985833
 13           619.971257
 14           533.120685  }

starting_topography_kwargs

{'method': 'splines',
 'region': (0, 40000, 0, 30000),
 'spacing': 1000,
 'constraints_df':          easting      northing  true_upward      upward     uncert    weight
 0    3899.714996  26468.618776   666.725645  668.402196  13.334513  0.005624
 1   30636.914105   6974.789900   406.379413  405.305719   8.127588  0.015138
 2   17659.550795  13659.470931   480.393921  486.547024   9.607878  0.010833
 3   28491.676758  27073.768551   517.183455  518.268507  10.343669  0.009347
 4   38163.601456   1697.178371   435.459979  430.794727   8.709200  0.013184
 5   21462.843076  17815.369689   454.924262  458.214230   9.098485  0.012080
 6   20042.577619  27603.626012   411.155008  421.877931   8.223100  0.014789
 7    3737.943068   7448.480613   582.902996  593.944123  11.658060  0.007358
 8   11200.681244  16357.717739   511.624141  504.423182  10.232483  0.009551
 9   19995.535031  26455.594497   450.505648  439.104058   9.010113  0.012318
 10  26810.739853   4728.744481   430.317216  424.953101   8.606344  0.013501
 11  31542.083372  15655.552259   424.023626  424.374090   8.480473  0.013905
 12  15475.763060  22011.476055   429.536764  409.563040   8.590735  0.013550
 13   3505.581182  19732.370745   612.137575  609.458963  12.242752  0.006672
 14  11949.532774  14097.080073   523.183432  510.146633  10.463669  0.009133,
 'dampings': array([1.00000000e-40, 2.53536449e-40, 6.42807312e-40, 1.62975083e-39,
        4.13201240e-39, 1.04761575e-38, 2.65608778e-38, 6.73415066e-38,
        1.70735265e-37, 4.32876128e-37, 1.09749877e-36, 2.78255940e-36,
        7.05480231e-36, 1.78864953e-35, 4.53487851e-35, 1.14975700e-34,
        2.91505306e-34, 7.39072203e-34, 1.87381742e-33, 4.75081016e-33,
        1.20450354e-32, 3.05385551e-32, 7.74263683e-32, 1.96304065e-31,
        4.97702356e-31, 1.26185688e-30, 3.19926714e-30, 8.11130831e-30,
        2.05651231e-29, 5.21400829e-29, 1.32194115e-28, 3.35160265e-28,
        8.49753436e-28, 2.15443469e-27, 5.46227722e-27, 1.38488637e-26,
        3.51119173e-26, 8.90215085e-26, 2.25701972e-25, 5.72236766e-25,
        1.45082878e-24, 3.67837977e-24, 9.32603347e-24, 2.36448941e-23,
        5.99484250e-23, 1.51991108e-22, 3.85352859e-22, 9.77009957e-22,
        2.47707636e-21, 6.28029144e-21, 1.59228279e-20, 4.03701726e-20,
        1.02353102e-19, 2.59502421e-19, 6.57933225e-19, 1.66810054e-18,
        4.22924287e-18, 1.07226722e-17, 2.71858824e-17, 6.89261210e-17,
        1.74752840e-16, 4.43062146e-16, 1.12332403e-15, 2.84803587e-15,
        7.22080902e-15, 1.83073828e-14, 4.64158883e-14, 1.17681195e-13,
        2.98364724e-13, 7.56463328e-13, 1.91791026e-12, 4.86260158e-12,
        1.23284674e-11, 3.12571585e-11, 7.92482898e-11, 2.00923300e-10,
        5.09413801e-10, 1.29154967e-09, 3.27454916e-09, 8.30217568e-09,
        2.10490414e-08, 5.33669923e-08, 1.35304777e-07, 3.43046929e-07,
        8.69749003e-07, 2.20513074e-06, 5.59081018e-06, 1.41747416e-05,
        3.59381366e-05, 9.11162756e-05, 2.31012970e-04, 5.85702082e-04,
        1.48496826e-03, 3.76493581e-03, 9.54548457e-03, 2.42012826e-02,
        6.13590727e-02, 1.55567614e-01, 3.94420606e-01, 1.00000000e+00]),
 'weights_col': 'weight'}

11. Extract parameters

parameters

{'Density contrast(s)': '[2731] kg/m3',
 'Reference level': '492.7997381853817 m',
 'Max iterations': 200,
 'L2 norm tolerance': '0.3',
 'Delta L2 norm tolerance': '1.008',
 'Deriv type': 'annulus',
 'Solver type': 'scipy least squares',
 'Solver damping': 0.02239318117126179,
 'Upper confining layer': 'Not enabled',
 'Lower confining layer': 'Not enabled',
 'Regularization weighting grid': 'Not enabled',
 'Time elapsed': '3 seconds',
 'Avg. iteration time': '1.15 seconds',
 'Final misfit RMSE / L2-norm': '0.0697 /0.2641 mGal',
 'Termination reason': ['l2-norm tolerance'],
 'Iteration times': [1.2237195130001055,
  1.4033744970001862,
  0.8336053599996376]}

density_contrast = float(parameters["Density contrast(s)"][1:-7])
zref = float(parameters["Reference level"][:-2])
spacing = utils.get_spacing(grav_df)
region = vd.get_region((grav_df.easting, grav_df.northing))
solver_damping = parameters["Solver damping"]

print(f"Density contrast: {density_contrast} kg/m3")
print(f"Reference level: {zref} m")
print(f"Spacing: {spacing} m")
print(f"Region: {region}")
print(f"Solver damping: {solver_damping}")

Density contrast: 2731.0 kg/m3
Reference level: 492.7997381853817 m
Spacing: 500.0 m
Region: (np.float64(0.0), np.float64(40000.0), np.float64(0.0), np.float64(30000.0))
Solver damping: 0.02239318117126179

# collect the results
topo_results, grav_results, parameters, elapsed_time = inversion_results

plotting.plot_inversion_results(
    grav_results,
    topo_results,
    parameters,
    region,
    iters_to_plot=2,
    plot_iter_results=False,
    plot_topo_results=True,
    plot_grav_results=True,
)

final_topography = topo_results.set_index(["northing", "easting"]).to_xarray().topo

_ = polar_utils.grd_compare(
    true_topography,
    final_topography,
    grid1_name="True topography",
    grid2_name="Inverted topography",
    robust=True,
    hist=True,
    inset=False,
    title="difference",
    grounding_line=False,
    reverse_cpt=True,
    cmap="rain",
    points=constraint_points,
    points_style="x.3c",
)

12. Absolute value of inversion error

inversion_error = np.abs(true_topography - final_topography)
fig = maps.plot_grd(
    inversion_error,
    hist=True,
    cmap="thermal",
    title="Absolute value of inversion error",
    robust=True,
    points=constraint_points,
    points_style="x.4c,",
    points_pen="3p",
    points_fill="white",
)
fig.show()

plot [WARNING]: .4c, not a valid number and may not be decoded properly.

13. Redo starting gravity and regional separation

use the optimal determined density contrast and zref values

# set the reference level from the prisms to 0
density_grid = xr.where(
    starting_topography >= zref,
    density_contrast,
    -density_contrast,
)
# create layer of prisms
starting_prisms = utils.grids_to_prisms(
    starting_topography,
    zref,
    density=density_grid,
)
grav_df["starting_gravity"] = starting_prisms.prism_layer.gravity(
    coordinates=(
        grav_df.easting,
        grav_df.northing,
        grav_df.upward,
    ),
    field="g_z",
    progressbar=False,
)

# regional separation
temp_reg_kwargs = copy.deepcopy(regional_grav_kwargs)

# temporarily set some kwargs
# temp_reg_kwargs["cv_kwargs"]["plot"]=True
temp_reg_kwargs["cv_kwargs"]["progressbar"] = True

grav_df = regional.regional_separation(
    grav_df=grav_df,
    **temp_reg_kwargs,
)

grav_df["staring_res"] = grav_df.res
grav_df["starting_reg"] = grav_df.reg

grav_df.describe()

INFO:invert4geom:using 12 startup trials

INFO:invert4geom:Trial with best score:
INFO:invert4geom:       trial number: 45
INFO:invert4geom:       parameter: {'damping': 0.049316247428762484}
INFO:invert4geom:       scores: [-0.3617634753499722]

	northing	easting	upward	upper_surface_grav	lower_surface_grav	observed_grav	gravity_anomaly	uncert	starting_gravity	true_res	misfit	reg	res	staring_res	starting_reg
count	4941.000000	4941.000000	4941.0	4941.000000	4941.000000	4941.000000	4941.000000	4.941000e+03	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000
mean	15000.000000	20000.000000	1001.0	-0.229678	0.004286	-0.225392	-0.225332	1.000000e-01	0.516151	-1.254871	-0.741482	-0.135895	-0.605588	-0.605588	-0.135895
std	8804.299419	11691.635129	0.0	7.103507	0.720692	7.247722	7.245591	1.064189e-10	7.527778	4.621594	5.072202	1.368514	4.571823	4.571823	1.368514
min	0.000000	0.000000	1001.0	-17.227928	-1.034071	-16.940400	-17.040892	1.000000e-01	-9.417770	-19.592049	-20.044860	-2.543056	-18.520552	-18.520552	-2.543056
25%	7500.000000	10000.000000	1001.0	-5.587227	-0.610336	-5.906763	-5.910496	1.000000e-01	-5.492789	-3.008015	-2.796283	-1.225539	-2.336660	-2.336660	-1.225539
50%	15000.000000	20000.000000	1001.0	-1.239580	-0.155380	-0.822890	-0.813296	1.000000e-01	-1.441579	-1.483188	-1.207644	-0.193214	-0.616694	-0.616694	-0.193214
75%	22500.000000	30000.000000	1001.0	3.490030	0.542849	3.738912	3.726649	1.000000e-01	5.240259	0.577570	1.435570	1.018149	1.262640	1.262640	1.018149
max	30000.000000	40000.000000	1001.0	19.754954	1.602413	19.932831	19.959187	1.000000e-01	25.010049	13.969991	15.345193	2.269982	14.670183	14.670183	2.269982

# plotting functions for uncertainty results
def uncert_plots(
    results,
    inversion_region,
    spacing,
    true_topography,
    constraint_points=None,
    weight_by=None,
):
    if (weight_by == "constraints") & (constraint_points is None):
        msg = "must provide constraint_points if weighting by constraints"
        raise ValueError(msg)

    stats_ds = uncertainty.merged_stats(
        results=results,
        plot=True,
        constraints_df=constraint_points,
        weight_by=weight_by,
        region=inversion_region,
    )

    try:
        mean = stats_ds.weighted_mean
        stdev = stats_ds.weighted_stdev
    except AttributeError:
        mean = stats_ds.z_mean
        stdev = stats_ds.z_stdev

    _ = polar_utils.grd_compare(
        true_topography,
        mean,
        region=vd.pad_region(inversion_region, -spacing),
        grid1_name="True topography",
        grid2_name="Inverted topography",
        robust=True,
        hist=True,
        inset=False,
        title="difference",
        grounding_line=False,
        reverse_cpt=True,
        cmap="rain",
        points=constraint_points,
        points_style="x.3c",
    )
    _ = polar_utils.grd_compare(
        np.abs(true_topography - mean),
        stdev,
        region=vd.pad_region(inversion_region, -spacing),
        grid1_name="True error",
        grid2_name="Stochastic uncertainty",
        cmap="thermal",
        robust=True,
        hist=True,
        inset=False,
        title="difference",
        grounding_line=False,
        points=constraint_points,
        points_style="x.3c",
        points_fill="white",
    )
    return stats_ds

14. Uncertainty analysis

# get mean spacing between nearest constraints
mean_spacing = np.mean(
    vd.median_distance(
        (constraint_points.easting, constraint_points.northing),
        k_nearest=1,
    )
)
mean_spacing

np.float64(5599.156508868462)

# use optimal eq source parameters for regional separation in uncertainty analysis
# re-load the study from the saved pickle file
with pathlib.Path(f"{regional_grav_kwargs.get('cv_kwargs').get('fname')}.pickle").open(
    "rb"
) as f:
    study = pickle.load(f)
reg_eq_damping = min(study.best_trials, key=lambda t: t.values[0]).params["damping"]

reg_eq_depth = 4.5 * mean_spacing

new_regional_grav_kwargs = dict(
    method="constraints",
    grid_method="eq_sources",
    constraints_df=constraint_points,
    damping=reg_eq_damping,
    depth=reg_eq_depth,
    block_size=None,
)
reg_eq_damping, reg_eq_depth

(0.049316247428762484, np.float64(25196.204289908077))

# add best damping parameter to kwargs
kwargs["solver_damping"] = solver_damping

# kwargs to reuse for all uncertainty analyses
uncert_kwargs = dict(
    grav_df=grav_df,
    density_contrast=density_contrast,
    zref=zref,
    starting_prisms=starting_prisms,
    starting_topography=starting_topography,
    regional_grav_kwargs=new_regional_grav_kwargs,
    **kwargs,
)

14.1. Solver damping component

To estimate the uncertainty of the choice of solver damping value, we will load the damping cross validation and calculate the standard deviation of the trials (excluding outliers)

# load study
with pathlib.Path("../tmp/combining_it_all_damping_CV_study.pickle").open("rb") as f:
    study = pickle.load(f)

study_df = study.trials_dataframe().drop(columns=["user_attrs_results"])
study_df.sort_values("value", inplace=True)

# calculate zscores of values
study_df["value_zscore"] = sp.stats.zscore(study_df["value"])

# drop outliers (values with zscore > |2|)
study_df2 = study_df[(np.abs(study_df.value_zscore) < 2)]

# pick damping standard deviation based on optimization
stdev = np.log10(study_df2.params_damping).std()
print(f"calculated stdev: {stdev}")
# stdev = .1
print(f"using stdev: {stdev}")

calculated stdev: 0.6606916830690623
using stdev: 0.6606916830690623

fig = plotting.plot_cv_scores(
    study_df.value,
    study_df.params_damping,
    param_name="Damping",
    logx=True,
    logy=True,
)
ax = fig.axes[0]

best = float(study_df2.params_damping.iloc[0])
upper = float(10 ** (np.log10(best) + stdev))
lower = float(10 ** (np.log10(best) - stdev))

y_lims = ax.get_ylim()
ax.vlines(best, ymin=y_lims[0], ymax=y_lims[1], color="r")
ax.vlines(upper, ymin=y_lims[0], ymax=y_lims[1], label="+/- std")
ax.vlines(lower, ymin=y_lims[0], ymax=y_lims[1])

x_lims = ax.get_xlim()
ax.set_xlim(
    min(x_lims[0], lower),
    max(x_lims[1], upper),
)
ax.legend()
print("best:", best, "\nstd:", stdev, "\n+1std:", upper, "\n-1std:", lower)

best: 0.02239318117126179
std: 0.6606916830690623
+1std: 0.10251973545005608
-1std: 0.004891297863456191

solver_dict = {
    "solver_damping": {
        "distribution": "normal",
        "loc": np.log10(solver_damping),  # mean of base 10 exponent
        "scale": stdev,  # standard deviation of base 10 exponent
        "log": True,
    },
}
fname = "../tmp/uncertainty_damping"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_damping_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=10,
    parameter_dict=solver_dict,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_damping_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:Sampled 'solver_damping' parameter values; mean: 0.055707817615901226, min: 0.0018338959126522312, max: 0.2734367635095143
INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_damping' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_damping'

14.2. Z-ref and density components

To estimate the uncertainty of the choice of Z-ref and density value, we will sample the values from random distributions.

# choose standard deviations for the density contrast and zref
density_stdev = 100
zref_stdev = 50

# to re-load the study from the saved pickle file
with pathlib.Path("../tmp/combining_it_all_density_and_zref_cv_study.pickle").open(
    "rb"
) as f:
    study = pickle.load(f)

plotting.plot_2_parameter_cv_scores_uneven(
    study,
    param_names=(
        "params_zref",
        "params_density_contrast",
    ),
    plot_param_names=(
        "Reference level (m)",
        "Density contrast (kg/m$^3$)",
    ),
)

ax = plt.gca()

# plot zref stdev
upper = zref + zref_stdev
lower = zref - zref_stdev
y_lims = ax.get_ylim()
ax.vlines(upper, ymin=y_lims[0], ymax=y_lims[1], label="+/- std")
ax.vlines(lower, ymin=y_lims[0], ymax=y_lims[1])
x_lims = ax.get_xlim()
ax.set_xlim(
    min(x_lims[0], lower),
    max(x_lims[1], upper),
)
print("best:", zref, "\nstd:", zref_stdev, "\n+1std:", upper, "\n-1std:", lower)

# plot density contrast stdev
upper = density_contrast + density_stdev
lower = density_contrast - density_stdev
x_lims = ax.get_xlim()
ax.hlines(upper, xmin=x_lims[0], xmax=x_lims[1])
ax.hlines(lower, xmin=x_lims[0], xmax=x_lims[1])
y_lims = ax.get_ylim()
ax.set_ylim(
    min(y_lims[0], lower),
    max(y_lims[1], upper),
)
ax.legend()
print(
    "best:",
    density_contrast,
    "\nstd:",
    density_stdev,
    "\n+1std:",
    upper,
    "\n-1std:",
    lower,
)

best: 492.7997381853817
std: 50
+1std: 542.7997381853818
-1std: 442.7997381853817
best: 2731.0
std: 100
+1std: 2831.0
-1std: 2631.0

density_dict = {
    "density_contrast": {
        "distribution": "normal",
        "loc": density_contrast,
        "scale": density_stdev,
    },
}
fname = "../tmp/uncertainty_density"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_density_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=10,
    parameter_dict=density_dict,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_density_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:Sampled 'density_contrast' parameter values; mean: 2731.0, min: 2566.514637304853, max: 2895.485362695147
INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_density' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_density'

zref_dict = {
    "zref": {
        "distribution": "normal",
        "loc": zref,
        "scale": zref_stdev,
    },
}
fname = "../tmp/uncertainty_zref"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_zref_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=10,
    parameter_dict=zref_dict,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_zref_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:Sampled 'zref' parameter values; mean: 492.79973818538167, min: 410.5570568378081, max: 575.0424195329554
INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_zref' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_zref'

14.3. Constraints component

fname = "../tmp/uncertainty_constraints"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_constraints_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    sample_constraints=True,
    constraints_df=constraint_points,
    starting_topography_kwargs=starting_topography_kwargs,
    runs=10,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_constraints_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_constraints' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_constraints'

14.4. Gravity component

grav_df.describe()

	northing	easting	upward	upper_surface_grav	lower_surface_grav	observed_grav	gravity_anomaly	uncert	starting_gravity	true_res	misfit	reg	res	staring_res	starting_reg
count	4941.000000	4941.000000	4941.0	4941.000000	4941.000000	4941.000000	4941.000000	4.941000e+03	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000	4941.000000
mean	15000.000000	20000.000000	1001.0	-0.229678	0.004286	-0.225392	-0.225332	1.000000e-01	0.516151	-1.254871	-0.741482	-0.135895	-0.605588	-0.605588	-0.135895
std	8804.299419	11691.635129	0.0	7.103507	0.720692	7.247722	7.245591	1.064189e-10	7.527778	4.621594	5.072202	1.368514	4.571823	4.571823	1.368514
min	0.000000	0.000000	1001.0	-17.227928	-1.034071	-16.940400	-17.040892	1.000000e-01	-9.417770	-19.592049	-20.044860	-2.543056	-18.520552	-18.520552	-2.543056
25%	7500.000000	10000.000000	1001.0	-5.587227	-0.610336	-5.906763	-5.910496	1.000000e-01	-5.492789	-3.008015	-2.796283	-1.225539	-2.336660	-2.336660	-1.225539
50%	15000.000000	20000.000000	1001.0	-1.239580	-0.155380	-0.822890	-0.813296	1.000000e-01	-1.441579	-1.483188	-1.207644	-0.193214	-0.616694	-0.616694	-0.193214
75%	22500.000000	30000.000000	1001.0	3.490030	0.542849	3.738912	3.726649	1.000000e-01	5.240259	0.577570	1.435570	1.018149	1.262640	1.262640	1.018149
max	30000.000000	40000.000000	1001.0	19.754954	1.602413	19.932831	19.959187	1.000000e-01	25.010049	13.969991	15.345193	2.269982	14.670183	14.670183	2.269982

grav_df.res

0      -0.542281
1      -1.131335
2      -0.864570
3      -0.748263
4      -0.658321
          ...
4936   -5.912816
4937   -6.166574
4938   -6.271337
4939   -6.202335
4940   -5.538048
Name: res, Length: 4941, dtype: float64

fname = "../tmp/uncertainty_grav"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_grav_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=10,
    sample_gravity=True,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_grav_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_grav' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_grav'

14.5. Regional gravity component

reg_eq_depth, reg_eq_damping

(np.float64(25196.204289908077), 0.049316247428762484)

regional_misfit_parameter_dict = {
    "depth": {
        "distribution": "normal",
        "loc": reg_eq_depth,  # mean
        "scale": reg_eq_depth / 10,  # standard deviation
    },
    "damping": {
        "distribution": "normal",
        "loc": np.log10(reg_eq_damping),  # mean base 10 exponent
        # "scale": np.abs(np.log10(reg_eq_damping)),  # standard deviation of exponent
        "scale": 4,
        "log": True,
    },
}

regional_misfit_stats, _ = uncertainty.regional_misfit_uncertainty(
    runs=100,
    parameter_dict=regional_misfit_parameter_dict,
    true_regional=grav_grid.lower_surface_grav,
    # weight_by="constraints",
    # weight_by="residual",
    # weight_by=None,
    grav_df=grav_df,
    **new_regional_grav_kwargs,
)

INFO:invert4geom:Sampled 'depth' parameter values; mean: 25196.204289908077, min: 18706.0921550931, max: 31686.316424723052
INFO:invert4geom:Sampled 'damping' parameter values; mean: 10199825.379978541, min: 2.452860269028061e-12, max: 991533146.4921261

grav_grid = grav_df.set_index(["northing", "easting"]).to_xarray()
grav_grid["reg_uncert"] = regional_misfit_stats.z_stdev
grav_df_new = vd.grid_to_table(grav_grid)
grav_df_new

	northing	easting	test	upward	upper_surface_grav	lower_surface_grav	observed_grav	gravity_anomaly	uncert	starting_gravity	true_res	misfit	reg	res	staring_res	starting_reg	reg_uncert
0	0.0	0.0	False	1001.0	9.042402	-0.041223	9.001179	9.017793	0.1	12.082565	-2.349335	-3.064772	-2.522491	-0.542281	-0.542281	-2.522491	3.278707
1	0.0	500.0	True	1001.0	9.380272	-0.046312	9.333959	9.331242	0.1	12.986009	-2.531914	-3.654768	-2.523433	-1.131335	-1.131335	-2.523433	3.090681
2	0.0	1000.0	False	1001.0	9.568091	-0.052122	9.515969	9.506800	0.1	12.894242	-2.657964	-3.387443	-2.522873	-0.864570	-0.864570	-2.522873	2.896165
3	0.0	1500.0	True	1001.0	9.298147	-0.059918	9.238228	9.269420	0.1	12.538464	-2.625601	-3.269044	-2.520780	-0.748263	-0.748263	-2.520780	2.696680
4	0.0	2000.0	False	1001.0	8.893385	-0.068787	8.824597	8.892680	0.1	12.068134	-2.514374	-3.175453	-2.517133	-0.658321	-0.658321	-2.517133	2.494100
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
4936	30000.0	38000.0	False	1001.0	2.982330	0.607664	3.589994	3.623280	0.1	7.479534	-4.259751	-3.856254	2.056561	-5.912816	-5.912816	2.056561	4.726586
4937	30000.0	38500.0	True	1001.0	2.931601	0.580987	3.512587	3.472411	0.1	7.606017	-4.529030	-4.133606	2.032968	-6.166574	-6.166574	2.032968	4.780391
4938	30000.0	39000.0	False	1001.0	2.946999	0.545480	3.492479	3.401948	0.1	7.664627	-4.573146	-4.262678	2.008659	-6.271337	-6.271337	2.008659	4.829202
4939	30000.0	39500.0	True	1001.0	2.873747	0.495910	3.369657	3.335928	0.1	7.554555	-4.249337	-4.218627	1.983708	-6.202335	-6.202335	1.983708	4.873132
4940	30000.0	40000.0	False	1001.0	2.784613	0.441061	3.225674	3.282663	0.1	6.862522	-3.788240	-3.579859	1.958189	-5.538048	-5.538048	1.958189	4.912319

4941 rows × 17 columns

# update grav_df
uncert_kwargs["grav_df"] = grav_df_new

fname = "../tmp/uncertainty_regional"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_regional_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=10,
    regional_misfit_parameter_dict=regional_misfit_parameter_dict,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_regional_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:Sampled 'depth' parameter values; mean: 25196.204289908077, min: 21051.797488741522, max: 29340.611091074632
INFO:invert4geom:Sampled 'damping' parameter values; mean: 18795.883330960092, min: 1.2988993965843748e-08, max: 187242.54294446588
INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_regional' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 10
saving results to pickle files with prefix: '../tmp/uncertainty_regional'

14.6. Total uncertainty

fname = "../tmp/uncertainty_full"

# delete files if already exist
for p in pathlib.Path().glob(f"{fname}*"):
    p.unlink(missing_ok=True)

uncert_results = uncertainty.full_workflow_uncertainty_loop(
    fname=fname,
    runs=20,
    sample_gravity=True,
    sample_constraints=True,
    constraints_df=constraint_points,
    parameter_dict=solver_dict | density_dict | zref_dict,
    regional_misfit_parameter_dict=regional_misfit_parameter_dict,
    starting_topography_kwargs=starting_topography_kwargs,
    **uncert_kwargs,
)

stats_ds = uncert_plots(
    uncert_results,
    region,
    spacing,
    true_topography,
    constraint_points=constraint_points,
    weight_by="constraints",
)

INFO:invert4geom:Sampled 'solver_damping' parameter values; mean: 0.060818624964864565, min: 0.001135491506553461, max: 0.44161894657496065
INFO:invert4geom:Sampled 'density_contrast' parameter values; mean: 2730.9999999999995, min: 2535.0036015459946, max: 2926.9963984540054
INFO:invert4geom:Sampled 'zref' parameter values; mean: 492.7997381853817, min: 394.801538958379, max: 590.7979374123845
INFO:invert4geom:Sampled 'depth' parameter values; mean: 25196.20428990808, min: 20257.83899437473, max: 30134.569585441426
INFO:invert4geom:Sampled 'damping' parameter values; mean: 172064.18326446132, min: 7.130731513548589e-10, max: 3410719.1609078436
INFO:invert4geom:No pickle files starting with '../tmp/uncertainty_full' found, creating new files

INFO:invert4geom:starting stochastic uncertainty analysis at run 0 of 20
saving results to pickle files with prefix: '../tmp/uncertainty_full'

14.7. Comparing results

results = [
    uncert_results,
    uncert_regional_results,
    uncert_grav_results,
    uncert_constraints_results,
    uncert_density_results,
    uncert_damping_results,
]
names = [
    "full",
    "regional",
    "grav",
    "constraints",
    "density",
    "damping",
]
# get cell-wise stats for each ensemble
stats = []
for r in results:
    ds = uncertainty.merged_stats(
        results=r,
        plot=False,
        constraints_df=constraint_points,
        weight_by="constraints",
        # region=vd.pad_region(region, -spacing),
    )
    stats.append(ds)

# get the standard deviation of the ensemble of ensembles
stdevs = []
for i, s in enumerate(stats):
    stdevs.append(s.weighted_stdev.rename(f"{names[i]}_stdev"))

merged = xr.merge(stdevs)
merged

<xarray.Dataset> Size: 62kB
Dimensions:            (northing: 31, easting: 41)
Coordinates:
  * northing           (northing) float64 248B 0.0 1e+03 2e+03 ... 2.9e+04 3e+04
  * easting            (easting) float64 328B 0.0 1e+03 2e+03 ... 3.9e+04 4e+04
Data variables:
    full_stdev         (northing, easting) float64 10kB 51.98 58.67 ... 60.38
    regional_stdev     (northing, easting) float64 10kB 37.88 30.18 ... 62.01
    grav_stdev         (northing, easting) float64 10kB 2.297 1.83 ... 1.864
    constraints_stdev  (northing, easting) float64 10kB 11.16 9.923 ... 41.75
    density_stdev      (northing, easting) float64 10kB 0.1882 0.7501 ... 2.177
    damping_stdev      (northing, easting) float64 10kB 0.7195 1.011 ... 15.44

xarray.Dataset

• Dimensions:
  • northing: 31
  • easting: 41
• Coordinates: (2)
  • northing
    (northing)
    float64
    0.0 1e+03 2e+03 ... 2.9e+04 3e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000.])

  • easting
    (easting)
    float64
    0.0 1e+03 2e+03 ... 3.9e+04 4e+04

    array([    0.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000.,  7000.,  8000.,
            9000., 10000., 11000., 12000., 13000., 14000., 15000., 16000., 17000.,
           18000., 19000., 20000., 21000., 22000., 23000., 24000., 25000., 26000.,
           27000., 28000., 29000., 30000., 31000., 32000., 33000., 34000., 35000.,
           36000., 37000., 38000., 39000., 40000.])

• Data variables: (6)
  • full_stdev
    (northing, easting)
    float64
    51.98 58.67 39.57 ... 50.56 60.38

    array([[51.98424369, 58.66819901, 39.56606085, ..., 21.07919731,
            21.14563055, 33.66247   ],
           [60.47170504, 29.76524111, 38.49568085, ..., 11.20783505,
            10.49729871, 18.73280448],
           [43.66468469, 42.21300021, 28.74170413, ..., 10.64297809,
            10.3974995 , 19.4113318 ],
           ...,
           [26.30711055, 19.15152739, 16.64739967, ..., 25.05815813,
            28.20501806, 37.46572548],
           [26.61038441, 17.03328588, 20.64354523, ..., 31.59410097,
            27.5855634 , 42.87598651],
           [37.52010777, 24.66139352, 23.09631389, ..., 49.81793113,
            50.56128903, 60.37586063]], shape=(31, 41))

  • regional_stdev
    (northing, easting)
    float64
    37.88 30.18 26.97 ... 54.25 62.01

    array([[37.88161875, 30.17949666, 26.9715422 , ...,  5.97441215,
             6.27232081,  8.69407848],
           [32.21541019, 23.94806914, 21.21093661, ...,  2.94569405,
             2.42072896,  5.025056  ],
           [31.45523929, 23.58898832, 20.99994828, ...,  2.60479397,
             2.54838683,  4.19779659],
           ...,
           [16.69735107,  9.66166958,  6.12670853, ..., 30.43158569,
            31.06687124, 40.97677159],
           [15.03877287,  8.06228705,  4.49641192, ..., 35.20313391,
            35.61036427, 46.20638612],
           [13.81917139,  6.42725384,  2.21136405, ..., 54.03238447,
            54.25478056, 62.01487034]], shape=(31, 41))

  • grav_stdev
    (northing, easting)
    float64
    2.297 1.83 1.433 ... 3.386 1.864

    array([[2.29666469, 1.82988588, 1.43261486, ..., 3.21542788, 3.33770351,
            2.83577783],
           [1.84820495, 1.87069653, 1.76719364, ..., 1.90295336, 3.44233328,
            2.31951759],
           [1.63498072, 1.91735682, 2.66488596, ..., 3.21074535, 2.14493836,
            2.79689135],
           ...,
           [1.04568776, 1.62670085, 1.41280262, ..., 1.8247736 , 2.85142799,
            2.70466368],
           [1.43354854, 1.27000493, 1.65577356, ..., 2.4290159 , 1.37157015,
            2.34006217],
           [1.99864248, 0.63309974, 1.6863621 , ..., 2.13524423, 3.38589774,
            1.86430744]], shape=(31, 41))

  • constraints_stdev
    (northing, easting)
    float64
    11.16 9.923 9.91 ... 32.18 41.75

    array([[11.15511382,  9.9225396 ,  9.90990461, ..., 16.88455655,
            16.76091597, 23.24389498],
           [ 9.21200322,  7.71120685,  7.69003927, ..., 10.76288723,
             9.79138872, 16.8838483 ],
           [ 8.82741022,  7.46086529,  7.47131898, ..., 10.99723544,
            10.67117192, 16.7639268 ],
           ...,
           [12.83955936,  9.83007366,  9.09074204, ..., 24.04201314,
            22.79793267, 32.77255713],
           [12.82126044,  9.70705253,  9.11746882, ..., 22.50617151,
            20.53234587, 31.55138763],
           [14.64747546, 11.38817688, 10.46921729, ..., 32.85365352,
            32.18077752, 41.75167252]], shape=(31, 41))

  • density_stdev
    (northing, easting)
    float64
    0.1882 0.7501 ... 1.809 2.177

    array([[0.18822361, 0.75013762, 0.6420553 , ..., 1.76032401, 1.49376853,
            2.34251718],
           [0.84089359, 1.32100449, 1.00268999, ..., 0.55201735, 0.18615004,
            1.09488924],
           [0.57120728, 0.98207758, 0.55225209, ..., 0.60038068, 0.32025828,
            1.20005564],
           ...,
           [1.08691835, 0.31156784, 0.54433415, ..., 0.67791933, 0.83132039,
            1.56448247],
           [0.80281631, 0.78811659, 0.89618398, ..., 0.50504219, 0.57461333,
            1.41567492],
           [1.59139497, 0.08330434, 0.45273228, ..., 1.64797236, 1.80869414,
            2.17730765]], shape=(31, 41))

  • damping_stdev
    (northing, easting)
    float64
    0.7195 1.011 1.458 ... 16.61 15.44

    array([[ 0.71953067,  1.01054024,  1.45790426, ...,  3.7656478 ,
             3.80156241,  3.40953163],
           [ 0.9020595 ,  0.86170723,  1.54384256, ...,  4.13707859,
             4.09566923,  3.74215423],
           [ 1.6206997 ,  1.56871927,  1.45823686, ...,  4.06491409,
             4.05873329,  3.63583263],
           ...,
           [ 6.95495365,  3.42935783,  2.37496715, ..., 16.3614466 ,
            16.50044764, 15.42715786],
           [ 6.57938533,  6.04607051,  2.75745453, ..., 16.86154954,
            18.67108707, 15.83453655],
           [15.17622766,  6.82171546,  4.05106829, ..., 16.13621826,
            16.61098054, 15.43779734]], shape=(31, 41))

• Indexes: (2)
  • northing
    PandasIndex

    PandasIndex(Index([    0.0,  1000.0,  2000.0,  3000.0,  4000.0,  5000.0,  6000.0,  7000.0,
            8000.0,  9000.0, 10000.0, 11000.0, 12000.0, 13000.0, 14000.0, 15000.0,
           16000.0, 17000.0, 18000.0, 19000.0, 20000.0, 21000.0, 22000.0, 23000.0,
           24000.0, 25000.0, 26000.0, 27000.0, 28000.0, 29000.0, 30000.0],
          dtype='float64', name='northing'))

  • easting
    PandasIndex

    PandasIndex(Index([    0.0,  1000.0,  2000.0,  3000.0,  4000.0,  5000.0,  6000.0,  7000.0,
            8000.0,  9000.0, 10000.0, 11000.0, 12000.0, 13000.0, 14000.0, 15000.0,
           16000.0, 17000.0, 18000.0, 19000.0, 20000.0, 21000.0, 22000.0, 23000.0,
           24000.0, 25000.0, 26000.0, 27000.0, 28000.0, 29000.0, 30000.0, 31000.0,
           32000.0, 33000.0, 34000.0, 35000.0, 36000.0, 37000.0, 38000.0, 39000.0,
           40000.0],
          dtype='float64', name='easting'))

• Attributes: (0)

titles = [
    "True error",
    "Total uncertainty",
    "Uncertainty from regional estimation",
    "Uncertainty from gravity",
    "Uncertainty from constraints",
    "Uncertainty from density",
    "Uncertainty from damping",
]
grids = list(merged.data_vars.values())

grids.insert(0, np.abs(stats[0].weighted_mean - true_topography))

cpt_lims = polar_utils.get_min_max(
    grids[0],
    robust=True,
)

fig_height = 9
for i, g in enumerate(grids):
    xshift_amount = 1
    if i == 0:
        fig = None
        origin_shift = "initialize"
    elif i == 3:
        origin_shift = "both"
        xshift_amount = -2.5
    else:
        origin_shift = "x"

    fig = maps.plot_grd(
        grid=g,
        fig_height=fig_height,
        title=titles[i],
        title_font="16p,Helvetica,black",
        cmap="thermal",
        cpt_lims=cpt_lims,
        robust=True,
        cbar_label=f"standard deviation (m), mean: {int(np.nanmean(g))}",
        hist=True,
        hist_bin_num=50,
        cbar_yoffset=1,
        fig=fig,
        origin_shift=origin_shift,
        xshift_amount=xshift_amount,
    )
    if i == 0:
        # plot profile location, and endpoints on map
        start = [region[0], region[3]]
        stop = [region[1], region[2]]
        fig.plot(
            vd.line_coordinates(start, stop, size=100),
            pen="2p,black",
        )
        fig.text(
            x=start[0],
            y=start[1],
            text="A",
            fill="white",
            font="12p,Helvetica,black",
            justify="CM",
            clearance="+tO",
            no_clip=True,
        )
        fig.text(
            x=stop[0],
            y=stop[1],
            text="A' ",
            fill="white",
            font="12p,Helvetica,black",
            justify="CM",
            clearance="+tO",
            no_clip=True,
        )
    fig.plot(
        x=constraint_points.easting,
        y=constraint_points.northing,
        style="x.2c",
        pen="1.5p,white",
        fill="white",
    )
    fig.text(
        position="TL",
        text=f"{string.ascii_lowercase[i]}",
        fill="white",
        pen=True,
        font="16p,Helvetica,black",
        offset="j.6/.2",
        clearance="+tO",
        no_clip=True,
    )
fig.show()

data_dict = profiles.make_data_dict(
    names=titles,
    grids=grids,
    colors=[
        "red",
        "black",
        "blue",
        "magenta",
        "cyan",
        "green",
        "purple",
    ],
)

fig, df_data = profiles.plot_data(
    "points",
    start=[region[0], region[3]],
    stop=[region[1], region[2]],
    num=10000,
    fig_height=4,
    fig_width=15,
    data_dict=data_dict,
    data_legend_loc="jTR+jTL",
    data_legend_box="+gwhite",
    data_buffer=0.01,
    data_frame=["neSW", "xafg+lDistance (m)", "yag+luncertainty (m)"],
    # data_pen_style=[None,None,"4_2:2p"],
    # data_pen_thickness=[1, 1.5, 1],
    share_yaxis=True,
    start_label="A",
    end_label="A' ",
)
fig.show()