Added 100 member test ensemble

analysis/analyze_scalar.py

@@ -504,12 +504,6 @@ def run_sensitivity_analysis(
         type=str,
         default="grounding_line_flux",
     )
-    parser.add_argument(
-        "--fudge_factor",
-        help="""Observational uncertainty multiplier. Default=3""",
-        type=float,
-        default=3.0,
-    )
     parser.add_argument(
         "--notebook",
         help="""Use when running in a notebook to display a nicer progress bar. Default=False.""",
@@ -574,7 +568,6 @@ def run_sensitivity_analysis(
     basin_files = options.FILES
     engine = options.engine
     filter_range = options.filter_range
-    fudge_factor = options.fudge_factor
     notebook = options.notebook
     parallel = options.parallel
     reference_date = options.reference_date
@@ -591,6 +584,8 @@ def run_sensitivity_analysis(
     params = list(params_short_dict.keys())
     obs_cmap = config["Plotting"]["obs_cmap"]
     sim_cmap = config["Plotting"]["sim_cmap"]
+    grace_fudge_factor = config["Importance Sampling"]["grace_fudge_factor"]
+    mankoff_fudge_factor = config["Importance Sampling"]["mankoff_fudge_factor"]
 
     result_dir = Path(options.result_dir)
     data_dir = result_dir / Path("posteriors")
@@ -615,7 +610,7 @@ def run_sensitivity_analysis(
         "hatch.linewidth": 0.25,
     }
 
-    plt.rcParams.update(rcparams)
+    mpl.rcParams.update(rcparams)
 
     simulated_ds = prepare_simulations(
         basin_files, config, reference_date, parallel=parallel, engine=engine
@@ -721,7 +716,7 @@ def run_sensitivity_analysis(
             obs_std_vars=obs_std_vars_mankoff,
             sim_vars=sim_vars_mankoff,
             filter_range=filter_range,
-            fudge_factor=fudge_factor,
+            fudge_factor=mankoff_fudge_factor,
             params=params,
         )
     )
@@ -734,6 +729,7 @@ def run_sensitivity_analysis(
             filter_var=filter_var,
             filter_range=filter_range,
             fig_dir=fig_dir,
+            fudge_factor=mankoff_fudge_factor,
             config=config,
         )
 
@@ -750,7 +746,7 @@ def run_sensitivity_analysis(
             obs_mean_vars=obs_mean_vars_grace,
             obs_std_vars=obs_std_vars_grace,
             sim_vars=sim_vars_grace,
-            fudge_factor=fudge_factor,
+            fudge_factor=grace_fudge_factor,
             filter_range=filter_range,
             params=params,
         )
@@ -763,6 +759,7 @@ def run_sensitivity_analysis(
             simulated_posterior_grace.sel({"filtered_by": filter_var})[sim_plot_vars],
             filter_var=filter_var,
             filter_range=filter_range,
+            fudge_factor=grace_fudge_factor,
             fig_dir=fig_dir,
             config=config,
         )
@@ -805,7 +802,7 @@ def run_sensitivity_analysis(
     if "frontal_melt.routing.parameter_a" in prior_posterior.columns:
         prior_posterior["frontal_melt.routing.parameter_a"] *= 10**4
     if "ocean.th.gamma_T" in prior_posterior.columns:
-        prior_posterior["ocean.th.gamma_T"] *= 10**5
+        prior_posterior["ocean.th.gamma_T"] *= 10**4
     if "calving.vonmises_calving.sigma_max" in prior_posterior.columns:
         prior_posterior["calving.vonmises_calving.sigma_max"] *= 10**-3
     prior_posterior_sorted = sort_columns(prior_posterior, params_sorted_list)

data/07_prepare_calving.py

@@ -22,7 +22,7 @@
 Seasonal calving.
 """
 
-from argparse import ArgumentParser
+
 from pathlib import Path
 from typing import Union
 
@@ -32,7 +32,9 @@
 import xarray as xr
 
 
-def smoothed_step_function(time, amplification_factor: float = 1.0):
+def smoothed_step_function(
+    time: pd.DatetimeIndex, amplification_factor: float = 1.0, step_year: int = 2000
+) -> np.ndarray:
     """
     Generate a smoothed step function based on the day of the year.
 
@@ -45,6 +47,8 @@ def smoothed_step_function(time, amplification_factor: float = 1.0):
         A pandas DatetimeIndex representing the time series.
     amplification_factor : float, optional
         The factor by which the function value is amplified during the stay-at-one period, by default 1.0.
+    step_year : int, optional
+        The year at which the amplification factor is applied, by default 2000.
 
     Returns
     -------
@@ -77,7 +81,11 @@ def smoothed_step_function(time, amplification_factor: float = 1.0):
             & (day_of_year > end_ramp_up)
             & (day_of_year < start_ramp_down)
         )
-        values[stay_at_one_mask] = amplification_factor
+        if year >= step_year:
+            values[stay_at_one_mask] = amplification_factor
+        else:
+            values[stay_at_one_mask] = 1
+
         # Ramp down from 1 to 0
         ramp_down_mask = (
             (years == year)
@@ -90,33 +98,10 @@ def smoothed_step_function(time, amplification_factor: float = 1.0):
     return values
 
 
-# set up the option parser
-parser = ArgumentParser()
-parser.add_argument("FILE", nargs="*")
-parser.add_argument(
-    "--year_high",
-    type=float,
-    help="Start when high values are applied.",
-    default=2001,
-)
-parser.add_argument(
-    "--calving_low",
-    type=float,
-    help="Low value.",
-    default=1.0,
-)
-parser.add_argument(
-    "--calving_high",
-    type=float,
-    help="High value.",
-    default=1.5,
-)
-
-options = parser.parse_args()
-amplification_factors = [-1.0, 1.0, 1.05, 1.10, 1.20, 1.50, 2.00]
+amplification_factors = [-1.0, 1.0, 1.05, 1.10, 1.20]
 start_year = 1900
 end_year = 2025
-year_high = 2000
+step_year = 2000
 
 time = xr.date_range(str(start_year), str(end_year), freq="D")
 time_centered = time[:-1] + (time[1:] - time[:-1]) / 2
@@ -133,7 +118,9 @@ def smoothed_step_function(time, amplification_factor: float = 1.0):
     print(f"Processing {filename}")
     # Apply the smoothed step function to the time coordinate
     ds = da.to_dataset().copy()
-    data = smoothed_step_function(time_centered)
+    data = smoothed_step_function(
+        time_centered, amplification_factor=amplification_factor, step_year=step_year
+    )
     if amplification_factor == -1.0:
         data *= 0
     ds["frac_calving_rate"].values = data

pism_ragis/data/ragis_config.toml

@@ -3,7 +3,7 @@
 'calving.vonmises_calving.sigma_max' = '$\sigma_{\mathrm{max}}$ (kPa)'
 'calving.rate_scaling.file' = 'Calving'
 'geometry.front_retreat.prescribed.file' = 'Retreat Method'
-'ocean.th.gamma_T' = '$\gamma$ (10$^{-5}$ 1)'
+'ocean.th.gamma_T' = '$\gamma$ (10$^{-4}$ 1)'
 'surface.given.file' = 'Climate Forcing'
 'ocean.th.file' = 'Ocean Forcing'
 'frontal_melt.routing.parameter_a' = '$a$ (10$^{-4}$ m$^{-\alpha}$ day$^{\alpha-1}$ Celsius$^{-\beta}$)'
@@ -101,5 +101,10 @@ gis = {D = "grounding_line_flux", MB = "mass_balance", SMB = "surface_mass_balan
 sim_alpha = 0.5
 sim_cmap = ["#CC6677", "#882255"]
 obs_alpha = 1.0
-obs_cmap = ["0.8", "0.7"]
+obs_cmap = ["0.8", "0.9"]
 hist_cmap = ["#a6cee3", "#1f78b4"]
+
+['Importance Sampling']
+
+grace_fudge_factor = 10
+mankoff_fudge_factor = 3

pism_ragis/plotting.py

@@ -136,6 +136,7 @@ def plot_basins(
     filter_var: str,
     filter_range: List[int] = [1990, 2019],
     fig_dir: Union[str, Path] = "figures",
+    fudge_factor: float = 3.0,
     plot_range: List[int] = [1980, 2020],
     config: Dict = {},
 ):
@@ -160,6 +161,9 @@ def plot_basins(
         A list containing the start and end years for filtering, by default [1990, 2019].
     fig_dir : Union[str, Path], optional
         The directory where figures will be saved, by default "figures".
+    fudge_factor : float, optional
+        A multiplicative factor applied to the observed standard deviation to widen the likelihood function,
+        allowing for greater tolerance in the matching process, by default 3.0.
     plot_range : List[int], optional
         A list containing the start and end years for plotting, by default [1980, 2020].
     config : Dict, optional
@@ -192,6 +196,7 @@ def plot_basins(
                 filter_var=filter_var,
                 filter_range=filter_range,
                 fig_dir=fig_dir,
+                fudge_factor=fudge_factor,
                 obs_alpha=obs_alpha,
                 sim_alpha=sim_alpha,
             )
@@ -268,6 +273,7 @@ def plot_obs_sims(
     filter_var: str,
     filter_range: List[int] = [1990, 2019],
     fig_dir: Union[str, Path] = "figures",
+    fudge_factor: float = 3.0,
     reference_year: float = 1986.0,
     sim_alpha: float = 0.4,
     obs_alpha: float = 1.0,
@@ -294,6 +300,9 @@ def plot_obs_sims(
         Range of years for filtering, by default [1990, 2019].
     fig_dir : Union[str, Path], optional
         Directory to save the figures, by default "figures".
+    fudge_factor : float, optional
+        A multiplicative factor applied to the observed standard deviation to widen the likelihood function,
+        allowing for greater tolerance in the matching process, by default 3.0.
     reference_year : float, optional
         Reference year for cumulative mass balance, by default 1986.0.
     sim_alpha : float, optional
@@ -334,6 +343,11 @@ def plot_obs_sims(
         "mass_flux_uncertainty"
     ]
 
+    if filter_var == "grounding_line_flux":
+        m = -1
+    else:
+        m = 1
+
     with mpl.rc_context({"font.size": fontsize}):
 
         fig, axs = plt.subplots(
@@ -379,6 +393,18 @@ def plot_obs_sims(
                 lw=0,
             )
 
+        if filter_var in obs.data_vars:
+            axs[m].fill_between(
+                obs["time"],
+                obs[filter_var] - fudge_factor * obs[filter_var + "_uncertainty"],
+                obs[filter_var] + fudge_factor * obs[filter_var + "_uncertainty"],
+                alpha=0.5,
+                edgecolor="k",
+                facecolor="w",
+                hatch="//",
+                lw=0.25,
+            )
+
         sim_cis = []
         if sim_prior is not None:
             sim_prior = sim_prior[
@@ -459,22 +485,22 @@ def plot_obs_sims(
 
         if sim_posterior is not None:
             y_min, y_max = axs[-1].get_ylim()
-            scaler = y_min + (y_max - y_min) * 0.05
+            scaler = y_min + (y_max - y_min) * 0.0505
             obs_filtered = obs.sel(
                 time=slice(f"{filter_range[0]}", f"{filter_range[-1]}")
             )
             filter_range_ds = obs_filtered[mass_cumulative_varname]
             filter_range_ds *= 0
             filter_range_ds += scaler
             _ = filter_range_ds.plot(
-                ax=axs[-1], lw=1, ls="solid", color="k", label="Filtering Range"
+                ax=axs[m], lw=1, ls="solid", color="k", label="Filtering Range"
             )
             x_s = (
                 filter_range_ds.time.values[0]
                 + (filter_range_ds.time.values[-1] - filter_range_ds.time.values[0]) / 2
             )
             y_s = scaler
-            axs[-1].text(
+            axs[m].text(
                 x_s,
                 y_s,
                 "Filtering Range",

uq/draw_ragis_samples.py

@@ -54,7 +54,7 @@
     "ragis": {
         "uq": {
             "calving.vonmises_calving.sigma_max": uniform(loc=350_000, scale=300_000),
-            "calving.rate_scaling.file": randint(0, 7),
+            "calving.rate_scaling.file": randint(0, 5),
             "ocean.th.gamma_T": uniform(loc=0.75e-4, scale=0.75e-4),
             "ocean_file": randint(0, len(gcms)),
             "climate_file": randint(0, len(rcms)),
@@ -316,7 +316,7 @@ def convert_samples(unif_sample):
 
 
 dist_sample, n_samples = convert_samples(unif_sample)
-dist_median_sample = convert_samples(np.median(unif_sample, axis=0, keepdims=True))
+dist_median_sample, _ = convert_samples(np.median(unif_sample, axis=0, keepdims=True))
 
 dist_median_sample[0, 1] = "seasonal_calving_id_2_1900_2025.nc"
 dist_median_sample[0, 3] = gcms[0]
@@ -356,7 +356,7 @@ def convert_samples(unif_sample):
 
 median_outfile = outfile.parent / Path(f"median_{outfile.name}")
 median_df = pd.DataFrame(
-    dist_median_sample, columns=keys, index=["MEDIAN-FREE", "MEDIAN-PRESCRIBED"]
+    dist_median_sample, columns=keys_prior, index=["MEDIAN-FREE", "MEDIAN-PRESCRIBED"]
 )
 median_df.to_csv(median_outfile, index=True, index_label="id")