Example 8: PyTracerLab Benchmarking

In this example, we benchmark simulation results from PyTracerLab against equivalent results obtained from TracerLPM. We consider the tracer input data given in Example 3 of the TracerLPM documentation (Jurgens et al., 2012). We generated simulation results using different model units in TracerLPM and compare those results against those obtained with PyTracerLab.

import PyTracerLab.model as ism
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib as mpl
from datetime import datetime
import pandas as pd
import seaborn as sns

1. Load Tracer Input Signal

# load input series
# this would be the tracer concentration in precipitation or recharge in a
# practical problem
file_name = "TracerLPM_benchmark_input_yearly.csv"
data = np.genfromtxt(
    file_name,
    delimiter=",",
    dtype=["<U7", float],
    encoding="utf-8",
    skip_header=1
)

timestamps = np.array([datetime.strptime(row[0], r"%Y") for row in data])
input_series = np.array([row[1] for row in data], dtype=float)

# load benchmark data
file_name = "TracerLPM_benchmark_simulations_1.csv"
data = np.genfromtxt(
    file_name,
    delimiter=",",
    dtype=["<U7", float, float, float, float],
    encoding="utf-8",
    skip_header=1
)

timestamps_benchmark = np.array([datetime.strptime(row[0], r"%Y") for row in data])
data_benchmark_1 = np.array([[row[1], row[2], row[3], row[4]] for row in data], dtype=float)

# load benchmark data case 2
file_name = "TracerLPM_benchmark_simulations_2.csv"
data = np.genfromtxt(
    file_name,
    delimiter=",",
    dtype=["<U7", float, float, float, float],
    encoding="utf-8",
    skip_header=1
)

data_benchmark_2 = np.array([[row[1], row[2], row[3], row[4]] for row in data], dtype=float)

# load benchmark data case 2
file_name = "TracerLPM_benchmark_simulations_3.csv"
data = np.genfromtxt(
    file_name,
    delimiter=",",
    dtype=["<U7", float, float, float, float],
    encoding="utf-8",
    skip_header=1
)

data_benchmark_3 = np.array([[row[1], row[2], row[3], row[4]] for row in data], dtype=float)

2. Get PyTracerLab Results

2.1 Piston Flow

Case 1

