feat: added scatter summary plot per horizon for all indicators (#753)

* feat: added scatter summary plot per horizon for all indicators

* refactored main

* fix: b2a societal cost add, with colours

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* release notes

* fix: add option for not collection years

* fix: remove b2a non-central, add stats values box

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix: wildcard was only working if scenarios are enabled

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: AndreasHD11

doc/release_notes.md

@@ -31,6 +31,8 @@
 
 * Change dispatch of biomass and biogas generators in CBA by (a) setting the buses' marginal prices as the generators' marginal costs and (b) removing energy budget constraints ([#719](https://github.com/open-energy-transition/open-tyndp/pull/719)).
 
+* Add CBA per horizon summary plots for each indicator benchmarking TYNDP and Open-TYNDP ([#753](https://github.com/open-energy-transition/open-tyndp/pull/753)).
+
 **Bugfixes and Compatibility**
 
 * Rename bus for `t339` project (Tyrrhenian) from ITSI to ITVI ([#751](https://github.com/open-energy-transition/open-tyndp/pull/751)).

rules/cba.smk

@@ -556,6 +556,36 @@ rule summarize_indicators_per_project:
         "../scripts/cba/summarize_indicators.py"
 
 
+def summary_benchmark_indicators(w):
+    """
+    Returns Indicator CSVs as inputs for the per-horizon summary benchmark plot.
+    If collection scenarios, returns the weighted-average ensemble indicators CSV as inputs for plotting.
+    """
+    run = get_run_name(w)
+    if run in cba_collection_scenarios(w):
+        return expand(
+            rules.average_indicators_per_project_and_planning_horizon.output.indicators,
+            planning_horizons=[w.planning_horizons],
+            cba_project=cba_projects(w),
+            run=[run],
+        )
+    return expand(
+        rules.collect_indicators.output.indicators,
+        planning_horizons=[w.planning_horizons],
+        run=[run],
+    )
+
+
+rule plot_summary_projects_benchmark:
+    input:
+        indicators=summary_benchmark_indicators,
+    output:
+        plot_file=RESULTS
+        + "cba/ensemble_plots/summary_benchmark_{planning_horizons}.png",
+    script:
+        "../scripts/cba/plot_benchmark_indicators.py"
+
+
 rule summarize_all_indicators:
     input:
         indicators=lambda w: expand(
@@ -691,6 +721,13 @@ def collect_cba_scenario_inputs(w):
             run=cba_scenarios(w),
         )
     )
+    inputs.extend(
+        expand(
+            rules.plot_summary_projects_benchmark.output.plot_file,
+            planning_horizons=config_provider("cba", "planning_horizons")(w),
+            run=cba_scenarios(w),
+        )
+    )
 
     run = get_run_name(w)
     if run in cba_collection_scenarios(w):
@@ -743,6 +780,13 @@ def cba_ensemble_inputs(w):
             run=runs,
         )
     )
+    inputs.extend(
+        expand(
+            rules.plot_summary_projects_benchmark.output.plot_file,
+            planning_horizons=config["cba"]["planning_horizons"],
+            run=runs,
+        )
+    )
     return inputs
 
 

scripts/cba/plot_benchmark_indicators.py

@@ -352,6 +352,152 @@ def plot_project_benchmarks(
     plt.close(fig)
 
 
+def plot_summary_projects_benchmark(
+    df: pd.DataFrame,
+    output_path: Path,
+    planning_horizon: str | None = None,
+    area_subtitle: str | None = None,
+    eps: float = 1e-6,
+) -> None:
+    """Plot one summary subplot per indicator for each horizon"""
+    model_df = df[df["source"] == "Open-TYNDP"]
+    if model_df.empty:
+        logger.info("No Open-TYNDP data to plot")
+        return
+    weighted_df = model_df[model_df["cyear"] == "weighted-average"]
+    if not weighted_df.empty:
+        model_df = weighted_df
+
+    average_df = pd.concat(
+        [model_df, df[df["source"] == "TYNDP 2024"]], ignore_index=True
+    )
+
+    indicators = sorted(
+        average_df.loc[average_df["source"] == "Open-TYNDP", "indicator"]
+        .dropna()
+        .unique()
+    )
+    project_ids = sorted(
+        average_df.loc[average_df["source"] == "Open-TYNDP", "project_id"]
+        .dropna()
+        .unique()
+    )
+
+    plot_items = {}
+    for indicator in indicators:
+        pairs = []
+        for project_id in project_ids:
+            project_df = average_df[average_df["project_id"] == project_id]
+            if indicator == "B2a_societal_cost_variation":
+                model_val = select_value_by_subindex(
+                    project_df, indicator, "Open-TYNDP", "central"
+                )
+                bench_val = select_value_by_subindex(
+                    project_df, indicator, "TYNDP 2024", "central"
+                )
+                if model_val is not None and bench_val is not None:
+                    pairs.append((bench_val, model_val))
+            else:
+                model = benchmark_range(project_df, indicator, source="Open-TYNDP")
+                bench = benchmark_range(project_df, indicator, source="TYNDP 2024")
+                if model is not None and bench is not None:
+                    pairs.append((bench[1], model[1]))
+        if pairs:
+            plot_items[indicator] = pairs
+    if not plot_items:
+        logger.info("Incomplete benchmark data to plot")
+        return
+
+    ncols = min(4, len(plot_items))
+    nrows = (len(plot_items) + ncols - 1) // ncols
+    fig, axes = plt.subplots(
+        nrows=nrows,
+        ncols=ncols,
+        figsize=(3.6 * ncols, 3.3 * nrows),
+        squeeze=False,
+    )
+
+    for ax, (indicator, pairs) in zip(axes.flatten(), plot_items.items()):
+        xs = [p[0] for p in pairs]
+        ys = [p[1] for p in pairs]
+        colors = "tab:blue"
+        ax.scatter(
+            xs,
+            ys,
+            s=20,
+            color=colors,
+            alpha=0.6,
+            edgecolor="white",
+            linewidth=0.5,
+        )
+        ax.axline(
+            (0, 0), slope=1, color="black", linestyle="--", linewidth=1, alpha=0.5
+        )
+
+        combined = [abs(v) for v in (*xs, *ys) if v != 0]
+        if combined:
+            abs_max = max(combined)
+            abs_min = min(combined)
+            if abs_min > 0 and abs_max / abs_min >= 1e3:
+                linthresh = max(abs_min, 1.0)
+                ax.set_xscale("symlog", linthresh=linthresh)
+                ax.set_yscale("symlog", linthresh=linthresh)
+
+        units = average_df.loc[
+            (average_df["indicator"] == indicator)
+            & (average_df["source"] == "Open-TYNDP"),
+            "units",
+        ].dropna()
+        unit_label = units.iloc[0] if not units.empty else ""
+        ax.set_title(
+            f"{indicator} ({unit_label})" if unit_label else indicator, fontsize=9
+        )
+        ax.set_xlabel("TYNDP 2024")
+        ax.set_ylabel("Open-TYNDP")
+        ax.axhline(0, color="gray", linewidth=0.5, alpha=0.4)
+        ax.axvline(0, color="gray", linewidth=0.5, alpha=0.4)
+        ax.grid(alpha=0.3)
+
+        benchmark_df = pd.DataFrame(pairs, columns=["TYNDP 2024", "Open-TYNDP"])
+        errors = (
+            (benchmark_df["Open-TYNDP"] - benchmark_df["TYNDP 2024"]).abs()
+            / (
+                (benchmark_df["Open-TYNDP"].abs() + benchmark_df["TYNDP 2024"].abs())
+                / 2
+                + eps
+            )
+            * 100
+        )
+        values = (
+            f"n = {len(pairs)}\n"
+            f"sMAPE = {errors.mean():.1f}%\n"
+            f"sMdAPE = {errors.median():.1f}%"
+        )
+        ax.text(
+            0.05,
+            0.95,
+            values,
+            transform=ax.transAxes,
+            fontsize=6,
+            verticalalignment="top",
+            bbox=dict(boxstyle="round", facecolor="white", alpha=0.8),
+        )
+
+    for ax in axes.flatten()[len(plot_items) :]:
+        ax.axis("off")
+
+    title = "Benchmark of indicators across all projects"
+    if planning_horizon:
+        title += f" ({planning_horizon})"
+    fig.suptitle(title, y=0.995)
+    if area_subtitle:
+        fig.text(0.5, 0.965, area_subtitle, ha="center", va="center", fontsize=9)
+
+    fig.tight_layout(rect=[0, 0, 1, 0.94])
+    fig.savefig(output_path, dpi=400)
+    plt.close(fig)
+
+
 def create_plots(
     indicators_file: str | Path,
     output_path: str | Path,
@@ -435,6 +581,20 @@ def create_plots(
     set_scenario_config(snakemake)
 
     planning_horizon = snakemake.wildcards.get("planning_horizons")
-    output_target = snakemake.output.get("plot_file") or snakemake.output.plot_dir
     area = snakemake.config.get("cba", {}).get("area")
-    create_plots(snakemake.input.indicators, output_target, planning_horizon, area)
+
+    if "cba_project" in snakemake.wildcards.keys():
+        output_target = snakemake.output.get("plot_file") or snakemake.output.plot_dir
+        create_plots(snakemake.input.indicators, output_target, planning_horizon, area)
+    elif not snakemake.input.indicators:
+        logger.warning(
+            "No indicators input files for summary plot: %s", snakemake.output.plot_file
+        )
+    else:
+        df = pd.concat(map(pd.read_csv, snakemake.input.indicators), ignore_index=True)
+        plot_summary_projects_benchmark(
+            df,
+            Path(snakemake.output.plot_file),
+            planning_horizon,
+            format_area_subtitle(area),
+        )