Exploritory work with grace and knockout tests - moving to main

test_var.py

@@ -0,0 +1,8 @@
+import sys
+import numpy as np
+sys.path.append("work/13_knockout_test")
+from knockout_test import compute_gmsl_posterior, posterior_full
+
+print("Computing GMSL posterior full:")
+exp_mm, std_mm = compute_gmsl_posterior(posterior_full)
+print(f"Exp: {exp_mm}, Std: {std_mm}")

work/08_joint_inversion/joint_inversion.py

@@ -73,8 +73,6 @@
 
 plot(s[2], symmetric=True)
 
-plt.show()
-
 # %%
 
 ssh_altimetry = GridPoints.ocean_altimetry(fp, 5.0, 66.0)
@@ -364,6 +362,3 @@ def progress_callback(xk):
 ax6.set_title(
     "f) Posterior Expectation (Inferred from Data)"
 )
-
-
-plt.show()

work/08_joint_inversion/joint_precon_factored.py

@@ -5,6 +5,7 @@
     BlockLinearOperator,
     CGMatrixSolver,
     EigenSolver,
+    EuclideanSpace,
     GaussianMeasure,
     HilbertSpaceDirectSum,
     LinearBayesianInversion,
@@ -20,13 +21,16 @@
 )
 from tqdm import tqdm
 
+from project import colors
 from project.factored_forward_operator import (
     build_factored_forward_operator,
     build_factored_forward_operator_precon,
 )
 from project.operators import (
     ice_thickness_to_estimated_gmsl_operator,
 )
+from pygeoinf_extras import standard_dev
+from pygeoinf_extras.plots import plot_bivariate_corner
 from pyslfp_extras.altimetry import GridPoints
 from pyslfp_extras.ice_thickness import IceSheetChange
 from pyslfp_extras.ocean_dynamics import OceanDynamics
@@ -350,7 +354,7 @@ def progress_callback(xk):
     bayesian_inversion.model_posterior_measure(
         data,
         CGMatrixSolver(
-            callback=progress_callback, maxiter=300
+            callback=progress_callback, maxiter=500
         ),
         preconditioner=precon_inverse_normal_operator,
     )
@@ -684,3 +688,83 @@ def progress_callback(xk):
 fig10.savefig(
     f"{dir}/joint_precon_total_load_posterior.pdf", dpi=600
 )
+
+# %%
+# =============================================================================
+# Ice vs firn GMSL covariance (corner plot)
+# =============================================================================
+
+ice_gmsl_op = ice.ice_thickness_to_gmsl_operator
+firn_gmsl_op = ice.firn_thickness_to_gmsl_operator
+
+# Row operators mapping the full 3-component model space to scalar GMSL (mm)
+ice_gmsl_row = RowLinearOperator(
+    [
+        1000 * ice_gmsl_op,
+        ice.firn_thickness.domain.zero_operator(
+            codomain=ice_gmsl_op.codomain
+        ),
+        odt.height_measure.domain.zero_operator(
+            codomain=ice_gmsl_op.codomain
+        ),
+    ]
+)
+firn_gmsl_row = RowLinearOperator(
+    [
+        ice.ice_thickness.domain.zero_operator(
+            codomain=firn_gmsl_op.codomain
+        ),
+        1000 * firn_gmsl_op,
+        odt.height_measure.domain.zero_operator(
+            codomain=firn_gmsl_op.codomain
+        ),
+    ]
+)
+
+# Scalar posterior measures for each component
+ice_gmsl_post = model_posterior_measure.affine_mapping(
+    operator=ice_gmsl_row
+)
+firn_gmsl_post = model_posterior_measure.affine_mapping(
+    operator=firn_gmsl_row
+)
+
+# Cross-covariance via polarization identity:
+# Cov(A, B) = 0.5 * (Var(A+B) - Var(A) - Var(B))
+sum_gmsl_post = model_posterior_measure.affine_mapping(
+    operator=ice_gmsl_row + firn_gmsl_row
+)
+var_ice = standard_dev(ice_gmsl_post) ** 2
+var_firn = standard_dev(firn_gmsl_post) ** 2
+var_sum = standard_dev(sum_gmsl_post) ** 2
+cross_cov = 0.5 * (var_sum - var_ice - var_firn)
+
+# Assemble 2x2 covariance matrix and mean on a flat EuclideanSpace(2)
+mu_ice = ice_gmsl_post.expectation[0]
+mu_firn = firn_gmsl_post.expectation[0]
+mean_2d = np.array([mu_ice, mu_firn])
+cov_2d = np.array(
+    [[var_ice, cross_cov], [cross_cov, var_firn]]
+)
+
+joint_gmsl_posterior_measure = (
+    GaussianMeasure.from_covariance_matrix(
+        EuclideanSpace(2), cov_2d, expectation=mean_2d
+    )
+)
+
+true_ice_gmsl = ice_gmsl_row(model_true)[0]
+true_firn_gmsl = firn_gmsl_row(model_true)[0]
+
+fig_cov, axes_cov = plot_bivariate_corner(
+    joint_gmsl_posterior_measure,
+    true_values=np.array([true_ice_gmsl, true_firn_gmsl]),
+    labels=["Ice GMSL (mm)", "Firn GMSL (mm)"],
+    title="Joint Posterior: Ice vs Firn GMSL Contributions",
+    figsize=(6.5, 6.5),
+    pdf_colors=[colors.ice, colors.firn],
+)
+fig_cov.savefig(
+    f"{dir}/joint_precon_ice_firn_gmsl_covariance.pdf",
+    dpi=600,
+)