independent verification

.gitignore

@@ -137,3 +137,8 @@ docsrc/notebooks/.DS_Store
 
 # test
 /test
+docsrc/notebooks/.jupyter/desktop-workspaces/default-37a8.jupyterlab-workspace
+docsrc/notebooks/.virtual_documents/lmr-real-pages2k-inpdt-verification.ipynb
+*.gz
+*.nc
+*.pkl

cfr/reconres.py

@@ -15,6 +15,9 @@
 import matplotlib.pyplot as plt
 from matplotlib import gridspec
 from .visual import CartopySettings
+from .reconjob import ReconJob
+import pandas as pd
+from . import utils,visual
 
 class ReconRes:
     ''' The class for reconstruction results '''
@@ -147,4 +150,106 @@ def plot_valid(self, recon_name_dict=None, target_name_dict=None,
                 fig[k], ax[k] = v.plot(label='recon', **valid_ts_kws)
             ax[k].set_ylabel(recon_name_dict[k])
 
-        return fig, ax
+        return fig, ax
+
+
+    def load_proxylabels(self, verbose=False):
+        """
+        Load proxy labels from the reconstruction results.
+        Proxy with "assim" means it is assimilated.
+        Proxy with "eval" means it is used for evaluation.
+        """
+        proxy_labels = []  # list of proxy labels
+        for path in self.paths:  # loop over all ensemble members
+            with xr.open_dataset(path) as ds_tmp:
+                proxy_labels.append(ds_tmp.attrs)  # dict for proxy labels
+
+        self.proxy_labels = proxy_labels
+        if verbose:
+            p_success(f">>> ReconRes.proxy_labels created")
+
+    def independent_verify(self, job_path, verbose=False, calib_period=(1850, 2000),min_verify_len=10):
+        """
+        Perform independent verification.
+        job_path (str): the path to the job.
+        verbose (bool, optional): print verbose information. Defaults to False.
+        """
+        # load the reconstructions for the "prior"
+        job = ReconJob()
+        job.load(job_path)
+        independent_info_list = []
+        for path_index ,path in enumerate(self.paths):
+            proxy_labels = self.proxy_labels[path_index]
+            job.load_clim(
+                tag="prior",
+                path_dict={
+                    "tas": path,
+                },
+                anom_period=(1951, 1980),
+            )
+            job.forward_psms(verbose=verbose)
+            if verbose:
+                p_success(f">>> Prior loaded from {path}")
+            # compare the pesudo-proxy records with the real records
+            calib_PDB = job.proxydb.filter(by="tag", keys=["calibrated"])
+            for i, (pname, proxy) in enumerate(calib_PDB.records.items()):
+                detail = proxy.psm.calib_details
+                attr_dict = {}
+                attr_dict['name'] = pname
+                attr_dict['seasonality'] = detail['seasonality']
+                if pname in proxy_labels['pids_assim']:
+                    attr_dict['assim'] = True
+                elif pname in proxy_labels['pids_eval']:
+                    attr_dict['assim'] = False
+                else:
+                    raise ValueError(f"Proxy {pname} is not labeled as assim or eval. Please check the proxy labels.")
+                reconstructed = pd.DataFrame(
+                    {
+                        "time": proxy.pseudo.time,
+                        "estimated": proxy.pseudo.value,
+                    }
+                )
+                real = pd.DataFrame(
+                    {
+                        "time": proxy.time,
+                        "observed": proxy.value,
+                    }
+                )
+                Df = real.dropna().merge(reconstructed, on="time", how="inner")
+                Df.set_index("time", drop=True, inplace=True)
+                Df.sort_index(inplace=True)
+                Df.astype(float)
+                masks = {
+                    "all": None,
+                    "in": (Df.index >= calib_period[0]) & (Df.index <= calib_period[1]), # in the calibration period
+                    "before": (Df.index < calib_period[0]), # before the calibration period
+                }
+                for mask_name, mask in masks.items():
+                    if mask is not None:
+                        Df_masked = Df[mask]
+                    else:
+                        Df_masked = Df
+                    if len(Df_masked) < min_verify_len:
+                        corr = np.nan
+                        ce = np.nan
+                    else:
+                        corr = Df_masked.corr().iloc[0, 1]
+                        ce = utils.coefficient_efficiency(
+                            Df_masked.observed.values, Df_masked.estimated.values
+                        )
+                    attr_dict[mask_name + '_corr'] = corr
+                    attr_dict[mask_name + '_ce'] = ce
+                independent_info_list.append(attr_dict)
+        independent_info_list = pd.DataFrame(independent_info_list)
+        self.independent_info_list = independent_info_list
+        if verbose:
+            p_success(f">>> Independent verification completed, results stored in ReconRes.independent_info_list")
+            p_success(f">>> Records Number: {len(independent_info_list)}")
+        return independent_info_list
+
+    def plot_independent_verify(self):
+        """
+        Plot the independent verification results.
+        """
+        fig, axs = visual.plot_independent_distribution(self.independent_info_list)
+        return fig, axs

cfr/visual.py

@@ -846,6 +846,37 @@ def in_notebook():
     return True
 
 
+def plot_independent_distribution(independent_info_list: pd.DataFrame,calib_period = [1880, 2000]):
+    """
+    Plot the distribution of independent test results 
+    """
+    fig = plt.figure(figsize=[20, 10])
+    gs = gridspec.GridSpec(2, 2)
+    gs.update(wspace=0.2, hspace=0.2)
+    bins = np.linspace(-1, 1, 41)
+    axs = {}
+    fs = 20
+    for ind_y, metric in enumerate(['corr','ce']):
+        i = 0
+        for label in [True, False]:
+            axs[str(label) + metric] = fig.add_subplot(gs[ind_y, i])
+            ax = axs[str(label) + metric]
+            table_use = independent_info_list[independent_info_list['assim'] == label]
+            ax.hist(
+                table_use[f'in_{metric}'], bins=bins, alpha=0.5, label=f'{calib_period[0]} to {calib_period[1]}', color='blue', density=True
+            )
+            ax.hist(table_use[f'before_{metric}'], bins=bins, alpha=0.5,
+                    label=f'before {calib_period[0]}', density=True, color='red')
+            ax.legend(loc='upper left', fontsize=fs)
+            ax.axvline(x=0, color='black', linestyle='--')
+            title = ['Assimilated Proxies', 'Non-assimilated Proxies'][i]
+            ax.set_title(title, fontsize=fs) if ind_y == 0 else None
+            ax.set_xlabel('Correlation', fontsize=fs) if metric == 'corr' else ax.set_xlabel('Correlation Efficiency', fontsize=fs)
+            ax.set_ylabel('Density', fontsize=fs)
+            i += 1
+    return fig, axs
+
+
 def showfig(fig, close=True):
     '''Show the figure