diff --git a/ema_workbench/analysis/regional_sa.py b/ema_workbench/analysis/regional_sa.py index e4d453ade..c17087a86 100644 --- a/ema_workbench/analysis/regional_sa.py +++ b/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() diff --git a/ema_workbench/analysis/scenario_discovery_util.py b/ema_workbench/analysis/scenario_discovery_util.py index a8f743f6c..4c362d26a 100644 --- a/ema_workbench/analysis/scenario_discovery_util.py +++ b/ema_workbench/analysis/scenario_discovery_util.py @@ -320,7 +320,7 @@ def plot_pair_wise_scatter( restricted_dims, diag="kde", upper="scatter", - lower="contour", + lower="hist", fill_subplots=True, ): """helper function for pair wise scatter plotting @@ -399,7 +399,9 @@ def plot_pair_wise_scatter( elif upper == "none": None else: - raise NotImplementedError(f"upper = {upper} not implemented.") + raise NotImplementedError( + f"upper = {upper} not implemented. Use either 'scatter', 'contour', 'hist' (bivariate histogram) or None plots for upper triangle." + ) # lower triangle if lower == "contour": @@ -417,7 +419,9 @@ def plot_pair_wise_scatter( elif lower == "none": raise ValueError(f"Lower triangle cannot be none.") else: - raise NotImplementedError(f"lower = {lower} not implemented.") + raise NotImplementedError( + f"lower = {lower} not implemented. Use either 'scatter', 'contour' or 'hist' (bivariate histogram) plots for lower triangle." + ) # diagonal if diag == "cdf": @@ -425,34 +429,35 @@ def plot_pair_wise_scatter( elif diag == "kde": grid.map_diag(sns.kdeplot, fill=False, common_norm=False, cut=0) else: - raise NotImplementedError(f"diag = {diag} not implemented.") + raise NotImplementedError( + f"diag = {diag} not implemented. Use either 'kde' (kernel density estimate) or 'cdf' (cumulative density function)." + ) # draw box pad = 0.1 cats = set(categorical_columns) for row, ylabel in zip(grid.axes, grid.y_vars): - ylim = boxlim[ylabel] - - if ylabel in cats: - height = (len(ylim[0]) - 1) + pad - y = -pad / 2 - # y = -0.2 - # height = len(ylim[0]) - 0.6 # 2 * 0.2 - else: - y = ylim[0] - height = ylim[1] - ylim[0] - for ax, xlabel in zip(row, grid.x_vars): if ylabel == xlabel: continue + xrange = ax.get_xlim()[1] - ax.get_xlim()[0] + yrange = ax.get_ylim()[1] - ax.get_ylim()[0] + + ylim = boxlim[ylabel] + + if ylabel in cats: + height = (len(ylim[0]) - 1) + pad * yrange + y = -yrange * pad / 2 + else: + y = ylim[0] + height = ylim[1] - ylim[0] + if xlabel in cats: xlim = boxlim.at[0, xlabel] - width = (len(xlim) - 1) + pad - x = -pad / 2 - # x = -0.2 - # width = len(xlim) - 0.6 # 2 * 0.2 + width = (len(xlim) - 1) + pad * xrange + x = -xrange * pad / 2 else: xlim = boxlim[xlabel] x = xlim[0] @@ -462,7 +467,7 @@ def plot_pair_wise_scatter( box = patches.Rectangle( xy, width, height, edgecolor="red", facecolor="none", lw=3, zorder=100 ) - if ax.has_data() == True: + if ax.has_data(): # keeps box from being drawn in upper triangle if empty ax.add_patch(box) else: ax.set_axis_off() @@ -508,7 +513,7 @@ def plot_pair_wise_scatter( upper = data[subplot.get_xlabel()].max() lower = data[subplot.get_xlabel()].min() - pad_rel = (upper - lower) * pad # padding relative to range of data points + pad_rel = (upper - lower) * 0.1 # padding relative to range of data points subplot.set_xlim(lower - pad_rel, upper + pad_rel) @@ -516,7 +521,7 @@ def plot_pair_wise_scatter( upper = data[subplot.get_ylabel()].max() lower = data[subplot.get_ylabel()].min() - pad_rel = (upper - lower) * pad # padding relative to range of data points + pad_rel = (upper - lower) * 0.1 # padding relative to range of data points subplot.set_ylim(lower - pad_rel, upper + pad_rel) diff --git a/ema_workbench/em_framework/points.py b/ema_workbench/em_framework/points.py index 88fc630be..317c2f866 100644 --- a/ema_workbench/em_framework/points.py +++ b/ema_workbench/em_framework/points.py @@ -108,6 +108,16 @@ def __init__(self, name, model_name, policy, scenario, experiment_id): self.model_name = model_name self.scenario = scenario + def __repr__(self): + return ( + f"Experiment(name={self.name!r}, model_name={self.model_name!r}, " + f"policy={self.policy!r}, scenario={self.scenario!r}, " + f"experiment_id={self.experiment_id!r})" + ) + + def __str__(self): + return f"Experiment {self.experiment_id} (model: {self.model_name}, policy: {self.policy.name}, scenario: {self.scenario.name})" + class ExperimentReplication(NamedDict): """helper class that combines scenario, policy, any constants, and diff --git a/ema_workbench/examples/example_lake_model.py b/ema_workbench/examples/example_lake_model.py index 87a5a3cb4..f0c013df8 100644 --- a/ema_workbench/examples/example_lake_model.py +++ b/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 diff --git a/ema_workbench/examples/regional_sa_flu.py b/ema_workbench/examples/regional_sa_flu.py new file mode 100644 index 000000000..fd0e4e0eb --- /dev/null +++ b/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()