# define list of result time series
pm_sims = []
pm_bms = []

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
pm_mtt_true = 15. # 15 years
m.add_unit(
    ism.PMUnit(mtt=pm_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="pm"
)

# simulate
pm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_1[:, 0],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=pm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
pm_sims.append(sim)
pm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
pm_mtt_true = 20. # 20 years
m.add_unit(
    ism.PMUnit(mtt=pm_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="pm"
)

# simulate
pm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_2[:, 0],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=pm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
pm_sims.append(sim)
pm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
pm_mtt_true = 30. # 20 years
m.add_unit(
    ism.PMUnit(mtt=pm_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="pm"
)

# simulate
pm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_3[:, 0],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=pm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
pm_sims.append(sim)
pm_bms.append(bm)

2.2 Exponential Model

Case 1

# define list of result time series
em_sims = []
em_bms = []

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
em_mtt_true = 15. # 15 years
m.add_unit(
    ism.EMUnit(mtt=em_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="em"
)

# simulate
em_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_1[:, 1],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=em_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
em_sims.append(sim)
em_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
em_mtt_true = 20. # 20 years
m.add_unit(
    ism.EMUnit(mtt=em_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="em"
)

# simulate
em_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_2[:, 1],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=em_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
em_sims.append(sim)
em_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
em_mtt_true = 30. # 20 years
m.add_unit(
    ism.EMUnit(mtt=em_mtt_true),
    fraction=1., # 100 percent of the overall response
    prefix="em"
)

# simulate
em_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_3[:, 1],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=em_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
em_sims.append(sim)
em_bms.append(bm)

2.3 Exponential Piston Flow

Case 1

# define list of result time series
epm_sims = []
epm_bms = []

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
epm_mtt_true = 15. # 15 years
epm_eta_true = 0.5 + 1. # epm_ratio + 1
m.add_unit(
    ism.EPMUnit(mtt=epm_mtt_true, eta=epm_eta_true),
    fraction=1., # 100 percent of the overall response
    prefix="epm"
)

# simulate
epm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_1[:, 2],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=epm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
epm_sims.append(sim)
epm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
epm_mtt_true = 20. # 15 years
epm_eta_true = 0.7 + 1. # epm_ratio + 1
m.add_unit(
    ism.EPMUnit(mtt=epm_mtt_true, eta=epm_eta_true),
    fraction=1., # 100 percent of the overall response
    prefix="epm"
)

# simulate
epm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_2[:, 2],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=epm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
epm_sims.append(sim)
epm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
epm_mtt_true = 30. # 15 years
epm_eta_true = 0.8 + 1. # epm_ratio + 1
m.add_unit(
    ism.EPMUnit(mtt=epm_mtt_true, eta=epm_eta_true),
    fraction=1., # 100 percent of the overall response
    prefix="epm"
)

# simulate
epm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_3[:, 2],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=epm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
epm_sims.append(sim)
epm_bms.append(bm)

2.4 Dispersion Model

Case 1

# define list of result time series
dm_sims = []
dm_bms = []

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
dm_mtt_true = 15. # 15 years
dm_dp_true = 0.5
m.add_unit(
    ism.DMUnit(mtt=dm_mtt_true, DP=dm_dp_true),
    fraction=1., # 100 percent of the overall response
    prefix="dm"
)

# simulate
dm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_1[:, 3],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=dm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
dm_sims.append(sim)
dm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
dm_mtt_true = 20. # 15 years
dm_dp_true = 0.7
m.add_unit(
    ism.DMUnit(mtt=dm_mtt_true, DP=dm_dp_true),
    fraction=1., # 100 percent of the overall response
    prefix="dm"
)

# simulate
dm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_2[:, 3],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=dm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
dm_sims.append(sim)
dm_bms.append(bm)

### define model (the true system; in practice we don't know this)
# get decay constant
# we assume a half life of 12.3 years
t_half = 12.3
lambda_ = np.log(2.0) / t_half

# create true observations using the model
# time step is 1 month
m = ism.Model(
    dt=1.0,
    lambda_=lambda_,
    input_series=input_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

# add an piston-flow unit
# define the true model parameters
dm_mtt_true = 30. # 15 years
dm_dp_true = 0.8
m.add_unit(
    ism.DMUnit(mtt=dm_mtt_true, DP=dm_dp_true),
    fraction=1., # 100 percent of the overall response
    prefix="dm"
)

# simulate
dm_sim = m.simulate()

# make pandas series
bm = pd.Series(
    data=data_benchmark_3[:, 3],
    index=timestamps_benchmark
)

sim = pd.Series(
    data=dm_sim,
    index=timestamps
).truncate(
    before=timestamps_benchmark[0],
    after=timestamps_benchmark[-1]
    )

# append to list
dm_sims.append(sim)
dm_bms.append(bm)

Create Global Figures

fig, ax = plt.subplots(
    2, 2,
    figsize=(10, 5),
    sharex="col",
    # sharey="row"
)

simulation_data = [pm_sims, em_sims, epm_sims, dm_sims]
benchmark_data = [pm_bms, em_bms, epm_bms, dm_bms]

names = ["Piston Flow Model", "Exponential Model", "Exponential Piston Flow Model", "Dispersion Model"]

colors_1 = ["red", "royalblue", "green"]
colors_2 = ["darkred", "navy", "darkgreen"]

# iterate over models / units
for i in range(4):
    unit_sims = simulation_data[i]
    unit_bms = benchmark_data[i]

    ax_ = ax.flatten()[i]

    for j in range(len(unit_sims)):
        # plot simulation
        ax_.plot(
            unit_sims[j].index,
            unit_sims[j].values,
            c=colors_1[j],
            lw=1.5,
            zorder=100
        )
        # plot benchmark
        ax_.scatter(
            unit_bms[j].index,
            unit_bms[j].values,
            c=colors_2[j],
            marker=".",
            s=20.,
            zorder=1000
        )
    # compose legend
    legend_elements = [
        Line2D([0], [0], color="k", lw=1, label="PyTracerLab"),
        Line2D([0], [0], marker=".", color="w",
               markerfacecolor="k", markeredgecolor="k", label="TracerLPM")
               ]
    ax_.legend(handles=legend_elements, loc="upper right")
    ax_.set_title(names[i])
    ax_.set_ylabel(r"Tritium $ [TU] $")
    ax_.grid(True, zorder=0, alpha=0.7)
    ax_.set_xlim(unit_sims[j].index[0], unit_sims[j].index[-1])
    ax_.set_ylim(0.)

    if i > 1:
        ax_.set_xlabel("Time")

plt.savefig("simulation_benchmark.png", dpi=400, bbox_inches="tight")

Calibration Benchmark

# load input series with 6 month resolution
# this would be the tracer concentration in precipitation or recharge in a
# practical problem
file_name = "TracerLPM_benchmark_input_6monthly.csv"
data = np.genfromtxt(
    file_name,
    delimiter=",",
    dtype=["<U7", float],
    encoding="utf-8",
    skip_header=1
)

timestamps = np.array([datetime.strptime(row[0], r"%Y-%m") for row in data])
input_series = np.array([row[1] for row in data], dtype=float)

# generate synthetic observations
# create model
t_half = 12.3 # years
lambda_ = np.log(2.0) / t_half

### define model (use the same structure / units as the true model)
# time step is 0.5 years
m_obs = ism.Model(
    dt=.5,
    lambda_=lambda_,
    input_series=input_series,
    target_series=None,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

dmfrac_1 = .2
dmfrac_2 = 1 - dmfrac_1

# add a dispersion model unit
# define the true model parameters
dm1_mtt_init = 2.
dm1_DP_init = .1
m_obs.add_unit(
    ism.DMUnit(mtt=dm1_mtt_init, DP=dm1_DP_init),
    fraction=dmfrac_1,
    bounds=[(0., 20.), (.001, 10.)],
    prefix="dm1"
)

# add a second dispersion model unit
# define the true model parameters
dm2_mtt_init = 25.
dm2_DP_init = .1
m_obs.add_unit(
    ism.DMUnit(mtt=dm2_mtt_init, DP=dm2_DP_init),
    fraction=dmfrac_2,
    bounds=[(10., 100.), (.001, 10.)],
    prefix="dm2"
)

# create a solver
sim_reference = m_obs.simulate()

# get observations from reference simulation
# define number of observations to get
n_obs = 10 # 20

# n_obs=10, seed=1234567

# set random seed
np.random.seed(1234) # 1234

# get random observations from the reference simulation
obs_idx = np.random.choice([i for i in range(100, len(input_series))], n_obs, replace=False)
# obs_idx = [i for i in range(120, 240)][::5][::3]
# n_obs = len(obs_idx)
obs_timestamps = timestamps[obs_idx]
obs_values = sim_reference[obs_idx]

# add noise to observations
# define observation noise level
obs_noise = 5
obs_values += np.random.normal(0, obs_noise, size=n_obs)

# make series we can later use in the model (has to be the same length as
# the input series, filled with NaN-values where we do not have any
# observations)
obs_series = np.full(len(input_series), np.nan)
obs_series[obs_idx] = obs_values

# plot noisy observations and reference simulation
fig, ax = plt.subplots(2, 1, figsize=(10, 7))

start = 100
# ax.plot(timestamps, input_series, label="Input")
ax[0].plot(timestamps[start:], sim_reference[start:], c="k", label="Reference Simulation")
# ax.plot(timestamps, input_series, c="grey", label="Reference Simulation")
ax[0].scatter(
    timestamps[start:], obs_series[start:],
    marker="o", facecolor="w",
    edgecolor="k", s=20,
    zorder=100, alpha=0.8,
    lw=1.
)

step_limit = 100
dt = 0.5
t_plot = [dt * i for i in range(step_limit)]
ttds = m_obs.get_ttds()
ax[1].plot(
    t_plot, 
    ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
    ttds["fractions"][1] * ttds["distributions"][1][:step_limit]
)

# ax.set_yscale("log")

[<matplotlib.lines.Line2D at 0x2506ff6c800>]

Original Model Hyperparameters

# those settings are reused for all MCMC runs of all models
n_samples = 2000
burn_in = 2000
thin = 1
n_cr = 4
sigma = np.sqrt(obs_noise)

# those settings are reused for all DE runs of all models
maxiter = 1000
popsize = 20

# generate synthetic observations
# create model
t_half = 12.3 # years
lambda_ = np.log(2.0) / t_half

### define model (use the same structure / units as the true model)
# time step is 0.5 years
m1 = ism.Model(
    dt=.5,
    lambda_=lambda_,
    input_series=input_series,
    target_series=obs_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

dmfrac_11 = .2
dmfrac_12 = 1 - dmfrac_11

# add a dispersion model unit
# define the true model parameters
dm1_mtt_init = 5.
dm1_DP_init = .5
m1.add_unit(
    ism.DMUnit(mtt=dm1_mtt_init, DP=dm1_DP_init),
    fraction=dmfrac_11,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm1"
)

# add a second dispersion model unit
# define the true model parameters
dm2_mtt_init = 50.
dm2_DP_init = .5
m1.add_unit(
    ism.DMUnit(mtt=dm2_mtt_init, DP=dm2_DP_init),
    fraction=dmfrac_12,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm2"
)

# create a solver
m1_solver = ism.Solver(m1)
res_1_x, res_1 = m1_solver.differential_evolution(
    maxiter=maxiter,
    popsize=popsize,
    sigma=sigma,
    random_state=42
)

res_mcmc_1 = m1_solver.dream_sample(
    n_samples=n_samples,
    n_chains=None,
    burn_in=burn_in,
    thin=thin,
    cr=[i / n_cr for i in range(1, n_cr + 1)],
    sigma=sigma,
    start=[1., 5., 1., 5.],
    return_sim=True,
    random_state=42
)

# plot chains to inspect convergence
fig, ax = plt.subplots(4, 1, figsize=(6, 6), sharex=True)

for i in range(res_mcmc_1["samples_chain"].shape[2]):
    for j in range(res_mcmc_1["samples_chain"].shape[0]):
        ax[i].plot(res_mcmc_1["samples_chain"][j, :, i], c="k", lw=.5, alpha=0.4)
    ax[i].set_ylabel(m1.param_keys(free_only=True)[i])
ax[-1].set_xlabel("Steps")

# print gelman rubin convergence statistics
print(res_mcmc_1["gelman_rubin"])

{'dm1.DP': 1.0351106346126266, 'dm1.mtt': 1.030493503812539, 'dm2.DP': 1.0199520327803513, 'dm2.mtt': 1.025483443900886}

Mixing Ratio Case 2

# generate synthetic observations
# create model
t_half = 12.3 # years
lambda_ = np.log(2.0) / t_half

### define model (use the same structure / units as the true model)
# time step is 0.5 years
m2 = ism.Model(
    dt=.5,
    lambda_=lambda_,
    input_series=input_series,
    target_series=obs_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

dmfrac_21 = .1
dmfrac_22 = 1 - dmfrac_21

# add a dispersion model unit
# define the true model parameters
dm1_mtt_init = 5.
dm1_DP_init = .5
m2.add_unit(
    ism.DMUnit(mtt=dm1_mtt_init, DP=dm1_DP_init),
    fraction=dmfrac_21,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm1"
)

# add a second dispersion model unit
# define the true model parameters
dm2_mtt_init = 50.
dm2_DP_init = .5
m2.add_unit(
    ism.DMUnit(mtt=dm2_mtt_init, DP=dm2_DP_init),
    fraction=dmfrac_22,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm2"
)

# create a solver
m2_solver = ism.Solver(m2)
res_2_x, res_2 = m2_solver.differential_evolution(
    maxiter=maxiter,
    popsize=popsize,
    sigma=sigma,
    random_state=42
)

res_mcmc_2 = m2_solver.dream_sample(
    n_samples=n_samples,
    n_chains=None,
    burn_in=burn_in,
    thin=thin,
    cr=[i / n_cr for i in range(1, n_cr + 1)],
    sigma=sigma,
    start=[1., 5., 1., 5.],
    return_sim=True,
    random_state=42
)

# plot chains to inspect convergence
fig, ax = plt.subplots(4, 1, figsize=(6, 6), sharex=True)

for i in range(res_mcmc_2["samples_chain"].shape[2]):
    for j in range(res_mcmc_2["samples_chain"].shape[0]):
        ax[i].plot(res_mcmc_2["samples_chain"][j, :, i], c="k", lw=.5, alpha=0.4)
    ax[i].set_ylabel(m2.param_keys(free_only=True)[i])
ax[-1].set_xlabel("Steps")

# print gelman rubin convergence statistics
print(res_mcmc_2["gelman_rubin"])

{'dm1.DP': 1.0333443421432507, 'dm1.mtt': 1.0103574959299573, 'dm2.DP': 1.0280525422613984, 'dm2.mtt': 1.0303220888517228}

Mixing Ratio Case 3

# generate synthetic observations
# create model
t_half = 12.3 # years
lambda_ = np.log(2.0) / t_half

### define model (use the same structure / units as the true model)
# time step is 0.5 years
m3 = ism.Model(
    dt=.5,
    lambda_=lambda_,
    input_series=input_series,
    target_series=obs_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

dmfrac_31 = .15
dmfrac_32 = 1 - dmfrac_31

# add a dispersion model unit
# define the true model parameters
dm1_mtt_init = 5.
dm1_DP_init = .5
m3.add_unit(
    ism.DMUnit(mtt=dm1_mtt_init, DP=dm1_DP_init),
    fraction=dmfrac_31,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm1"
)

# add a second dispersion model unit
# define the true model parameters
dm2_mtt_init = 50.
dm2_DP_init = .5
m3.add_unit(
    ism.DMUnit(mtt=dm2_mtt_init, DP=dm2_DP_init),
    fraction=dmfrac_32,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm2"
)

# create a solver
m3_solver = ism.Solver(m3)
res_3_x, res_3 = m3_solver.differential_evolution(
    maxiter=maxiter,
    popsize=popsize,
    sigma=sigma,
    random_state=42
)

res_mcmc_3 = m3_solver.dream_sample(
    n_samples=n_samples,
    n_chains=None,
    burn_in=burn_in,
    thin=thin,
    cr=[i / n_cr for i in range(1, n_cr + 1)],
    sigma=sigma,
    start=[1., 5., 1., 5.],
    return_sim=True,
    random_state=42
)

# plot chains to inspect convergence
fig, ax = plt.subplots(4, 1, figsize=(6, 6), sharex=True)

for i in range(res_mcmc_3["samples_chain"].shape[2]):
    for j in range(res_mcmc_3["samples_chain"].shape[0]):
        ax[i].plot(res_mcmc_3["samples_chain"][j, :, i], c="k", lw=.5, alpha=0.4)
    ax[i].set_ylabel(m3.param_keys(free_only=True)[i])
ax[-1].set_xlabel("Steps")

# print gelman rubin convergence statistics
print(res_mcmc_3["gelman_rubin"])

{'dm1.DP': 1.0138866562346596, 'dm1.mtt': 1.0103786244130015, 'dm2.DP': 1.0124930668853718, 'dm2.mtt': 1.0134066969566984}

Mixing Ratio Case 4

# generate synthetic observations
# create model
t_half = 12.3 # years
lambda_ = np.log(2.0) / t_half

### define model (use the same structure / units as the true model)
# time step is 0.5 years
m4 = ism.Model(
    dt=.5,
    lambda_=lambda_,
    input_series=input_series,
    target_series=obs_series,
    steady_state_input=6.,
    n_warmup_half_lives=10
)

dmfrac_41 = .3
dmfrac_42 = 1 - dmfrac_41

# add a dispersion model unit
# define the true model parameters
dm1_mtt_init = 5.
dm1_DP_init = .5
m4.add_unit(
    ism.DMUnit(mtt=dm1_mtt_init, DP=dm1_DP_init),
    fraction=dmfrac_41,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm1"
)

# add a second dispersion model unit
# define the true model parameters
dm2_mtt_init = 50.
dm2_DP_init = .5
m4.add_unit(
    ism.DMUnit(mtt=dm2_mtt_init, DP=dm2_DP_init),
    fraction=dmfrac_42,
    bounds=[(1., 100.), (.001, 5.)],
    prefix="dm2"
)

# create a solver
m4_solver = ism.Solver(m4)
res_4_x, res_4 = m4_solver.differential_evolution(
    maxiter=maxiter,
    popsize=popsize,
    sigma=sigma,
    random_state=42
)

res_mcmc_4 = m4_solver.dream_sample(
    n_samples=n_samples,
    n_chains=None,
    burn_in=burn_in,
    thin=thin,
    cr=[i / n_cr for i in range(1, n_cr + 1)],
    sigma=sigma,
    start=[1., 5., 1., 5.],
    return_sim=True,
    random_state=42
)

# plot chains to inspect convergence
fig, ax = plt.subplots(4, 1, figsize=(6, 6), sharex=True)

for i in range(res_mcmc_4["samples_chain"].shape[2]):
    for j in range(res_mcmc_4["samples_chain"].shape[0]):
        ax[i].plot(res_mcmc_4["samples_chain"][j, :, i], c="k", lw=.5, alpha=0.4)
    ax[i].set_ylabel(m4.param_keys(free_only=True)[i])
ax[-1].set_xlabel("Steps")

# print gelman rubin convergence statistics
print(res_mcmc_4["gelman_rubin"])

{'dm1.DP': 1.008582838160302, 'dm1.mtt': 1.0080960695184484, 'dm2.DP': 1.0070163323853205, 'dm2.mtt': 1.0118377270130168}

def nse(obs, sim):
    obs_ = obs.reshape(-1, 1)
    sim_ = sim.reshape(-1, 1)
    mask = ~np.isnan(obs_) & ~np.isnan(sim_)
    resid = sim_ - obs_
    dev = obs_ - np.nanmean(obs_)
    nse_ = 1 - np.sum(resid[mask] ** 2) / np.sum(dev[mask] ** 2)
    return nse_

print(nse(obs_series, res_1))

0.9994503948336053

fig, ax = plt.subplots(2, 1, figsize=(5, 6),
                       gridspec_kw={"height_ratios": [2, 1.2], "hspace": .5, "wspace": .1},
                       sharey="row", sharex="row")

# plot results
start = 100
end = 240

c1 = "royalblue"
c2 = "red"

ax[0].plot(
    timestamps[start:end],
    res_1[start:end],
    c=c2,
    ls=":",
    lw=3.,
    zorder=1000,
    label="Differential Evolution"
)

sim_data = res_mcmc_1["sims"]
s = sim_data.reshape((-1, sim_data.shape[-1]))
q = np.quantile(s, [0.01, 0.5, 0.99], axis=0)
ax[0].plot(
    timestamps[start:end],
    q[1, start:end],
    c=c1,
    ls="-",
    lw=1.5,
    label="DREAM"
)
ax[0].fill_between(
    timestamps[start:end],
    q[0, start:end], q[2, start:end],
    color=c1,
    alpha=0.2
)

ax[0].plot(
    timestamps[start:end],
    sim_reference[start:end],
    c="k",
    ls="-",
    lw=1.,
    zorder=1000,
    label="True Model"
)

ax[0].scatter(
    timestamps, obs_series,
    marker="o", facecolor="w",
    edgecolor="k", s=40,
    zorder=10000, alpha=0.8,
    lw=1.,
    label="Observations"
)

ax[0].set_title(r"Simulated and Observed Concentrations", fontsize=11)
ax[0].grid(True, alpha=0.7, zorder=0)
ax[0].set_xlim([timestamps[start], timestamps[end-1]])
ax[0].set_ylim(0)
ax[0].set_xlabel("Time")
ax[0].set_ylabel(r"Tritium $ [TU] $")
ax[0].legend()

# plot ttds
step_limit = 70
dt = 0.5
t_plot = [dt * i for i in range(step_limit)]

# DE results
for item in res_1_x.items():
    m1.set_param(key=item[0], value=item[1])
ttds = m1.get_ttds()
ax[1].plot(
    t_plot, 
    ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
    ttds["fractions"][1] * ttds["distributions"][1][:step_limit],
    lw=3.,
    c=c2,
    ls=":", zorder=1000
)

# MCMC results
# get samples from MCMC
s = res_mcmc_1["samples"]
# get random subset
s = s[np.random.choice(s.shape[0], 2000, replace=True)]

# iterate over samples
tt_dists = []
for sample in s:
    # set model parameters
    for item in zip(m1.param_keys(free_only=True), sample):
        m1.set_param(key=item[0], value=item[1])
    # get travel time distribution
    ttds = m1.get_ttds()
    tt_dists.append(
        ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
        ttds["fractions"][1] * ttds["distributions"][1][:step_limit]
    )
# plot
ax[1].plot(
    t_plot, 
    np.median(tt_dists, axis=0),
    alpha=1.,
    c=c1,
    lw=1.5,
    ls="-"
)
ax[1].fill_between(
    t_plot, 
    np.quantile(tt_dists, 0.99, axis=0),
    np.quantile(tt_dists, 0.01, axis=0),
    alpha=.2,
    color=c1
)

ttds = m_obs.get_ttds()
ax[1].plot(
    t_plot,
    ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
    ttds["fractions"][1] * ttds["distributions"][1][:step_limit],
    c="k", lw=1.,
)

ax[1].set_xlim(0., step_limit * dt - dt)
ax[1].grid(True, alpha=0.7, zorder=0)
ax[1].set_title(r"Inferred and True Travel Time Distributions", fontsize=11)
ax[1].set_xlabel(r"Travel Time $ [a] $")
ax[1].set_ylabel(r"Fraction $ [-] $")
# ax[1].set_yscale("log")
# ax[1].set_ylim(0.00045, .15)

plt.savefig("calibration_benchmark.png", dpi=400, bbox_inches="tight")

fig, ax = plt.subplots(2, 2, figsize=(10, 5),
                       gridspec_kw={"height_ratios": [2, 1.2], "hspace": .5, "wspace": .1},
                       sharey="row", sharex="row")

ax1 = ax[1, :]
ax = ax[0, :]

data_sims = [res_2, res_1, res_3, res_4]
colors = ["red", "green", "royalblue", "orange"] # mpl.colormaps["Set2"].colors # ["r", "g", "b", "m"]
linestyles = ["--", "-", "--", "--"] # ["-.", "--", "-", ":"]
lws = [1.5, 2.5, 1.5, 1.5]
start = 100
end = 240

for num, sims in enumerate(data_sims):
    ax[0].plot(
        timestamps[start:end],
        sims[start:end],
        c=colors[num],
        ls=linestyles[num],
        lw=lws[num]
    )

ax[0].scatter(
    timestamps, obs_series,
    marker="o", facecolor="w",
    edgecolor="k", s=40,
    zorder=100, alpha=0.8,
    lw=1.
)

custom_lines = [Line2D([0], [0], color=colors[0], lw=lws[0], ls=linestyles[0]),
                Line2D([0], [0], color=colors[1], lw=lws[1], ls=linestyles[1]),
                Line2D([0], [0], color=colors[2], lw=lws[2], ls=linestyles[2]),
                Line2D([0], [0], color=colors[3], lw=lws[3], ls=linestyles[3]),
                Line2D([0], [0], marker="o", color="w", markeredgewidth=1.,
                    markerfacecolor="w", markeredgecolor="k", alpha=0.8)]

ax[0].legend(
    custom_lines,
    [
    fr"{dmfrac_21:.0%} $ DM_1 $ (${res_2_x["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_22:.0%} $ DM_2 $ (${res_2_x["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, res_2):.2f}",
    fr"{dmfrac_11:.0%} $ DM_1 $ (${res_1_x["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_12:.0%} $ DM_2 $ (${res_1_x["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, res_1):.2f}",
    fr"{dmfrac_31:.0%} $ DM_1 $ (${res_3_x["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_32:.0%} $ DM_2 $ (${res_3_x["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, res_3):.2f}",
    fr"{dmfrac_41:.0%} $ DM_1 $ (${res_4_x["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_42:.0%} $ DM_2 $ (${res_4_x["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, res_4):.2f}",
    "Observations"
    ],
    loc="upper right",
    fontsize=6.5,
)
ax[0].set_title(r"Simulated and Observed Tritium Conc. - $ \mathbf{DE} $")
ax[0].grid(True, alpha=0.7, zorder=0)
ax[0].set_xlim([timestamps[start], timestamps[end-1]])
ax[0].set_xlabel("Time")
ax[0].set_ylabel(r"Tritium $ [TU] $")

mcmc_sims = [res_mcmc_2, res_mcmc_1, res_mcmc_3, res_mcmc_4]

median_sims = []
for num, sims in enumerate(mcmc_sims):
    # reshape simulations array
    chains = [1, 2, 3, 4, 5, 6, 7]
    sim_data = sims["sims"]
    s = sim_data.reshape((-1, sim_data.shape[-1]))
    q = np.quantile(s, [0.01, 0.5, 0.99], axis=0)
    median_sims.append(q[1])
    ax[1].plot(
        timestamps[start:end],
        q[1, start:end],
        c=colors[num],
        ls=linestyles[num],
        lw=lws[num]
    )
    ax[1].fill_between(
        timestamps[start:end],
        q[0, start:end], q[2, start:end],
        color=colors[num],
        alpha=0.2
    )

ax[1].scatter(
    timestamps, obs_series,
    marker="o", facecolor="w",
    edgecolor="k", s=40,
    zorder=100, alpha=0.8,
    lw=1.
)

ax[1].legend(
    custom_lines,
    [
    fr"{dmfrac_21:.0%} $ DM_1 $ (${res_mcmc_2["posterior_map"]["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_22:.0%} $ DM_2 $ (${res_mcmc_2["posterior_map"]["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, median_sims[0]):.2f}",
    fr"{dmfrac_11:.0%} $ DM_1 $ (${res_mcmc_1["posterior_map"]["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_12:.0%} $ DM_2 $ (${res_mcmc_1["posterior_map"]["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, median_sims[1]):.2f}",
    fr"{dmfrac_31:.0%} $ DM_1 $ (${res_mcmc_3["posterior_map"]["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_32:.0%} $ DM_2 $ (${res_mcmc_3["posterior_map"]["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, median_sims[2]):.2f}",
    fr"{dmfrac_41:.0%} $ DM_1 $ (${res_mcmc_4["posterior_map"]["dm1.mtt"]:.1f}a$), " + \
        fr"{dmfrac_42:.0%} $ DM_2 $ (${res_mcmc_4["posterior_map"]["dm2.mtt"]:.1f}a$), " + \
            fr"NSE = {nse(obs_series, median_sims[3]):.2f}",
    "Observations"
    ],
    loc="upper right",
    fontsize=6.5,
)
ax[1].set_title(r"Simulated and Observed Tritium Conc. - $ \mathbf{MCMC} $")
ax[1].grid(True, alpha=0.7, zorder=0)
ax[1].set_xlabel("Time")

# plot travel time distribution
step_limit = 70
dt = 0.5
t_plot = [dt * i for i in range(step_limit)]
# ax.plot(ttds["fractions"][0] * ttds["distributions"][0][:step_limit])
# ax.plot(ttds["fractions"][1] * ttds["distributions"][1][:step_limit])
for num, (ml, res_x) in enumerate(zip([m2, m1, m3, m4], [res_2_x, res_1_x, res_3_x, res_4_x])):
    # we have to set the model parameters again for each model because they
    # currently are at the last MCMC iteration
    for item in res_x.items():
        ml.set_param(key=item[0], value=item[1])
    ttds = ml.get_ttds()
    ax1[0].plot(
        t_plot, 
        ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
        ttds["fractions"][1] * ttds["distributions"][1][:step_limit],
        # marker=".",
        markeredgecolor=colors[num],
        lw=lws[num],
        c=colors[num],
        markerfacecolor="lightgrey",
        markersize=10,
        # clip_on=False,
        # zorder=1000
        ls=linestyles[num]
    )
ax1[0].set_xlim(0., step_limit * dt - dt)
ax1[0].grid(True, alpha=0.7, zorder=0)
ax1[0].set_title(r"Travel Time Distributions - $ \mathbf{DE} $")
ax1[0].set_xlabel(r"Travel Time $ [a] $")
ax1[0].set_ylabel(r"Fraction $ [-] $")
ax1[0].set_yscale("log")
ax1[0].set_ylim(0.00045, .15)

for num, (ml, s) in enumerate(zip([m2, m1, m3, m4], [res_mcmc_2, res_mcmc_1, res_mcmc_3, res_mcmc_4])):
    # get samples from MCMC
    s = s["samples"]
    # get random subset
    s = s[np.random.choice(s.shape[0], 2000, replace=True)]
    
    # iterate over samples
    tt_dists = []
    for sample in s:
        # set model parameters
        for item in zip(ml.param_keys(free_only=True), sample):
            ml.set_param(key=item[0], value=item[1])
        # get travel time distribution
        ttds = ml.get_ttds()
        tt_dists.append(
            ttds["fractions"][0] * ttds["distributions"][0][:step_limit] + \
            ttds["fractions"][1] * ttds["distributions"][1][:step_limit]
        )
    # plot
    ax1[1].plot(
        t_plot, 
        np.median(tt_dists, axis=0),
        alpha=1.,
        c=colors[num],
        lw=lws[num],
        ls=linestyles[num]
    )
    ax1[1].fill_between(
        t_plot, 
        np.quantile(tt_dists, 0.99, axis=0),
        np.quantile(tt_dists, 0.01, axis=0),
        alpha=.2,
        color=colors[num]
    )

ax1[1].grid(True, alpha=0.7, zorder=0)
ax1[1].set_title(r"Travel Time Distributions - $ \mathbf{MCMC} $")
ax1[1].set_xlabel(r"Travel Time $ [a] $")

plt.savefig("calibration_benchmark_4models.png", dpi=400, bbox_inches="tight")

# make pairplot of true model posterior

param_names_axes = [
    r"$ DM_1, \; DP \; [-] $",
    r"$ DM_1, \; t^* \; [a] $",
    r"$ DM_2, \; DP \; [-] $",
    r"$ DM_2, \; t^* \; [a] $"
]

# make dict from sample array
step = 1

# this below is actually not really necessary but we need to re-map the 
# parameter names in order for them to be used by the pairplot
sample_df = pd.DataFrame({
    param_names_axes[0]: res_mcmc_1["samples"][::step, 0],
    param_names_axes[1]: res_mcmc_1["samples"][::step, 1],
    param_names_axes[2]: res_mcmc_1["samples"][::step, 2],
    param_names_axes[3]: res_mcmc_1["samples"][::step, 3],
})

obs_df = pd.DataFrame({
    param_names_axes[0]: [m_obs.params[res_mcmc_1["param_names"][0]]["value"]],
    param_names_axes[1]: [m_obs.params[res_mcmc_1["param_names"][1]]["value"]],
    param_names_axes[2]: [m_obs.params[res_mcmc_1["param_names"][2]]["value"]],
    param_names_axes[3]: [m_obs.params[res_mcmc_1["param_names"][3]]["value"]],
})

map_df = pd.DataFrame({
    param_names_axes[0]: [res_mcmc_1["posterior_map"][res_mcmc_1["param_names"][0]]],
    param_names_axes[1]: [res_mcmc_1["posterior_map"][res_mcmc_1["param_names"][1]]],
    param_names_axes[2]: [res_mcmc_1["posterior_map"][res_mcmc_1["param_names"][2]]],
    param_names_axes[3]: [res_mcmc_1["posterior_map"][res_mcmc_1["param_names"][3]]],
})

sol_de_df = pd.DataFrame({
    param_names_axes[0]: [res_1_x[res_mcmc_1["param_names"][0]]],
    param_names_axes[1]: [res_1_x[res_mcmc_1["param_names"][1]]],
    param_names_axes[2]: [res_1_x[res_mcmc_1["param_names"][2]]],
    param_names_axes[3]: [res_1_x[res_mcmc_1["param_names"][3]]],
})

# p = sns.jointplot(sample_df, x="DP", y="MTT", kind="kde", height=3.5, cmap="binary",
#               fill=True, levels=7, linewidths=2., thresh=0.02,
#               marginal_kws={"color": "k", "fill": True})
cmap = mpl.colormaps["Blues"]
colors = cmap(np.linspace(0, 1, 32))
c_blue = colors[7]
p = sns.pairplot(sample_df, kind="hist", corner=True, height=1.5,
                 # plot_kws={"fill": True, "levels": 5, "cmap": "binary"},
                 plot_kws={"bins": 20, "cmap": "Blues"},
                 diag_kws={"bins": 20, "color": "k", "kde": True,
                           "facecolor": c_blue, "linewidths": .5},)
def off_diag(x, y, **kwargs):    
    plt.gca().scatter(
        [map_df[x.name]],
        [map_df[y.name]], s=100,
        color="r", lw=.5, zorder=1000000, marker="o", edgecolor="k")
    
    plt.gca().scatter(
        [sol_de_df[x.name]],
        [sol_de_df[y.name]], s=30,
        color="w", lw=.5, zorder=1000000, marker="v", edgecolor="k",)

    plt.gca().axvline(obs_df[x.name].values, color="k", lw=1., ls="--")
    plt.gca().axhline(obs_df[y.name].values, color="k", lw=1., ls="--")

def diag(x, **kwargs):    
    plt.gca().axvline(obs_df[x.name].values, color="k", lw=1., ls="--")

p.map_offdiag(off_diag)
p.map_diag(diag)

p.map_lower(
    sns.kdeplot,
    color="k",
    zorder=1000,
    levels=5,
    linewidths=0.5,
    thresh=0.01
)

# p.map_diag(
#     sns.kdeplot,
#     color="k",
#     zorder=10000,
#     linewidth=0.5
# )

p.savefig("joint_posterior_benchmark_1.png", dpi=600)

c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)
c:\Users\MRudolph\anaconda3\envs\isosimpy_app\Lib\site-packages\matplotlib\cbook.py:1719: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return math.isfinite(val)