Merge branch 'master' into MPIEvaluator

.pre-commit-config.yaml

@@ -4,16 +4,16 @@ ci:
 
 repos:
 -   repo: https://github.com/psf/black
-    rev: 23.9.1
+    rev: 23.10.1
     hooks:
       - id: black-jupyter
 -   repo: https://github.com/asottile/pyupgrade
-    rev: v3.13.0
+    rev: v3.15.0
     hooks:
     -   id: pyupgrade
         args: [--py39-plus]
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0  # Use the ref you want to point at
+    rev: v4.5.0  # Use the ref you want to point at
     hooks:
     -   id: trailing-whitespace
     -   id: check-toml

ema_workbench/analysis/regional_sa.py

@@ -82,45 +82,15 @@ def plot_discrete_cdf(ax, unc, x, y, xticklabels_on, ccdf):
             x_plot = [j * 1, j * 1 + 1]
             y_plot = [freq] * 2
 
-            ax.plot(x_plot, y_plot, c=cp[i + 1], label=i == 1)
-            ax.scatter(
-                x_plot[0],
-                y_plot[0],
-                edgecolors=cp[i + 1],
-                facecolors=cp[i + 1],
-                linewidths=1,
-                zorder=2,
-            )
-            ax.scatter(
-                x_plot[1],
-                y_plot[0],
-                edgecolors=cp[i + 1],
-                facecolors="white",
-                linewidths=1,
-                zorder=2,
-            )
+            ax.plot(x_plot, y_plot, c=cp[i + 1], label=i == 1, marker="o")
 
             # misnomer
             cum_freq_un = (j + 1) / n_cat
             if ccdf:
                 cum_freq_un = (len(freqs) - j - 1) / n_cat
 
-            ax.plot(x_plot, [cum_freq_un] * 2, lw=1, c="darkgrey", zorder=1, label="x==y")
-            ax.scatter(
-                x_plot[0],
-                cum_freq_un,
-                edgecolors="darkgrey",
-                facecolors="darkgrey",
-                linewidths=1,
-                zorder=1,
-            )
-            ax.scatter(
-                x_plot[1],
-                cum_freq_un,
-                edgecolors="darkgrey",
-                facecolors="white",
-                linewidths=1,
-                zorder=1,
+            ax.plot(
+                x_plot, [cum_freq_un] * 2, lw=1, c="darkgrey", zorder=1, label="x==y", marker="o"
             )
 
     ax.set_xticklabels([])
@@ -256,10 +226,18 @@ def plot_cdfs(x, y, ccdf=False):
     x = x.copy()
 
     try:
-        x = x.drop("scenario", axis=1)
+        x = x.drop("scenario_id", axis=1)
     except KeyError:
         pass
 
+    for entry in ["model", "policy"]:
+        if x.loc[:, entry].unique().shape != (1,):
+            continue
+        try:
+            x = x.drop(entry, axis=1)
+        except KeyError:
+            pass
+
     uncs = x.columns.tolist()
     cp = sns.color_palette()
 

ema_workbench/examples/example_lake_model.py

@@ -35,37 +35,47 @@ def lake_problem(
         decisions = [kwargs[str(i)] for i in range(100)]
     except KeyError:
         decisions = [0] * 100
+        print("No valid decisions found, using 0 water release every year as default")
 
-    Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
     nvars = len(decisions)
-    X = np.zeros((nvars,))
-    average_daily_P = np.zeros((nvars,))
     decisions = np.array(decisions)
-    reliability = 0.0
 
-    for _ in range(nsamples):
-        X[0] = 0.0
+    # Calculate the critical pollution level (Pcrit)
+    Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
 
-        natural_inflows = np.random.lognormal(
-            math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
-            math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
-            size=nvars,
+    # Generate natural inflows using lognormal distribution
+    natural_inflows = np.random.lognormal(
+        mean=math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
+        sigma=math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
+        size=(nsamples, nvars),
+    )
+
+    # Initialize the pollution level matrix X
+    X = np.zeros((nsamples, nvars))
+
+    # Loop through time to compute the pollution levels
+    for t in range(1, nvars):
+        X[:, t] = (
+            (1 - b) * X[:, t - 1]
+            + (X[:, t - 1] ** q / (1 + X[:, t - 1] ** q))
+            + decisions[t - 1]
+            + natural_inflows[:, t - 1]
         )
 
-        for t in range(1, nvars):
-            X[t] = (
-                (1 - b) * X[t - 1]
-                + X[t - 1] ** q / (1 + X[t - 1] ** q)
-                + decisions[t - 1]
-                + natural_inflows[t - 1]
-            )
-            average_daily_P[t] += X[t] / float(nsamples)
+    # Calculate the average daily pollution for each time step
+    average_daily_P = np.mean(X, axis=0)
 
-        reliability += np.sum(X < Pcrit) / float(nsamples * nvars)
+    # Calculate the reliability (probability of the pollution level being below Pcrit)
+    reliability = np.sum(X < Pcrit) / float(nsamples * nvars)
 
+    # Calculate the maximum pollution level (max_P)
     max_P = np.max(average_daily_P)
+
+    # Calculate the utility by discounting the decisions using the discount factor (delta)
     utility = np.sum(alpha * decisions * np.power(delta, np.arange(nvars)))
-    inertia = np.sum(np.absolute(np.diff(decisions)) < 0.02) / float(nvars - 1)
+
+    # Calculate the inertia (the fraction of time steps with changes larger than 0.02)
+    inertia = np.sum(np.abs(np.diff(decisions)) > 0.02) / float(nvars - 1)
 
     return max_P, utility, inertia, reliability
 

ema_workbench/examples/regional_sa_flu.py

@@ -0,0 +1,18 @@
+""" A simple example of performing regional sensitivity analysis
+
+
+"""
+import matplotlib.pyplot as plt
+
+from ema_workbench.analysis import regional_sa
+from ema_workbench import ema_logging, load_results
+
+fn = "./data/1000 flu cases with policies.tar.gz"
+results = load_results(fn)
+x, outcomes = results
+
+y = outcomes["deceased population region 1"][:, -1] > 1000000
+
+fig = regional_sa.plot_cdfs(x, y)
+
+plt.show()