Skip to content

Commit

Permalink
Include Least Squared Correlated TC in model tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@hectornieto
hectornieto committed Feb 12, 2021
1 parent 0c63fef commit 22eee76
Showing 1 changed file with 93 additions and 12 deletions.
105 changes: 93 additions & 12 deletions test/test_ctc.py
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,19 @@
import matplotlib.pyplot as plt
import numpy as np
from model_evaluation import triple_collocation as tc
from matplotlib.lines import Line2D

sampling_sizes = [50, 100, 500, 1000]
simulations = 100000
rho_steps = 30
rho_steps = 10
rhos = np.linspace(0.01, 0.99, rho_steps)
# [sigma_1, sigma_2, sigma_3]
case_dict = {"small_uncorrelated": [0.5, 0.25, 0.1],
"equal": [0.5, 0.5, 0.5],
"large_uncorrelated": [0.1, 0.25, 0.5]}

bs = np.array([1., 1., 1.])

for case, stds in case_dict.items():
fig, axs = plt.subplots(nrows=len(sampling_sizes), ncols=3,
figsize=(15, len(sampling_sizes) * 5),
Expand All @@ -23,6 +26,16 @@
valid_etc = []
std_etc = []
bias_etc = []
valid_lsetc = []
std_lsetc = []
bias_lsetc = []
valid_rho_xy = []
bias_rho_xy = []
std_rho_xy = []
valid_rho_xy_lsetc = []
bias_rho_xy_lsetc = []
std_rho_xy_lsetc = []

for rho in rhos:
print(f"Processing {case} with N={sampling_size} "
f"and error correlation={rho:4.2f}")
Expand All @@ -31,67 +44,135 @@
noise_tuple = np.random.multivariate_normal([0, 0],
cov_matrix,
size=true.shape)
x = true + noise_tuple[:, :, 0]
y = true + noise_tuple[:, :, 1]
x = bs[0] * true + noise_tuple[:, :, 0]
y = bs[1] * true + noise_tuple[:, :, 1]
z = true + np.random.normal(0, stds[2], size=true.shape)

q_hat = tc.covariance_matrix_vec(x, y, z)
stderr, rho_xy = tc.ctc(q_hat)
stderr_etc, rho_etc, snr_db, sensitivity = tc.etc(q_hat)
cov_matrix = tc.covariance_matrix_vec(x, y, z)
factors = tc.scaling_factors(cov_matrix)
y_b = tc.apply_calibration(x, y, factors[1])

cov_matrix_b = tc.covariance_matrix_vec(x, y_b, z)

stderr, rho_xy = tc.ctc(cov_matrix_b)
stderr = stderr / factors

stderr_lsetc, rho_xy_lsetc = tc.lsetc(cov_matrix_b)
stderr_lsetc = stderr_lsetc / factors

stderr_etc, rho_etc, snr_db, sensitivity = tc.etc(cov_matrix)

valid = np.sum(np.isfinite(stderr), axis=1)
valid_ctc.append(valid / simulations)
valid = np.sum(np.isfinite(stderr_etc), axis=1)
valid_etc.append(valid / simulations)
valid = np.sum(np.isfinite(stderr_lsetc), axis=1)
valid_lsetc.append(valid / simulations)
valid = np.sum(np.isfinite(rho_xy), axis=0)
valid_rho_xy.append(valid / simulations)
valid = np.sum(np.isfinite(rho_xy_lsetc), axis=0)
valid_rho_xy_lsetc.append(valid / simulations)

bias_ctc.append((np.nanmean(stderr, axis=1) - stds))
bias_etc.append((np.nanmean(stderr_etc, axis=1) - stds))
bias_ctc.append((np.nanmean(stderr, axis=1) - stds) / stds)
bias_etc.append((np.nanmean(stderr_etc, axis=1) - stds) / stds)
bias_lsetc.append((np.nanmean(stderr_lsetc, axis=1) - stds) / stds)
bias_rho_xy.append((np.nanmean(rho_xy, axis=0) - rho) / rho)
bias_rho_xy_lsetc.append((np.nanmean(rho_xy_lsetc, axis=0) - rho) / rho)

std_ctc.append(np.nanstd(stderr, axis=1))
std_etc.append(np.nanstd(stderr_etc, axis=1))
std_lsetc.append(np.nanstd(stderr_lsetc, axis=1))
std_rho_xy.append(np.nanstd(rho_xy, axis=0))
std_rho_xy_lsetc.append(np.nanstd(rho_xy_lsetc, axis=0))

valid_ctc = np.asarray(valid_ctc).T
valid_etc = np.asarray(valid_etc).T
valid_lsetc = np.asarray(valid_lsetc).T
std_ctc = np.asarray(std_ctc).T
std_etc = np.asarray(std_etc).T
std_lsetc = np.asarray(std_lsetc).T
std_rho_xy = np.asarray(std_rho_xy)
std_rho_xy_lsetc = np.asarray(std_rho_xy_lsetc)
bias_ctc = np.array(bias_ctc).T
bias_etc = np.array(bias_etc).T
bias_lsetc = np.array(bias_lsetc).T
bias_rho_xy = np.asarray(bias_rho_xy)
bias_rho_xy_lsetc = np.asarray(bias_rho_xy_lsetc)

# Fraction of valid retrievals
axs[i, 0].plot(rhos, valid_ctc[0], "k-", label="CTC x")
axs[i, 0].plot(rhos, valid_ctc[1], "b-", label="CTC y")
axs[i, 0].plot(rhos, valid_ctc[2], "r-", label="CTC z")
axs[i, 0].plot(rhos, valid_lsetc[0], "k--", label="LSE x")
axs[i, 0].plot(rhos, valid_lsetc[1], "b--", label="LSE y")
axs[i, 0].plot(rhos, valid_lsetc[2], "r--", label="LSE z")
axs[i, 0].plot(rhos, valid_etc[0], "k:", label="ETC x")
axs[i, 0].plot(rhos, valid_etc[1], "b:", label="ETC y")
axs[i, 0].plot(rhos, valid_etc[2], "r:", label="ETC z")
axs[i, 0].legend()
axs[i, 0].plot(rhos, valid_rho_xy, "g-",
label=r"CTC $\rho_{\delta_x,\delta_y}$")
axs[i, 0].plot(rhos, valid_rho_xy_lsetc, "g--",
label=r"LSE $\rho_{\delta_x,\delta_y}$")

axs[i, 0].set_ylabel(f"N={sampling_size}")
axs[i, 0].set_ylim((0, 1))

# Normalized Mean
axs[i, 1].plot(rhos, bias_ctc[0], "k-", label="CTC x")
axs[i, 1].plot(rhos, bias_ctc[1], "b-", label="CTC y")
axs[i, 1].plot(rhos, bias_ctc[2], "r-", label="CTC z")
axs[i, 1].plot(rhos, bias_lsetc[0], "k--", label="LSE x")
axs[i, 1].plot(rhos, bias_lsetc[1], "b--", label="LSE y")
axs[i, 1].plot(rhos, bias_lsetc[2], "r--", label="LSE z")
axs[i, 1].plot(rhos, bias_etc[0], "k:", label="ETC x")
axs[i, 1].plot(rhos, bias_etc[1], "b:", label="ETC y")
axs[i, 1].plot(rhos, bias_etc[2], "r:", label="ETC z")
axs[i, 1].set_ylim((None, None))
axs[i, 1].plot(rhos, bias_rho_xy, "g-", label=r"CTC $\rho_{\delta_x,\delta_y}$")
axs[i, 1].plot(rhos, bias_rho_xy_lsetc, "g--",
label=r"LSE $\rho_{\delta_x,\delta_y}$")
axs[i, 1].set_ylim((-1, 1))

# Normalized Uncertainity
axs[i, 2].plot(rhos, std_ctc[0], "k-", label="CTC x")
axs[i, 2].plot(rhos, std_ctc[1], "b-", label="CTC y")
axs[i, 2].plot(rhos, std_ctc[2], "r-", label="CTC z")
axs[i, 2].plot(rhos, std_lsetc[0], "k--", label="LSE x")
axs[i, 2].plot(rhos, std_lsetc[1], "b--", label="LSE y")
axs[i, 2].plot(rhos, std_lsetc[2], "r--", label="LSE z")
axs[i, 2].plot(rhos, std_etc[0], "k:", label="ETC x")
axs[i, 2].plot(rhos, std_etc[1], "b:", label="ETC y")
axs[i, 2].plot(rhos, std_etc[2], "r:", label="ETC z")
axs[i, 2].set_ylim((0, None))
axs[i, 2].plot(rhos, std_rho_xy, "g-",
label=r"CTC $\rho_{\delta_x,\delta_y}$")
axs[i, 2].plot(rhos, std_rho_xy_lsetc, "g--",
label=r"LSE $\rho_{\delta_x,\delta_y}$")
axs[i, 2].set_ylim((0, 0.5))

leg_handles = [Line2D([0], [0], c="k", ls="-", label="CTC x"),
Line2D([0], [0], c="b", ls="-", label="CTC y"),
Line2D([0], [0], c="r", ls="-", label="CTC z"),
Line2D([0], [0], c="g", ls="-",
label=r"CTC $\rho_{\delta_x,\delta_y}$"),
Line2D([0], [0], c="k", ls="--", label="LSE x"),
Line2D([0], [0], c="b", ls="--", label="LSE y"),
Line2D([0], [0], c="r", ls="--", label="LSE z"),
Line2D([0], [0], c="g", ls="--",
label=r"LSE $\rho_{\delta_x,\delta_y}$"),
Line2D([0], [0], c="k", ls=":", label="ETC x"),
Line2D([0], [0], c="b", ls=":", label="ETC y"),
Line2D([0], [0], c="r", ls=":", label="ETC z"),
]

axs[0, 0].set_title("Fraction of valid retrievals")
axs[0, 1].set_title("Normalized Mean")
axs[0, 2].set_title("Normalized Uncertainity")
axs[-1, 0].set_xlabel(r"$\rho_{12}$")
axs[-1, 1].set_xlabel(r"$\rho_{12}$")
axs[-1, 2].set_xlabel(r"$\rho_{12}$")

fig.suptitle(case)
fig.legend(handles=leg_handles, ncol=3)
plt.subplots_adjust(hspace=0)

plt.show()
plt.show()

0 comments on commit 22eee76

Please sign in to comment.