analysis_script.py

#!/usr/bin/python3.10
########################################################################################
# analysis_script.py - Analysis script for stranded-assets outputs.		   			   #
# 																					   #
# Author: Ben Winchester															   #
# Copyright: Ben Winchester, 2023		   											   #
# Date created: 21/09/2023															   #
# License: MIT																		   #
########################################################################################
"""
analysis_script.py - Analysis script for stranded-assets outputs.

Several optimisations are run using CLOVER [1] to produce results for the stranded-
assests paper. The results from these optimisations need to be plotted in a sensible
way.

1. Winchester B., Beath H., Nelson J., and Sandwell P., CLOVER: Continuous Lifetime
	Optimisation of Variable Electricity Resources. Available from:
	https://github.com/CLOVER-energy/CLOVER. [Accessed on 20/09/2023]

"""

import argparse
import collections
import itertools
import matplotlib.colors as mcolors
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import os
import pdb
import sys

from matplotlib import rc
from typing import Any

import json
import numpy as np
import pandas as pd
import seaborn as sns

from tqdm import tqdm

# DPI:
#   The DPI to use when saving figures.
DPI: int = 400

# INDEX:
#   Variable used for plotting to ensure non-overwriting behaviour when saving plots
INDEX: int = 14

# MONTHS:
#   A list of names of months
MONTHS: list[str] = [
    "Jan",
    "Feb",
    "Mar",
    "Apr",
    "May",
    "Jun",
    "Jul",
    "Aug",
    "Sep",
    "Oct",
    "Nov",
    "Dec",
]

PURPLE_GREEN_COLOURS: list[str] = [
    "#003D30",
    "#005745",
    "#00735C",
    "#009175",
    "#00AD8E",
    "#00CBA7",
    "#00EBC1",
    "#86FFDE",
    "#FFCCFE",
    "#FF92FD",
    "#FF3CFE",
    "#DA00FD",
    "#A700FC",
    "#8400CD",
    "#8400CD",
    "#65019F",
    "#450270",
]

# Plitting context
rc("font", **{"family": "sans-serif", "sans-serif": ["Arial"]})
sns.set_context("notebook")
sns.set_style("ticks")

plt.rcParams["pdf.fonttype"] = 42
plt.rcParams["ps.fonttype"] = 42


def _parse_args(args: list[Any]) -> argparse.Namespace:
    """
    Parse the command-line arguments.

    :param: args
            Unparsed command-line arguments.

    """

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "--grid-output-file",
        "-grid",
        "-g",
        help="The name of the output file containing combined results.",
        type=str,
    )

    parser.add_argument(
        "--unsubsidised-output-file",
        "-unsubf",
        help="The name of the output file containing UN-subsidised combined results.",
        type=str,
    )

    parser.add_argument(
        "--stranded-assets-output-file",
        "-strandf",
        help="The name of the stranded-assets exploratory output file.",
        type=str,
    )

    parser.add_argument(
        "--subsidised-output-file",
        "-subf",
        help="The name of the output file containing subsidised combined results.",
        type=str,
    )

    return parser.parse_args(args)


def _heatmap_plot(
    data: pd.DataFrame,
    variable_name: str,
    *,
    annot: bool = False,
    annot_kws: dict[str, Any] | None = None,
    ax: plt.Axes | None = None,
    center: float | None = None,
    cmap=None,
    cmap_colour: str | None = None,
    cmap_name: str | None = None,
    column_name: str = "grid_hours",
    contour: bool = False,
    default_value: float = 0,
    figsize: tuple[float, float] = (12, 8),
    fmt: str = ".2g",
    index_name: str = "service_hours",
    prioritise_self_generation: bool = True,
    vmax: float | None = None,
    vmin: float | None = None,
    ylabel: str | None = None,
) -> plt.Axes:
    """
    Plot the PV size heatmap

    :param: data
            The data to plot.

    :param: variable_name
        The name of the variable to plot.

    :param: annot
        Whether to annotate.

    :param: annot_kws
        Keyword arguments to pass through to the plotting function.

    :param: ax
        If specified, axes on which to plot the data.

    :param: center
        The midpoint for the data to plot.

    :param: cmap
        The colourmap to use.

    :param: cmap_name
        The name of the colourmap to use.

    :param: column_name
        The variable to use for the y-axis of the data plots.

    :param: contour
        Whether to overlay contours on the plot (True) or use heatmaps only (False).

    :param: default_value
        The value to use for data which has not been provided in the dataset.

    :param: figsize
        The size, in inches, to use for the figure, of the format (x, y) where x and y
        are both floats.

    :param: fmt
        The string formatting to use for the annotations, if specified, on the heatmaps.

    :param: index_name
        The variable name to use for the x-axis of the data plots.

    :param: prioritise_self_generation
        Whether to show results for self- or grid-prioritisation.

    :param: vmax
        The maximum value for the heatmap to show.

    :param: vmin
        The minimum value for the heatmap to show.

    :param: ylabel
        The label to use for the y-axis, if specified.

    """

    # Plot only for prioritising self-generation
    relevant_data = data[
        data["prioritise_self_generation"] == prioritise_self_generation
    ]

    plotting_data = relevant_data.pivot(
        index=index_name, columns=column_name, values=variable_name
    )  # .fillna(default_value)

    if cmap_colour is not None:
        cmap = sns.light_palette(cmap_colour, as_cmap=True)
        dark_cmap = sns.dark_palette(cmap_colour, reverse=True, as_cmap=True)
    elif cmap is not None:
        dark_cmap = None
    else:
        cmap = None
        dark_cmap = None

    if ax is None:
        fig, ax = plt.subplots(figsize=figsize)

    sns.heatmap(
        plotting_data,
        alpha=0 if annot else 1,
        annot=annot,
        annot_kws=(
            annot_kws.update({"size": 10})
            if (annot_kws is not None and annot)
            else {"size": 10} if annot else None
        ),
        ax=ax,
        cbar=not annot,
        cbar_kws=(
            None
            if annot
            else {
                "label": (
                    ylabel
                    if ylabel is not None
                    else variable_name.replace("_", " ").capitalize()
                )
            }
        ),
        center=center,
        cmap=cmap if cmap_name is None else cmap_name,
        fmt=fmt,
        vmin=vmin,
        vmax=vmax,
        square=True,
        linewidths=1,
    )

    if contour:
        ax.contour(
            plotting_data,
            cmap=(
                (
                    f"{cmap_name}_r"
                    if not cmap_name.endswith("_r")
                    else cmap_name.replace("_r", "")
                )
                if cmap_name is not None
                else dark_cmap
            ),
            vmin=vmin,
            vmax=vmax,
        )

    ax.set_xlabel(column_name.replace("_", " ").capitalize())
    ax.set_ylabel(index_name.replace("_", " ").capitalize())
    ax.invert_yaxis()

    return ax


def _plot_solar_heatmap() -> None:
    """
    Plot the solar generation heatmap.

    Plot the four-subplot plot with the solar generation, grid availabilitiy,
    community energy demand and enterprise/commercial energy demand

    """

    sns.set_style("ticks")

    # sns.set_context("poster")

    def _plot_twenty_year_heatmap(
        data_to_plot: pd.DataFrame,
        *,
        axis: plt.Axes | None = None,
        cmap_colour: str | None = None,
        cmap_name: str = "PuBu",
        label: str,
        title: str,
        vmax: float | None = None,
        vmin: float | None = None,
    ) -> None:
        """
        Plot the twenty-year heatmap.

        :param: data
            The data to plot

        :param: label
            The label to use.

        """

        # Calculate the yearly average performance
        _days_so_far: int = 0
        _monthly_data: dict[int, np.ndarray] = {}
        _yearly_average_data = np.reshape(
            np.reshape(data_to_plot, (20, 8760)).mean(axis=0), (365, 24)
        )

        for index, month_length in enumerate(
            [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
        ):
            _monthly_data[MONTHS[index]] = np.reshape(
                _yearly_average_data[_days_so_far : _days_so_far + month_length],
                (month_length, 24),
            ).mean(axis=0)
            _days_so_far += month_length

        # Generate the colour map
        if cmap_colour is not None:
            cmap = sns.light_palette("#E04606", as_cmap=True)
        else:
            cmap = cmap_name

        if axis is None:
            axis = sns.heatmap(
                pd.DataFrame(_monthly_data).transpose(),
                cmap=cmap,
                linewidths=1,
                # square=True,
                vmin=vmin,
                vmax=vmax,
                cbar_kws={"label": label},
            )

        else:
            sns.heatmap(
                pd.DataFrame(_monthly_data).transpose(),
                ax=axis,
                cmap=cmap,
                linewidths=1,
                square=True,
                vmin=vmin,
                vmax=vmax,
                cbar_kws={"label": label},
            )

        axis.set_xlabel("Hour of the day")
        axis.set_ylabel("Month")
        axis.set_title(title, fontweight="bold")

    # Create the subplot space
    fig, axes = plt.subplots(2, 2, figsize=(64 / 5, 32 / 5))
    fig.subplots_adjust(wspace=0.25)

    # Plot the solar heatmap.
    try:
        with open(
            os.path.join(
                os.path.expanduser("~"),
                "clover_locations",
                "UP",
                "auto_generated",
                "solar",
                "fixed_tilt_29_azim_180_solar_generation_20_years.csv",
            ),
            "r",
            encoding="UTF-8",
        ) as solar_twenty_year_file:
            solar_data = pd.read_csv(solar_twenty_year_file)
    except FileNotFoundError:
        with open(
            os.path.join(
                os.path.expanduser("~"),
                "fixed_tilt_29_azim_180_solar_generation_20_years.csv",
            ),
            "r",
            encoding="UTF-8",
        ) as solar_twenty_year_file:
            solar_data = pd.read_csv(solar_twenty_year_file)

    _plot_twenty_year_heatmap(
        solar_data["electricity"],
        axis=axes[0, 0],
        cmap_name="Reds",
        label="Electricity generation / kWh/kW$_\mathrm{p}$",
        vmin=0,
        vmax=1,
        title="Solar generation",
    )

    # Plot the grid-availability heatmap.
    with open(
        os.path.join(
            "grid_status_csvs",
            "UP_GRID_REGENERATED_13.csv",
        ),
        "r",
        encoding="UTF-8",
    ) as all_grid_availability_file:
        grid_data = pd.read_csv(all_grid_availability_file, index_col=0)

    _plot_twenty_year_heatmap(
        grid_data["13"],
        axis=axes[0, 1],
        cmap_name="Blues",
        # cmap_name="Blues"
        label="Grid-availability probability",
        title="National grid availability",
        vmin=0,
        vmax=1,
    )

    # Plot the load heatmap.
    try:
        with open(
            os.path.join(
                os.path.expanduser("~"),
                "clover_locations",
                "UP",
                "auto_generated",
                "load",
                "electric",
                "device_load",
                "total_load.csv",
            ),
            "r",
            encoding="UTF-8",
        ) as load_twenty_year_file:
            load_data = pd.read_csv(load_twenty_year_file)
    except FileNotFoundError:
        with open(
            os.path.join(
                "total_load.csv",
            ),
            "r",
            encoding="UTF-8",
        ) as load_twenty_year_file:
            load_data = pd.read_csv(load_twenty_year_file)

    _plot_twenty_year_heatmap(
        load_data["domestic"] + load_data["commercial"] + load_data["public"],
        axis=axes[1, 0],
        cmap_name="Greens",
        label="Average community load / kWh",
        title="Average total community energy demand",
    )

    # Plot fractional enterprise load
    _plot_twenty_year_heatmap(
        round(
            100
            * load_data["commercial"]
            / (load_data["domestic"] + load_data["commercial"] + load_data["public"])
        ),
        axis=axes[1, 1],
        cmap_name="Oranges",
        label="Fractional enterprise load",
        title="Fractional enterprise energy demand",
        vmin=0,
    )

    # Add annotation letters.
    (axis := axes[0, 0]).text(
        (_label_x_position := -0.12),
        (_label_y_position := 1.2),
        "a.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )
    (axis := axes[0, 1]).text(
        _label_x_position,
        _label_y_position,
        "b.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )
    (axis := axes[1, 0]).text(
        _label_x_position,
        _label_y_position,
        "c.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )
    (axis := axes[1, 1]).text(
        _label_x_position,
        _label_y_position,
        "d.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )

    plt.savefig(
        f"solar_heatmap_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )
    # plt.show()

    sns.set_context("notebook")


def _sanitise_data(
    data: pd.DataFrame,
    grid_only_simulations: pd.DataFrame | None = None,
    subsidised: bool = False,
) -> pd.DataFrame:
    """
    Sanitise the data from simulations that ran.

    :param: data
        The data to sanitise.

    :param: grid_only_simulations
        The data containing grid-only simulations.

    :param: subsidised
        Whether the grid cost should use the subsidised (True) or unsubsidised (False)
        values.

    """

    data = data[data["service_hours"] != 0]

    # Add the relevant grid-only data to the plots.
    if grid_only_simulations is not None:
        grid_only_simulations = grid_only_simulations[
            grid_only_simulations["grid_hours"] != 0
        ]
        grid_only_simulations["service_hours"] = grid_only_simulations["grid_hours"]

        # Add the simulation runs to the heatmaps
        for prioritise_self_generation in {True, False}:
            for grid_hours in range(1, 25):
                # Locate the matching result that ran as a simulation
                relevant_grid_simulation = grid_only_simulations.loc[
                    (grid_only_simulations["grid_hours"] == grid_hours)
                    & (
                        grid_only_simulations["prioritise_self_generation"]
                        == prioritise_self_generation
                    )
                ]

                # Ensure only usable columns are contained
                relevant_grid_simulation = relevant_grid_simulation[
                    data.columns.intersection(relevant_grid_simulation.columns)
                ]

                # Append service-hour data
                relevant_grid_simulations_with_service_hours: list[pd.DataFrame] = []
                for service_hours in range(1, grid_hours + 1):
                    _temp_row = relevant_grid_simulation.copy()
                    _temp_row["service_hours"] = service_hours
                    # _temp_row["unmet_energy_fraction"] = 0
                    _temp_row["subsidised"] = (subsidised,)
                    _temp_row["grid_fraction"] = 1 - _temp_row["renewables_fraction"]
                    relevant_grid_simulations_with_service_hours.append(_temp_row)

                # Add this result to the current data.
                data = pd.concat(
                    [
                        data,
                    ]
                    + relevant_grid_simulations_with_service_hours,
                    ignore_index=True,
                )

    return data.replace([np.inf, -np.inf, np.nan], np.nan)


def _setup_seaborn() -> None:
    """Setup seaborn for plotting nicely."""

    sns.set_context("notebook")
    sns.set_style("whitegrid")

    colorblind_palette = sns.color_palette(
        [
            "#E04606",  # Orange
            "#F09F52",  # Pale orange
            "#52C0AD",  # Pale green
            "#006264",  # Green
            "#D8247C",  # Pink
            "#EDEDED",  # Pale pink
            "#E7DFBE",  # Pale yellow
            "#FBBB2C",  # Yellow
        ]
    )

    sns.set_palette(colorblind_palette)


def main(args: list[Any]) -> None:
    """
    Main method.

    :param: args
        Unparsed command-line arguments.

    """

    parsed_args = _parse_args(args)

    #################################
    # Data parsing and sanitisation #
    #################################

    # Open the subsidised output file
    if (subsidised_output_filename := parsed_args.subsidised_output_file) is None:
        print(message := "Subsidised output file must be specified.")
        raise Exception(message)

    with open(
        subsidised_output_filename, "r", encoding="UTF-8"
    ) as subsidised_output_file:
        subsidised_data = pd.read_csv(subsidised_output_file, index_col=0)

    # Open the unsubsidised output file
    if (unsubsidised_output_filename := parsed_args.unsubsidised_output_file) is None:
        print(message := "Unsubsidised output file must be specified.")
        raise Exception(message)

    with open(
        unsubsidised_output_filename, "r", encoding="UTF-8"
    ) as unsubsidised_output_file:
        unsubsidised_data = pd.read_csv(unsubsidised_output_file, index_col=0)

    # Open the stranded-assets output file
    if (
        stranded_assets_output_filename := parsed_args.stranded_assets_output_file
    ) is None:
        print(message := "Standed-assets output file must be specified.")
        raise Exception(message)

    with open(
        stranded_assets_output_filename, "r", encoding="UTF-8"
    ) as stranded_assets_output_file:
        stranded_assets_data = pd.read_csv(stranded_assets_output_file, index_col=0)

    # Open the grid-only output file
    if (grid_output_filename := parsed_args.grid_output_file) is None:
        print(message := "Grid-only output file must be specified.")
        raise Exception(message)

    with open(grid_output_filename, "r", encoding="UTF-8") as output_file:
        grid_only_simulations = pd.read_csv(output_file, index_col=0)

    # Setup seaborn
    _setup_seaborn()

    # Plot box violin plots

    # Add service hours column
    # subsidised_data["service_hours"] = round(
    #     (24 * (1 - subsidised_data["blackouts"])), 0
    # )
    # subsidised_data["service_hours"] = subsidised_data["service_hours"].astype(int)
    # subsidised_data["service_hours_provided"] = subsidised_data["service_hours"]
    # unsubsidised_data["service_hours"] = round(
    #     24 * (1 - unsubsidised_data["blackouts"]), 0
    # )
    # unsubsidised_data["service_hours"] = unsubsidised_data["service_hours"].astype(int)
    # unsubsidised_data["service_hours_provided"] = unsubsidised_data["service_hours"]

    stranded_assets_data["grid_fraction"] = (
        1 - stranded_assets_data["renewables_fraction"]
    )
    subsidised_data["grid_fraction"] = 1 - subsidised_data["renewables_fraction"]
    unsubsidised_data["grid_fraction"] = 1 - unsubsidised_data["renewables_fraction"]

    # Sort the data based on the grid hours, service hours, and self-prioritisation
    stranded_assets_data = stranded_assets_data.sort_values(
        by=["prioritise_self_generation", "grid_hours", "service_hours"],
        ascending=[True, True, True],
        na_position="first",
    )
    subsidised_data = subsidised_data.sort_values(
        by=["prioritise_self_generation", "grid_hours", "service_hours"],
        ascending=[True, True, True],
        na_position="first",
    )
    unsubsidised_data = unsubsidised_data.sort_values(
        by=["prioritise_self_generation", "grid_hours", "service_hours"],
        ascending=[True, True, True],
        na_position="first",
    )

    # Remove duplicate stranded-assets results.
    stranded_assets_data = stranded_assets_data.drop_duplicates()

    # Reset indicies on the dataframes
    stranded_assets_data = stranded_assets_data.reset_index(drop=True)
    subsidised_data = subsidised_data.reset_index(drop=True)
    unsubsidised_data = unsubsidised_data.reset_index(drop=True)

    # Add columns describing the subsidisation state.
    stranded_assets_data[(_subsidised := "subsidised")] = True
    subsidised_data[_subsidised] = True
    unsubsidised_data[_subsidised] = False

    #####################################################################
    # Figure 3. LCUE and emissions premium for additional service hours #
    #####################################################################

    bar_plot_data = subsidised_data.copy()
    bar_plot_data = _sanitise_data(
        bar_plot_data, grid_only_simulations.copy(), subsidised=True
    )
    bar_plot_data = bar_plot_data[
        (bar_plot_data["prioritise_self_generation"] == False)
        & (bar_plot_data["service_hours"] >= 1)
    ]
    bar_plot_data = bar_plot_data.reset_index(drop=True)

    bar_plot_data["lcue_premium"] = None
    bar_plot_data["emissions_intensity_premium"] = None
    for index, _ in enumerate(bar_plot_data.iterrows()):
        if (this_grid_hours := bar_plot_data.loc[index, "grid_hours"]) >= (
            this_service_hours := bar_plot_data.loc[index, "service_hours"]
        ):
            bar_plot_data.loc[index, "lcue_premium"] = 0.0
            bar_plot_data.loc[index, "emissions_intensity_premium"] = 0.0
            continue
        try:
            bar_plot_data.loc[index, "lcue_premium"] = float(
                100
                * (
                    bar_plot_data.loc[index, "lcue"]
                    - bar_plot_data[
                        (bar_plot_data["grid_hours"] == int(this_grid_hours))
                        & (bar_plot_data["service_hours"] == int(this_grid_hours))
                    ]["lcue"]
                ).to_list()[0]
                / (this_service_hours - this_grid_hours)
            )
            bar_plot_data.loc[index, "emissions_intensity_premium"] = float(
                (
                    bar_plot_data.loc[index, "emissions_intensity"]
                    - bar_plot_data[
                        (bar_plot_data["grid_hours"] == int(this_grid_hours))
                        & (bar_plot_data["service_hours"] == int(this_grid_hours))
                    ]["emissions_intensity"]
                ).to_list()[0]
                / (this_service_hours - this_grid_hours)
            )
        except IndexError:
            continue

    bar_plot_data["lcue_premium"] = bar_plot_data["lcue_premium"].astype(float)
    bar_plot_data["emissions_intensity_premium"] = bar_plot_data[
        "emissions_intensity_premium"
    ].astype(float)

    def _colours_from_values(values, palette_name):
        """
        Function sourced from: https://stackoverflow.com/questions/36271302/changing-color-scale-in-seaborn-bar-plot

        """
        # normalize the values to range [0, 1]
        normalized = (values - min(values)) / (max(values) - min(values))
        # convert to indices
        indices = np.round(normalized * (len(values) - 1)).astype(np.int32)
        # use the indices to get the colors
        palette = sns.color_palette(palette_name, len(values))
        return np.array(palette).take(indices, axis=0)

    fig, axes = plt.subplots(1, 2, figsize=(48 / 5, 32 / 5))
    sns.set_style("ticks")

    bar_plot_row = bar_plot_data[
        (bar_plot_data["service_hours"] == bar_plot_data["grid_hours"] + 2)
        & (bar_plot_data["grid_hours"] >= 8)
    ]

    sns.barplot(
        x=bar_plot_row["grid_hours"],
        y=bar_plot_row["lcue_premium"],
        # palette=(_palette:=_colours_from_values(bar_plot_row["lcue_premium"], "PiYG_r")),
        palette=(
            _palette := sns.color_palette(
                "PiYG_r",
                as_cmap=True,
            )
        ),
        hue=bar_plot_row["lcue_premium"],
        ax=(axis := axes[0]),
        hue_norm=(
            _normalize := mcolors.TwoSlopeNorm(
                vcenter=0,
                vmin=min(bar_plot_row["lcue_premium"]),
                vmax=max(bar_plot_row["lcue_premium"]),
            )
        ),
    )

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("LCUE premium for additional\nservice hours / cents/kWh/hour")

    axis.set_xlim(-0.5, 14.5)
    axis.set_ylim(-7, 5)

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.get_legend().remove()

    bar_plot_row = bar_plot_data[
        (bar_plot_data["service_hours"] == bar_plot_data["grid_hours"] + 4)
        & (bar_plot_data["grid_hours"] >= 8)
    ]

    sns.barplot(
        x=bar_plot_row["grid_hours"],
        y=bar_plot_row["lcue_premium"],
        palette=_palette,
        hue=bar_plot_row["lcue_premium"],
        ax=(axis := axes[1]),
        hue_norm=_normalize,
    )

    axis.set_xlim(-0.5, 12.5)
    axis.set_ylim(-7, 5)

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("")

    axis.get_legend().remove()

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    plt.savefig(
        f"lcue_premium_for_two_and_four_hours_small_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    line_plot_hours: dict[int, Any] = {}
    for hour in range(1, 11):
        line_plot_hours[hour] = bar_plot_data[
            (bar_plot_data["service_hours"] == bar_plot_data["grid_hours"] + hour)
            & (bar_plot_data["grid_hours"] >= 8)
        ]

    sns.set_palette(
        sns.color_palette(
            [
                "#8C3777",
                "#9E4671",
                "#B0546A",
                "#BC5950",
                "#C75D36",
                "#CC7E4F",
                "#D09E67",
                "#C7B27D",
                "#BEC593",
                "#CACEC6",
            ]
        )
    )

    fig, axes = plt.subplots(
        4,
        1,
        height_ratios=[2, 2, 0.05, 2],
        figsize=(48 / 5, 64 / 5),
        layout="constrained",
    )
    axes[2].set_visible(False)
    axes[3].set_visible(False)
    sns.set_style("ticks")

    plt.subplot(4, 1, 1)
    palettes = [
        sns.cubehelix_palette(rot=-0.2, n_colors=10, reverse=False),
        sns.cubehelix_palette(rot=-0.6, n_colors=10, reverse=False),
    ]
    palette = sns.color_palette(palettes[0].as_hex() + palettes[1].as_hex())
    markers = ["h", "H"] * 5

    sns.set_palette(palette)
    for hour, line_plot_data in line_plot_hours.items():
        sns.scatterplot(
            x=line_plot_data["grid_hours"],
            y=line_plot_data["lcue_premium"],
            marker=markers[hour - 1],
            alpha=(_alpha := 0.8),
            color=f"C{hour - 1}",
            s=400,
            linewidth=0,
            label=hour,
            ax=(axis := plt.gca()),
            palette=palette,
        )
        sns.lineplot(
            x=line_plot_data["grid_hours"],
            y=line_plot_data["lcue_premium"],
            color=f"C{hour - 1}",
            linewidth=2,
            alpha=_alpha,
            ax=axis,
            palette=palette,
        )

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("LCUE premium for additional service\nhours / cents/kWh/hour")

    axis.axhline(y=0, color="grey", linestyle="--")

    axis.set_xlim(7.5, 23.5)
    axis.set_ylim(-13, 8)

    axis.get_legend().remove()

    axis.margins(tight=True)

    axis.text(
        (_label_x_position := -0.1),
        (_label_y_position := 1.1),
        "a.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    norm = plt.Normalize(0.5, 10.5)
    scalar_mappable = plt.cm.ScalarMappable(
        cmap=mcolors.LinearSegmentedColormap.from_list(
            "Custom", palettes[0].as_hex(), 10
        ),
        norm=norm,
    )

    colorbar = axis.figure.colorbar(
        scalar_mappable,
        ax=axis,
        label="Additional service hours provided beyond\nthose supplied by the grid",
        pad=(_pad := 0.025),
    )

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    plt.subplot(4, 1, 2)
    for hour, line_plot_data in line_plot_hours.items():
        sns.scatterplot(
            x=line_plot_data["grid_hours"],
            y=line_plot_data["emissions_intensity_premium"],
            marker=markers[hour - 1],
            alpha=(_alpha := 0.8),
            color=f"C{10 + hour - 1}",
            s=400,
            linewidth=0,
            label=hour,
            ax=(axis := plt.gca()),
            palette=palette,
        )
        sns.lineplot(
            x=line_plot_data["grid_hours"],
            y=line_plot_data["emissions_intensity_premium"],
            color=f"C{10 + hour - 1}",
            linewidth=2,
            alpha=_alpha,
            ax=axis,
            palette=palette,
        )

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel(
        "Emissions intensity premium for additional\nservice hours / gCO2/kWh/hour"
    )

    axis.axhline(y=0, color="grey", linestyle="--")

    axis.set_xlim(7.5, 23.5)

    axis.get_legend().remove()

    axis.margins(tight=True)

    axis.text(
        _label_x_position,
        _label_y_position,
        "b.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    norm = plt.Normalize(0.5, 10.5)
    scalar_mappable = plt.cm.ScalarMappable(
        cmap=mcolors.LinearSegmentedColormap.from_list(
            "Custom", palettes[1].as_hex(), 10
        ),
        norm=norm,
    )

    colorbar = axis.figure.colorbar(
        scalar_mappable,
        ax=axis,
        label="Additional service hours provided beyond\nthose supplied by the grid",
        pad=_pad,
    )

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    plt.subplot(4, 2, 7)
    bar_plot_row = bar_plot_data[
        (bar_plot_data["service_hours"] == bar_plot_data["grid_hours"] + 2)
        & (bar_plot_data["grid_hours"] >= 8)
    ]

    sns.barplot(
        x=[int(entry) for entry in bar_plot_row["grid_hours"]],
        y=bar_plot_row["lcue_premium"],
        # palette=(_palette:=_colours_from_values(bar_plot_row["lcue_premium"], "PiYG_r")),
        palette=(
            _palette := sns.color_palette(
                "PiYG_r",
                as_cmap=True,
            )
        ),
        hue=bar_plot_row["lcue_premium"],
        ax=(axis := plt.gca()),
        hue_norm=(
            _normalize := mcolors.TwoSlopeNorm(
                vcenter=0,
                vmin=min(bar_plot_row["lcue_premium"]),
                vmax=max(bar_plot_row["lcue_premium"]),
            )
        ),
    )

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("LCUE premium for additional\nservice hours / cents/kWh/hour")

    axis.set_xlim(-0.5, 14.5)
    axis.set_ylim(-7, 5)

    axis.text(
        (_label_x_position := -0.075),
        (_label_y_position := 1.05),
        "c.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.get_legend().remove()

    plt.subplot(4, 2, 8)
    bar_plot_row = bar_plot_data[
        (bar_plot_data["service_hours"] == bar_plot_data["grid_hours"] + 4)
        & (bar_plot_data["grid_hours"] >= 8)
    ]

    sns.barplot(
        x=[int(entry) for entry in bar_plot_row["grid_hours"]],
        y=bar_plot_row["lcue_premium"],
        # palette=(_palette:=_colours_from_values(bar_plot_row["lcue_premium"], "PiYG_r")),
        palette=(_palette),
        hue=bar_plot_row["lcue_premium"],
        ax=(axis := plt.gca()),
        hue_norm=(_normalize),
    )

    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("LCUE premium for additional\nservice hours / cents/kWh/hour")

    axis.set_xlim(-0.5, 12.5)
    axis.set_ylim(-7, 5)

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.text(
        _label_x_position,
        _label_y_position,
        "d.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    axis.get_legend().remove()

    plt.savefig(
        f"lcue_premium_for_all_hours_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    #################################################################################
    # Alternative Figure 3. Sparkplot for each individual hour, rather than scatter #
    #################################################################################

    fig, axes = plt.subplots(len(line_plot_hours), 1, figsize=(32 / 5, 48 / 5))

    min_premium = min(
        bar_plot_data[bar_plot_data["grid_hours"] >= 8]["lcue_premium"].dropna()
    )
    max_premium = max(
        bar_plot_data[bar_plot_data["grid_hours"] >= 8]["lcue_premium"].dropna()
    )

    for hour, line_plot_data in line_plot_hours.items():
        sns.lineplot(
            x=line_plot_data["grid_hours"],
            y=line_plot_data["lcue_premium"],
            label=f"{hour} hours",
            ax=(axis := axes[hour - 1]),
            color=f"C{hour-1}",
        )
        axis.axhline(c="grey", alpha=0.5)
        axis.axes.get_xaxis().get_label().set_visible(False)
        axis.axes.get_yaxis().get_label().set_visible(False)
        axis.set_xlim(8, 24)
        axis.set_ylim(min_premium, max_premium)
        axis.fill_between(
            line_plot_data["grid_hours"],
            [0] * len(line_plot_data["lcue_premium"]),
            line_plot_data["lcue_premium"],
            hatch="//",
            color=f"C{hour-1}",
            alpha=0.3,
            where=line_plot_data["lcue_premium"] > 0,
        )
        axis.fill_between(
            line_plot_data["grid_hours"],
            line_plot_data["lcue_premium"],
            [0] * len(line_plot_data["lcue_premium"]),
            hatch="\\\\",
            color=f"C{hour-1}",
            alpha=0.3,
            where=line_plot_data["lcue_premium"] < 0,
        )
        plt.setp(axis.get_xticklabels(), visible=False)
        plt.setp(axis.get_yticklabels(), visible=False)
        plt.setp(axis.get_xticklines(), visible=False)
        plt.setp(axis.get_yticklines(), visible=False)
        plt.setp(axis.spines.values(), visible=False)
        axis.get_legend().remove()

    # bottom sparkline
    plt.setp(axes[-1].get_xticklabels(), visible=True)
    plt.setp(axes[-1].get_xticklines(), visible=True)
    axes[-1].xaxis.tick_bottom()  # but onlyt the lower x ticks not x ticks at the top

    plt.tight_layout()

    plt.savefig(
        f"lcue_premium_spark_plot_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    ##################################################################
    # Figure 1. System size, LCUE and emissions-intensity heat maps. #
    ##################################################################

    plotting_data = subsidised_data.copy()
    plotting_data = _sanitise_data(
        plotting_data, grid_only_simulations.copy(), subsidised=True
    )
    plotting_data = plotting_data[
        (plotting_data["grid_hours"] % 2 == 0)
        & (8 <= plotting_data["grid_hours"])
        & (plotting_data["service_hours"] % 2 == 0)
        & (plotting_data["service_hours"] >= 8)
        & (plotting_data["prioritise_self_generation"] == False)
    ]
    plotting_data = plotting_data.reset_index(drop=True)

    plotting_data["lcue_premium"] = None
    for index, _ in enumerate(plotting_data.iterrows()):
        if (this_grid_hours := plotting_data.loc[index, "grid_hours"]) > (
            this_service_hours := plotting_data.loc[index, "service_hours"]
        ):
            plotting_data.loc[index, "lcue_premium"] = 0.0
            continue
        if this_grid_hours == this_service_hours:
            plotting_data.loc[index, "lcue_premium"] = 0.0
            continue
        plotting_data.loc[index, "lcue_premium"] = float(
            100
            * (
                plotting_data.loc[index, "lcue"]
                - plotting_data[
                    (plotting_data["grid_hours"] == this_grid_hours)
                    & (plotting_data["service_hours"] == this_grid_hours)
                ]["lcue"]
            ).to_list()[0]
            / (this_service_hours - this_grid_hours)
        )

    plotting_data["lcue_premium"] = plotting_data["lcue_premium"].astype(float)

    _, axes = plt.subplots(2, 2, figsize=(48 / 5, 48 / 5))

    sns.heatmap(
        plotting_data.pivot(
            index="service_hours", columns="grid_hours", values="initial_pv_size"
        ),
        annot=True,
        ax=(axis := axes[0, 0]),
        cmap="Blues",
        fmt=".2g",
        linewidths=1,
        square=True,
        cbar_kws={"label": "Initial PV capacity / kWp"},
    )
    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("Minimum service hours")
    axis.set_title("PV capacity")
    axis.margins(x=0, y=0)
    axis.invert_yaxis()
    axis.text(
        (_label_x_position := -0.08),
        (_label_y_position := 1.1),
        "a.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    sns.heatmap(
        plotting_data.pivot(
            index="service_hours", columns="grid_hours", values="initial_storage_size"
        ),
        annot=True,
        ax=(axis := axes[0, 1]),
        cmap="Greens",
        fmt=".2g",
        linewidths=1,
        square=True,
        cbar_kws={"label": "Initial PV capacity / kWp"},
    )
    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("Minimum service hours")
    axis.set_title("Storage capacity")
    axis.margins(x=0, y=0)
    axis.invert_yaxis()
    axis.text(
        _label_x_position,
        _label_y_position,
        "b.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )

    sns.heatmap(
        plotting_data.pivot(index="service_hours", columns="grid_hours", values="lcue"),
        annot=True,
        ax=(axis := axes[1, 0]),
        cmap="PuBu",
        fmt=".2g",
        linewidths=1,
        square=True,
        cbar_kws={"label": "LCUE / $/kWh"},
    )
    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("Minimum service hours")
    axis.set_title("LCUE")
    axis.margins(x=0, y=0)
    axis.invert_yaxis()
    axis.text(
        _label_x_position,
        _label_y_position,
        "c.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )

    sns.heatmap(
        plotting_data.pivot(
            index="service_hours", columns="grid_hours", values="emissions_intensity"
        ),
        annot=True,
        ax=(axis := axes[1, 1]),
        cmap="BuPu",
        fmt=".2g",
        linewidths=1,
        square=True,
        cbar_kws={"label": "Emissions intensity / gCO2eq/kWh"},
    )
    axis.set_xlabel("Hours of national grid availability")
    axis.set_ylabel("Minimum service hours")
    axis.set_title("Emissions intensity")
    axis.margins(x=0, y=0)
    axis.invert_yaxis()
    axis.text(
        _label_x_position,
        _label_y_position,
        "d.",
        transform=axis.transAxes,
        fontsize=_fontsize,
        fontweight="bold",
        va="top",
        ha="right",
    )

    axis.add_patch(
        plt.Polygon(
            (
                exceeded_ndc_columns := [
                    [8.95, 8.95],
                    [6, 8.95],
                    [6, 8],
                    [5, 8],
                    [5, 5],
                    [4, 5],
                    [4, 4],
                    [3, 4],
                    [3, 3],
                    [2, 3],
                    [2, 2],
                    [1, 2],
                    [1, 1],
                    [0.05, 1],
                    [0.05, 0.05],
                    [8.95, 0.05],
                ]
            ),
            # color="grey",
            color="#ffa600",
            fill=False,
            lw=2,
            label="Assets stranded by the arrival of the grid",
            # hatch="\\\\//",
            # alpha=0.4,
        )
    )

    plt.margins(x=0, y=0)

    plt.savefig(
        f"stranded_assets_heatmaps_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    #########################################################################
    # Figure 6. Electricity supply-and-demand heatmaps used for data input. #
    #########################################################################

    # Plot the solar generation
    _plot_solar_heatmap()

    ######################################
    # Figure 7. Impacts of grid subsidy. #
    ######################################

    # NOTE: Additional figures are included here beyond that shown in the S.M.
    #

    # Plot grid subsidisation figures
    fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))
    sns.violinplot(
        (
            data_change := (
                unsubsidised_data[
                    unsubsidised_data["prioritise_self_generation"] == False
                ].loc[:, ["initial_pv_size", "initial_storage_size"]]
                - subsidised_data[
                    subsidised_data["prioritise_self_generation"] == False
                ].loc[:, ["initial_pv_size", "initial_storage_size"]]
            )
        ),
        cut=0,
        fill=False,
        # linecolor="C0",
    )
    sns.stripplot(
        data_change,
        alpha=0.1,
        # color="C0",
        jitter=False,
        marker="D",
        size=20,
    )
    plt.xlabel("Component")
    ax.set_xticklabels(["PV panels", "Batteries"])
    plt.ylabel("Change in capacity / kW$_\mathrm{p}$ (PV), kWh (batteries)")
    plt.savefig(
        f"change_in_capacity_with_grid_subsidy_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))
    sns.violinplot(
        (
            data_change := (
                unsubsidised_data[
                    unsubsidised_data["prioritise_self_generation"] == True
                ].loc[:, ["initial_pv_size", "initial_storage_size"]]
                - subsidised_data[
                    subsidised_data["prioritise_self_generation"] == True
                ].loc[:, ["initial_pv_size", "initial_storage_size"]]
            )
        ),
        cut=0,
        fill=False,
        # linecolor="C0",
    )
    sns.stripplot(
        data_change,
        alpha=0.1,
        # color="C0",
        jitter=False,
        marker="D",
        size=20,
    )
    plt.xlabel("Component")
    ax.set_xticklabels(["PV panels", "Batteries"])
    plt.ylabel("Change in capacity / kW$_\mathrm{p}$ (PV), kWh (batteries)")
    plt.savefig(
        f"change_in_capacity_with_grid_subsidy_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    fig, ax = plt.subplots()
    sns.violinplot(
        (
            data_change := (
                unsubsidised_data[
                    unsubsidised_data["prioritise_self_generation"] == False
                ].loc[:, ["lcue"]]
                - subsidised_data[
                    subsidised_data["prioritise_self_generation"] == False
                ].loc[:, ["lcue"]]
            )
        ),
        cut=0,
        fill=False,
        # linecolor="C0",
    )
    sns.stripplot(
        data_change,
        alpha=0.1,
        # color="C0",
        jitter=False,
        marker="D",
        size=20,
    )
    plt.xlabel("System runs")
    ax.set_xticklabels([""])
    plt.ylabel("Change in LCUE / USD/kWh")
    plt.savefig(
        f"change_in_lcue_with_grid_subsidy_single_violin_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))
    sns.violinplot(
        data=(
            combined_frame := pd.concat([subsidised_data, unsubsidised_data], axis=0)
        ),
        x="grid_hours",
        y="lcue",
        bw_adjust=2,
        cut=0,
        fill=False,
        hue=_subsidised,
        inner="quart",
        # linecolor="C0",
        palette={True: "C3", False: "C0"},
        split=True,
    )
    sns.stripplot(
        data=combined_frame,
        x="grid_hours",
        y="lcue",
        alpha=0.2,
        # color="C0",
        dodge=True,
        hue=_subsidised,
        jitter=False,
        marker="D",
        palette={True: "C3", False: "C0"},
        size=5,
    )
    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("LCUE / USD/kWh")
    plt.savefig(
        f"lcue_with_grid_subsidy_violins_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    unsubsidised_data["percentage_lcue_change"] = (
        100
        * (unsubsidised_data["lcue"] - subsidised_data["lcue"])
        / subsidised_data["lcue"]
    )

    #################################################
    # Figure 7. Impacts of grid subsidy — diamonds. #
    #################################################

    stranded_assets_data["grid_hours"] = stranded_assets_data["grid_hours"].astype(int)
    stranded_assets_data["service_hours"] = stranded_assets_data[
        "service_hours"
    ].astype(int)
    unsubsidised_data["grid_hours"] = unsubsidised_data["grid_hours"].astype(int)
    unsubsidised_data["service_hours"] = unsubsidised_data["service_hours"].astype(int)

    sns.set_style("ticks", {"axes.grid": False})
    fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))
    unsubsidised_data["Grid fraction"] = unsubsidised_data["grid_fraction"]
    sns.stripplot(
        data=unsubsidised_data,
        x="grid_hours",
        y="percentage_lcue_change",
        alpha=0.2,
        # color="C0",
        # dodge=True,
        hue="Grid fraction",
        jitter=False,
        marker="D",
        # palette={True: "C3", False: "C0"},
        size=15,
    )
    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Percentage change in LCUE compared to a subsidised grid")
    plt.savefig(
        f"change_in_lcue_with_grid_subsidy_diamonds_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    #####################################################################
    # Data sanitisation for the impact of self- vs grid-prioritisation. #
    #####################################################################

    # Sanitise with the zero-hours information and grid-only information
    unsubsidised_data = _sanitise_data(unsubsidised_data)
    subsidised_data = _sanitise_data(
        subsidised_data, grid_only_simulations, subsidised=True
    )

    stranded_assets_data["service_hours_provided"] = 24 * (
        1 - stranded_assets_data["blackouts"]
    )
    subsidised_data["service_hours_provided"] = 24 * (1 - subsidised_data["blackouts"])

    # Compute the increase in LCUE compared to optimised systems
    for index, row in tqdm(enumerate(stranded_assets_data.iterrows())):
        # Compute and save the change in LCUE
        optimised_lcue = subsidised_data[
            (
                subsidised_data["grid_hours"]
                == stranded_assets_data.loc[index, "grid_hours"]
            )
            & (
                subsidised_data["service_hours"]
                == stranded_assets_data.loc[index, "service_hours"]
            )
            & (
                subsidised_data["prioritise_self_generation"]
                == stranded_assets_data.loc[index, "prioritise_self_generation"]
            )
        ]["lcue"]
        try:
            stranded_assets_data.loc[index, "change_in_lcue"] = float(
                stranded_assets_data.loc[index, "lcue"] - optimised_lcue
            )
        except Exception:
            pass

        # Compute and save the change in LCUE as a result of changing prioritisation
        this_lcue = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    == subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["lcue"]
        )
        other_lcue = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    != subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["lcue"]
        )

        # If prioritising the grid, then we use this value as the grid value
        if subsidised_data.loc[index, "prioritise_self_generation"] == False:
            grid_to_self_lcue_percentage_change = (
                100 * (other_lcue - this_lcue) / this_lcue
            )
        else:
            grid_to_self_lcue_percentage_change = (
                100 * (this_lcue - other_lcue) / other_lcue
            )
        try:
            subsidised_data.loc[index, "grid_to_self_lcue_percentage_change"] = float(
                grid_to_self_lcue_percentage_change
            )
        except Exception:
            pass
        # Compute and save the change in LCUE as a result of changing prioritisation
        this_total_system_cost = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    == subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["total_system_cost"]
        )
        other_total_system_cost = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    != subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["total_system_cost"]
        )

        # If prioritising the grid, then we use this value as the grid value
        if subsidised_data.loc[index, "prioritise_self_generation"] == False:
            grid_to_self_total_system_cost_percentage_change = (
                100
                * (other_total_system_cost - this_total_system_cost)
                / this_total_system_cost
            )
        else:
            grid_to_self_total_system_cost_percentage_change = (
                100
                * (this_total_system_cost - other_total_system_cost)
                / other_total_system_cost
            )
        try:
            subsidised_data.loc[
                index, "grid_to_self_total_system_cost_percentage_change"
            ] = float(grid_to_self_total_system_cost_percentage_change)
        except Exception:
            pass
        # Compute and save the change in LCUE as a result of changing prioritisation
        this_service_hours_provided = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    == subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["service_hours_provided"]
        )
        other_service_hours_provided = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    != subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["service_hours_provided"]
        )

        # If prioritising the grid, then we use this value as the grid value
        if subsidised_data.loc[index, "prioritise_self_generation"] == False:
            grid_to_self_service_hours_provided_percentage_change = (
                other_service_hours_provided - this_service_hours_provided
            ) / this_service_hours_provided
        else:
            grid_to_self_service_hours_provided_percentage_change = (
                this_service_hours_provided - other_service_hours_provided
            ) / other_service_hours_provided
        try:
            subsidised_data.loc[
                index, "grid_to_self_service_hours_provided_percentage_change"
            ] = float(grid_to_self_service_hours_provided_percentage_change)
        except Exception:
            pass
        this_unmet_energy_fraction = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    == subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["unmet_energy_fraction"]
        )
        other_unmet_energy_fraction = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    != subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["unmet_energy_fraction"]
        )

        # If prioritising the grid, then we use this value as the grid value
        if subsidised_data.loc[index, "prioritise_self_generation"] == False:
            grid_to_self_unmet_energy_fraction_change = 100 * (
                other_unmet_energy_fraction - this_unmet_energy_fraction
            )
        else:
            grid_to_self_unmet_energy_fraction_change = 100 * (
                this_unmet_energy_fraction - other_unmet_energy_fraction
            )
        try:
            subsidised_data.loc[index, "grid_to_self_unmet_energy_fraction_change"] = (
                float(grid_to_self_unmet_energy_fraction_change)
            )
        except Exception:
            pass

        this_emissions_intensity = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    == subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["emissions_intensity"]
        )
        other_emissions_intensity = float(
            subsidised_data[
                (
                    subsidised_data["grid_hours"]
                    == subsidised_data.loc[index, "grid_hours"]
                )
                & (
                    subsidised_data["service_hours"]
                    == subsidised_data.loc[index, "service_hours"]
                )
                & (
                    subsidised_data["prioritise_self_generation"]
                    != subsidised_data.loc[index, "prioritise_self_generation"]
                )
            ]["emissions_intensity"]
        )

        # If prioritising the grid, then we use this value as the grid value
        if subsidised_data.loc[index, "prioritise_self_generation"] == False:
            grid_to_self_emissions_intensity_percentage_change = (
                100
                * (other_emissions_intensity - this_emissions_intensity)
                / this_emissions_intensity
            )
        else:
            grid_to_self_emissions_intensity_percentage_change = (
                100
                * (this_emissions_intensity - other_emissions_intensity)
                / other_emissions_intensity
            )
        try:
            subsidised_data.loc[
                index, "grid_to_self_emissions_intensity_percentage_change"
            ] = float(grid_to_self_emissions_intensity_percentage_change)
        except Exception:
            pass

        # Compute and save the change in service hours provided
        if (
            stranded_assets_data.loc[index, "grid_hours"]
            >= stranded_assets_data.loc[index, "service_hours"]
        ):
            baseline_service_hours = stranded_assets_data.loc[index, "grid_hours"]
        else:
            baseline_service_hours = subsidised_data[
                subsidised_data["grid_hours"]
                == stranded_assets_data.loc[index, "grid_hours"]
            ][
                subsidised_data["service_hours"]
                == stranded_assets_data.loc[index, "service_hours"]
            ][
                subsidised_data["prioritise_self_generation"]
                == stranded_assets_data.loc[index, "prioritise_self_generation"]
            ][
                "service_hours"
            ]
        try:
            stranded_assets_data.loc[index, "additional_service_hours"] = float(
                24 * (1 - stranded_assets_data.loc[index, "blackouts"])
                - baseline_service_hours
            )
            stranded_assets_data.loc[index, "service_hours_provided"] = float(
                stranded_assets_data.loc[index, "additional_service_hours"]
                + stranded_assets_data.loc[index, "service_hours"]
            )
        except Exception:
            pass

        # Compute and save the change in emissions intensity
        optimised_emissions_intensity = subsidised_data[
            subsidised_data["grid_hours"]
            == stranded_assets_data.loc[index, "grid_hours"]
        ][
            subsidised_data["service_hours"]
            == stranded_assets_data.loc[index, "service_hours"]
        ][
            subsidised_data["prioritise_self_generation"]
            == stranded_assets_data.loc[index, "prioritise_self_generation"]
        ][
            "emissions_intensity"
        ]
        try:
            stranded_assets_data.loc[index, "change_in_emissions_intensity"] = float(
                stranded_assets_data.loc[index, "emissions_intensity"]
                - optimised_emissions_intensity
            )
        except Exception:
            pass
        # Compute and save the change in unmet energy fraction
        optimised_unmet_energy_fraction = subsidised_data[
            subsidised_data["grid_hours"]
            == stranded_assets_data.loc[index, "grid_hours"]
        ][
            subsidised_data["service_hours"]
            == stranded_assets_data.loc[index, "service_hours"]
        ][
            subsidised_data["prioritise_self_generation"]
            == stranded_assets_data.loc[index, "prioritise_self_generation"]
        ][
            "unmet_energy_fraction"
        ]
        try:
            stranded_assets_data.loc[index, "change_in_unmet_energy_fraction"] = float(
                stranded_assets_data.loc[index, "unmet_energy_fraction"]
                - optimised_unmet_energy_fraction
            )
        except Exception:
            pass

    ###############################################
    # LCUE of systems undergoing asset stranding. #
    ###############################################

    _heatmap_plot(
        stranded_assets_data,
        "lcue",
        cmap_name="Reds",
        default_value=0,
        prioritise_self_generation=False,
        ylabel="LCUE / USD/kWh",
    )
    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    # plt.show()

    ##########################################################################
    # Figure 8. Service hours provided when systems undergo asset stranding. #
    ##########################################################################

    # Figure 8a. — Heatmap
    #

    # Plot stranded-assets figure
    ax = _heatmap_plot(
        (_plotting_data := stranded_assets_data),
        "service_hours_provided",
        center=12,
        cmap_name="PuBu",
        default_value=0,
        figsize=(64 / 5, 24 / 5),
        prioritise_self_generation=False,
        ylabel="Service hours provided by assets installed\nfor 13 hours of grid availability",
    )
    _heatmap_plot(
        _plotting_data,
        "service_hours_provided",
        annot=True,
        ax=ax,
        center=12,
        cmap_name="PuBu",
        default_value=0,
        figsize=(48 / 5, 32 / 5),
        prioritise_self_generation=False,
        ylabel="Service hours provided by assets installed\nfor 13 hours of grid availability",
    )

    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    plt.savefig(
        f"service_hours_asset_stranding_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    # Figure 8b. — Line plot.

    line_plot_service_hours: dict[int, Any] = {}

    stranded_assets_data["optimum_system_service_hours_provided"] = None
    for service_hours in range(14, 25):
        for grid_hours in range(0, 25):
            # Skip entries not included within the subsidised data set.
            if grid_hours >= service_hours:
                continue
            for self_generation in [True, False]:
                stranded_assets_data.loc[
                    (
                        stranded_assets_data["prioritise_self_generation"]
                        == self_generation
                    )
                    & (stranded_assets_data["grid_hours"] == grid_hours)
                    & (stranded_assets_data["service_hours"] == service_hours),
                    "optimum_system_service_hours_provided",
                ] = subsidised_data.loc[
                    (subsidised_data["grid_hours"] == grid_hours)
                    & (subsidised_data["prioritise_self_generation"] == self_generation)
                    & (subsidised_data["service_hours"] == service_hours),
                    "service_hours_provided",
                ].iloc[
                    0
                ]

    for hour in range(14, 25):
        line_plot_service_hours[hour] = stranded_assets_data[
            (stranded_assets_data["service_hours"] == hour)
            & (stranded_assets_data["prioritise_self_generation"] == False)
        ]

    sns.set_palette(
        (
            palette := sns.color_palette(
                [
                    "#003D30",
                    "#005745",
                    "#00735C",
                    "#009175",
                    "#00CBA7",
                    "#86FFDE",
                    "#FFCCFE",
                    "#FF92FD",
                    "#A700FC",
                    "#65019F",
                    "#450270",
                ]
            )
        )
    )

    fig = plt.figure(figsize=(64 / 5, 24 / 5))
    sns.set_style("ticks")

    for hour, line_plot_data in line_plot_service_hours.items():
        sns.lineplot(
            x=line_plot_data["grid_hours"],
            y=(
                line_plot_data["service_hours_provided"]
                - line_plot_data["optimum_system_service_hours_provided"]
            ),
            label=f"{hour} hours",
        )

    plt.xlabel("Hours of national grid availability")
    plt.ylabel("Difference in service hours provided compared to optimum system")

    plt.xticks(list(range(0, 25, 4)))
    plt.yticks(list(range(-14, 6, 2)))

    plt.axhline(y=0, color="grey", linestyle="--")
    plt.axvline(x=13, color="grey", linestyle="--")

    plt.xlim(0, 24)

    plt.savefig(
        f"additional_service_hours_for_all_hours_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    # Line plot attempt 2
    sns.set_palette(
        (
            palette := sns.color_palette(
                [
                    "#003D30",
                    "#005745",
                    "#00735C",
                    "#009175",
                    "#00CBA7",
                    "#86FFDE",
                    "#FFCCFE",
                    "#FF92FD",
                    "#A700FC",
                    "#65019F",
                    "#450270",
                ]
            )
        )
    )

    fig = plt.figure(figsize=(58 / 5, 24 / 5))
    sns.set_style("ticks")

    sns.lineplot(
        stranded_assets_data[
            stranded_assets_data["prioritise_self_generation"] == False
        ],
        x="grid_hours",
        y="additional_service_hours",
        hue="service_hours",
        palette=palette,
    )

    plt.xlabel("Hours of national grid availability")
    plt.ylabel("Difference in service hours provided\ncompared to optimum system")

    plt.xticks(list(range(0, 25, 4)))
    plt.yticks(list(range(-14, 6, 2)))

    plt.axhline(y=0, color="grey", linestyle="--")
    plt.axvline(x=13, color="grey", linestyle="--")

    plt.xlim(0, 24)

    plt.legend(title="Minimum service when sized")

    plt.savefig(
        f"additional_service_hours_for_all_hours_reattempt_{INDEX}.pdf",
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
        transparent=True,
    )

    # plt.show()

    ax = _heatmap_plot(
        stranded_assets_data[stranded_assets_data["grid_hours"] >= 13],
        "additional_service_hours",
        center=0,
        cmap_name="PiYG",
        default_value=0,
        prioritise_self_generation=False,
        ylabel="Change in average service hours provided",
    )
    _heatmap_plot(
        stranded_assets_data[stranded_assets_data["grid_hours"] >= 13][
            stranded_assets_data["service_hours"] >= 13
        ],
        "service_hours_provided",
        annot=True,
        annot_kws={"color": "b"},
        ax=ax,
        center=19,
        cmap_name="Blues",
        default_value=0,
        figsize=(48 / 5, 32 / 5),
        prioritise_self_generation=False,
        # ylabel="Change in average service hours provided",
    )

    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    plt.savefig(
        f"additional_service_hours_of_stranded_assets_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    ax = _heatmap_plot(
        stranded_assets_data[stranded_assets_data["grid_hours"] >= 13],
        "additional_service_hours",
        center=0,
        cmap_name="PiYG",
        default_value=0,
        prioritise_self_generation=False,
        ylabel="Change in average service hours provided",
    )
    _heatmap_plot(
        stranded_assets_data[stranded_assets_data["grid_hours"] >= 13][
            stranded_assets_data["service_hours"] >= 13
        ],
        "service_hours_provided",
        annot=True,
        annot_kws={"color": "b"},
        ax=ax,
        center=19,
        cmap_name="Blues_r",
        default_value=0,
        figsize=(48 / 5, 32 / 5),
        prioritise_self_generation=False,
        # ylabel="Change in average service hours provided",
    )

    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    plt.savefig(
        f"additional_service_hours_of_stranded_assets_1_reversed_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    #######################################################################
    # Figure 4. The change in LCUE of systems undergoing asset stranding. #
    #######################################################################

    stranded_assets_data["grid_hours"] = stranded_assets_data["grid_hours"].astype(int)
    stranded_assets_data["service_hours"] = stranded_assets_data[
        "service_hours"
    ].astype(int)
    stranded_assets_data["service_hours_provided"] = stranded_assets_data[
        "service_hours_provided"
    ].astype(int)

    stranded_assets_data["change_in_lcue"] = [
        float(f"{entry:.2g}") if float(f"{entry:.2g}") != int(entry) else int(entry)
        for entry in stranded_assets_data["change_in_lcue"].fillna(0)
    ]

    ax = _heatmap_plot(
        (
            _plotting_data := stranded_assets_data[
                stranded_assets_data["grid_hours"]
                < stranded_assets_data["service_hours"]
            ]
        ),
        "change_in_lcue",
        center=0,
        cmap_name="PiYG_r",
        default_value=0,
        figsize=(80 / 5, 32 / 5),
        prioritise_self_generation=False,
        ylabel="Excess LCUE compared to optimisation / USD/kWh",
    )
    _heatmap_plot(
        _plotting_data,
        "change_in_lcue",
        annot=True,
        ax=ax,
        center=0,
        cmap_name="PiYG_r",
        default_value=0,
        figsize=(80 / 5, 32 / 5),
        fmt=".2f",
        prioritise_self_generation=False,
        ylabel="Excess LCUE compared to optimisation / USD/kWh",
    )

    ax.add_patch(
        patches.Rectangle((13.1, 0.05), 0.8, 10.9, edgecolor="grey", fill=False, lw=2)
    )

    ax.add_patch(
        plt.Polygon(
            (
                stranded_assets_columns := [
                    [23.95, 10.95],
                    [23.95, 10.05],
                    [22.95, 10.05],
                    [22.95, 9.05],
                    [21.95, 9.05],
                    [21.95, 8.05],
                    [20.95, 8.05],
                    [20.95, 7.05],
                    [19.95, 7.05],
                    [19.95, 6.05],
                    [18.95, 6.05],
                    [18.95, 5.05],
                    [17.95, 5.05],
                    [17.95, 4.05],
                    [15.05, 4.05],
                    [15.05, 5.05],
                    [14.05, 5.05],
                    [14.05, 10.95],
                ]
            ),
            # color="grey",
            color="#8e0152",
            fill=False,
            lw=2,
            label="Assets stranded by the arrival of the grid",
            # hatch="\\\\//",
            # alpha=0.4,
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                not_simulated_columns := [
                    [24.95, 10.95],
                    [24.05, 10.95],
                    [24.05, 9.95],
                    [23.05, 9.95],
                    [23.05, 8.95],
                    [22.05, 8.95],
                    [22.05, 7.95],
                    [21.05, 7.95],
                    [21.05, 6.95],
                    [20.05, 6.95],
                    [20.05, 5.95],
                    [19.05, 5.95],
                    [19.05, 4.95],
                    [18.05, 4.95],
                    [18.05, 3.95],
                    [17.05, 3.95],
                    [17.05, 2.95],
                    [16.05, 2.95],
                    [16.05, 1.95],
                    [15.05, 1.95],
                    [15.05, 0.95],
                    [14.05, 0.95],
                    [14.05, 0.05],
                    [24.95, 0.05],
                ]
            ),
            color="grey",
            # color="#8e0152",
            fill=True,
            lw=2,
            label="Assets stranded by the arrival of the grid",
            hatch="\\\\//",
            alpha=0.2,
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                grid_supplemented_minigrid_columns := [
                    [14.05, 4.95],
                    [14.95, 4.95],
                    [14.95, 3.95],
                    [16.95, 3.95],
                    [16.95, 3.05],
                    [15.95, 3.05],
                    [15.95, 2.05],
                    [14.95, 2.05],
                    [14.95, 1.05],
                    [14.05, 1.05],
                ]
            ),
            # color="#7f3b08",
            # color="#E04606",
            color="#276419",
            fill=False,
            lw=2,
            label="Increase in grid reliability supplements the minigrid",
            # hatch="\\\\",
            # alpha=0.4,
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                greater_utilisation_columns := [
                    [0.05, 10.95],
                    [12.95, 10.95],
                    # [12.95, 10.05],
                    # [9.95, 10.05],
                    # [9.95, 8.95],
                    # [12.95, 8.95],
                    [12.95, 1.05],
                    [5.05, 1.05],
                    [5.05, 1.95],
                    [6.05, 1.95],
                    [6.05, 3.05],
                    [2.05, 3.05],
                    [2.05, 4.05],
                    [0.05, 4.05],
                    [0.05, 5.95],
                    [1.05, 5.95],
                    [1.05, 7.05],
                    [0.05, 7.05],
                ]
            ),
            # color="#2d004b",
            color="#276419",
            fill=False,
            lw=2,
            label="Decrease in grid reliability results in greater asset utilisation",
            # hatch="\\/",
            # alpha=0.4,
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                undersized_assets_columns := [
                    [12.95, 0.05],
                    [12.95, 0.95],
                    [4.95, 0.95],
                    [4.95, 2.05],
                    [5.95, 2.05],
                    [5.95, 2.95],
                    [1.95, 2.95],
                    [1.95, 3.95],
                    [0.05, 3.95],
                    [0.05, 0.05],
                ]
            ),
            # color="#276419",
            color="#8e0152",
            # color="#e5a618",
            fill=False,
            lw=2,
            label="Decrease in grid reliability results in fewer hours of electricity met so higher LCUE",
            # hatch="\\/..",
            # alpha=0.4,
        )
    )
    ax.add_patch(
        plt.Polygon(
            (
                undersized_assets_columns := [
                    [0.05, 6.05],
                    [0.05, 6.95],
                    [0.95, 6.95],
                    [0.95, 6.05],
                ]
            ),
            color="#8e0152",
            # color="#276419",
            # color="#e5a618",
            fill=False,
            lw=2,
            label="Decrease in grid reliability results in fewer hours of electricity met so higher LCUE",
            # hatch="\\/..",
            # alpha=0.4,
        )
    )

    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    plt.savefig(
        f"change_in_lcue_asset_stranding_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    ###########################################################################
    # Additional figure — Change in emissions intensity when assets stranded. #
    ###########################################################################

    stranded_assets_data["grid_hours"] = stranded_assets_data["grid_hours"].astype(int)
    stranded_assets_data["service_hours"] = stranded_assets_data[
        "service_hours"
    ].astype(int)
    stranded_assets_data["service_hours_provided"] = stranded_assets_data[
        "service_hours_provided"
    ].astype(int)

    stranded_assets_data["change_in_emissions_intensity"] = [
        (
            float(f"{entry:.2g}")
            if float(f"{entry:.2g}") != int(float(f"{entry:.2g}"))
            else int(float(f"{entry:.2g}"))
        )
        for entry in stranded_assets_data["change_in_emissions_intensity"].fillna(0)
    ]

    ax = _heatmap_plot(
        stranded_assets_data,
        "change_in_emissions_intensity",
        center=0,
        cmap_name="PuOr_r",
        default_value=0,
        figsize=(80 / 5, 32 / 5),
        prioritise_self_generation=False,
        ylabel="Change in emissions intensity / kgCO$_\mathrm{2}$eq/kWh",
    )
    _heatmap_plot(
        stranded_assets_data,
        "change_in_emissions_intensity",
        annot=True,
        ax=ax,
        center=0,
        cmap_name="PuOr_r",
        default_value=0,
        fmt=".3g",
        prioritise_self_generation=False,
        ylabel="Change in emissions intensity / kgCO$_\mathrm{2}$eq/kWh",
    )

    ax.add_patch(
        patches.Rectangle((13.1, 0.05), 0.8, 10.9, edgecolor="C3", fill=False, lw=2)
    )

    ax.add_patch(
        plt.Polygon(
            (
                stranded_emissions := [
                    [14.05, 10.95],
                    [24.95, 10.95],
                    [24.95, 0.05],
                    [23.05, 0.05],
                    [23.05, 3.05],
                    [22.05, 3.05],
                    [22.05, 5.05],
                    [21.05, 5.05],
                    [21.05, 7.05],
                    [19.95, 7.05],
                    [19.95, 6.05],
                    [18.95, 6.05],
                    [18.95, 5.05],
                    [14.05, 5.05],
                    [14.05, 10.95],
                ]
            ),
            color="#b35806",
            fill=False,
            lw=2,
            label="Increase in grid availability results in asset utilisation decreasing, "
            "and, hence, an effective increase in emissions intensity.",
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                greater_static_emissions_impact_columns := [
                    [3.05, 0.05],
                    [3.05, 7.95],
                    [12.95, 7.95],
                    [12.95, 0.05],
                ]
            ),
            color="#b35806",
            fill=False,
            lw=2,
            label="The lower provision of service hours of the smaller system results "
            "in fewer hours of electricity so a greater impact resulting from the static "
            "emissions.",
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                introduction_of_renewables_to_grid_columns := [
                    [14.05, 0.05],
                    [22.95, 0.05],
                    [22.95, 2.95],
                    [21.95, 2.95],
                    [21.95, 4.95],
                    [20.95, 4.95],
                    [20.95, 6.95],
                    [20.05, 6.95],
                    [20.05, 5.95],
                    [19.05, 5.95],
                    [19.05, 4.95],
                    [14.05, 4.95],
                ]
            ),
            color="#542788",
            fill=False,
            lw=2,
            label="Reneable PV-storage assets provide renewable electricity alongside the "
            "grid for non-grid hours, reducing the effective emissions intensity.",
        )
    )

    ax.add_patch(
        plt.Polygon(
            (
                greater_utilisation_of_installed_assets := [
                    [0.05, 0.05],
                    [0.05, 10.95],
                    [12.95, 10.95],
                    [12.95, 8.05],
                    [2.95, 8.05],
                    [2.95, 0.05],
                ]
            ),
            color="#542788",
            fill=False,
            lw=2,
            label="Smaller systems not sized for very-low grid availability or very-high "
            "thresholds of minimum service are more greatly utilised than optimised assets "
            "and hence better offset the embedded static emissions.",
        )
    )

    plt.xlabel("Hours of grid availability per day")
    plt.ylabel("Minimum average hours of service per day")
    plt.savefig(
        f"change_in_emissions_asset_stranding_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )
    # plt.show()

    ##################################################################
    # Additional figure — Change in unmet energy of stranded assets. #
    ##################################################################

    _heatmap_plot(
        stranded_assets_data,
        "change_in_unmet_energy_fraction",
        center=0,
        cmap_name="PiYG_r",
        default_value=0,
        prioritise_self_generation=False,
        ylabel="Change in unmet energy fraction",
    )

    _heatmap_plot(
        subsidised_data,
        "unmet_energy_fraction",
        # cmap_colour="#006264",
        cmap_name="Oranges",
        default_value=0,
        vmin=0,
        vmax=1,
        ylabel="Minimum average hours of service per day",
    )

    # The y limits for the plot should be scaled differently based on the range of
    # utilisations that are likely to be plotted.
    # plt_ylim: float = 0.8

    # Generate plots for both types of grid prioritisation
    for grid_prioritisation, plt_ylim, initial_hours_string in zip(
        ("GP", "SG", "GP", "SG"),
        (1.05, 1.05, 1.05, 1.05),
        ("_from_24", "_from_24", "_from_16", "_from_16"),
    ):
        # for grid_prioritisation, plt_ylim, initial_hours_string in zip(
        #     ("GP", "SG"), (1.05, 1.05), ("_from_16", "_from_16")
        # ):

        ####################################################
        # Figure 5. Change in asset utilisation over time. #
        ####################################################

        import pdb

        pdb.set_trace()

        initial: int = 13
        data_map: dict[int, pd.DataFrame] = {}
        info_file_map: dict[Any, Any] = {}
        for final, percent, less_more in zip(
            (15, 18, 21, 24, 13, 13),
            (20, 50, 80, 100, 50, 0),
            ("less", "less", "less", "less", "more", None),
        ):
            change_and_percent = (
                f"{percent}_{less_more}_demand" if percent != 0 else "sim"
            )
            with open(
                os.path.join(
                    os.path.expanduser("~"),
                    "clover_locations",
                    "UP",
                    "outputs",
                    "simulation_outputs",
                    f"{initial}_to_{final}_{change_and_percent}{initial_hours_string}_{grid_prioritisation}",
                    "simulation_output_1.csv",
                ),
                "r",
            ) as simulation_results_file:
                data_map[(final, percent, less_more)] = pd.read_csv(
                    simulation_results_file
                )
            with open(
                os.path.join(
                    os.path.expanduser("~"),
                    "clover_locations",
                    "UP",
                    "outputs",
                    "simulation_outputs",
                    f"{initial}_to_{final}_{change_and_percent}{initial_hours_string}_{grid_prioritisation}",
                    "info_file.json",
                ),
                "r",
            ) as info_file:
                info_file_map[(final, percent, less_more)] = json.load(info_file)

        # Compute the utilisations as daily averages
        utilisation_map: dict[int, dict[int, float] | pd.DataFrame] = (
            collections.defaultdict(dict)
        )
        for key, result in data_map.items():
            for day in range(int(len(result) / 24)):
                utilisation_map[key][day] = (
                    sum(
                        result["Renewables energy used (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                    + sum(
                        result["Storage energy supplied (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                ) / (
                    sum(result["PV energy supplied (kWh)"][24 * day : 24 * (day + 1)])
                    + sum(result["Dumped energy (kWh)"][24 * day : 24 * (day + 1)])
                )
            utilisation_map[key] = pd.DataFrame.from_dict(
                utilisation_map[key], orient="index"
            )

        # Open the standard simulation
        # # with open(
        # #     os.path.join(
        # #         os.path.expanduser("~"),
        # #         "clover_locations",
        # #         "UP",
        # #         "outputs",
        # #         "simulation_outputs",
        # #         "13_to_13_sim",
        # #         "simulation_output_1.csv",
        # #     ),
        # #     "r",
        # # ) as simulation_results_file:
        # #     grid_13 = pd.read_csv(simulation_results_file)
        # #     grid_13["utilisation"] = grid_13["Renewables energy supplied (kWh)"] / (
        # #         grid_13["PV energy supplied (kWh)"] + grid_13["Dumped energy (kWh)"]
        # #     ).astype(float)

        # # # Compute its daily utilisation
        # # for day in range(int(len(grid_13) / 24)):
        # #     utilisation_map[13][day] = (
        # #         sum(grid_13["Renewables energy used (kWh)"][24 * day : 24 * (day + 1)])
        # #         + sum(grid_13["Storage energy supplied (kWh)"][24 * day : 24 * (day + 1)])
        # #     ) / (
        # #         sum(grid_13["PV energy supplied (kWh)"][24 * day : 24 * (day + 1)])
        # #         + sum(grid_13["Dumped energy (kWh)"][24 * day : 24 * (day + 1)])
        # #     )

        # utilisation_map[13] = pd.DataFrame.from_dict(utilisation_map[13], orient="index")

        # Plot a rolling average of the utilisations
        sns.set_palette(
            sns.color_palette(
                ["#E5D4E8", "#C3A4CF", "#986EAC", "#742881", "#5CAE63", "#666666"]
            )
        )
        sns.set_context("notebook")
        sns.set_style("ticks")

        fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))

        for index, result in enumerate(utilisation_map.values()):
            plt.plot(
                (result.index / (8760 / 24)),
                result[0].rolling(window=(_window := 30)).mean(),
                "--" if index == 5 else "-",
                color=f"C{index}",
            )

        plt.xlabel("Year of simulation")
        plt.ylabel("Mini-grid asset utilisation")

        plt.xlim(0, 8)
        plt.ylim(0, plt_ylim)

        plt.tick_params(right=True, top=True, labelrotation=0)

        plt.savefig(
            f"utilisation_changing_grid_reliability{initial_hours_string}_{grid_prioritisation}_{INDEX}.pdf",
            dpi=DPI,
            bbox_inches="tight",
            pad_inches=0,
            transparent=True,
        )

        plt.savefig(
            f"utilisation_changing_grid_reliability{initial_hours_string}_{grid_prioritisation}_{INDEX}.png",
            dpi=DPI,
            bbox_inches="tight",
            pad_inches=0,
            transparent=True,
        )

        ###################################################################
        # Thesis Figure. Change in asset utilisation with fixed user base #
        ###################################################################

        initial: int = 13
        data_map: dict[int, pd.DataFrame] = {}
        for final, percent, less_more in zip(
            (13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24), [0] * 12, [None] * 12
        ):
            change_and_percent = (
                f"{percent}_percent_{less_more}_demand" if percent != 0 else "sim"
            )
            with open(
                os.path.join(
                    os.path.expanduser("~"),
                    "clover_locations",
                    "UP",
                    "outputs",
                    "simulation_outputs",
                    f"{initial}_to_{final}_{change_and_percent}{initial_hours_string}_{grid_prioritisation}",
                    "simulation_output_1.csv",
                ),
                "r",
            ) as simulation_results_file:
                data_map[(final, percent, less_more)] = pd.read_csv(
                    simulation_results_file
                )

        # Compute the utilisations as daily averages
        utilisation_map: dict[int, dict[int, float] | pd.DataFrame] = (
            collections.defaultdict(dict)
        )
        for key, result in data_map.items():
            for day in range(int(len(result) / 24)):
                utilisation_map[key][day] = (
                    sum(
                        result["Renewables energy used (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                    + sum(
                        result["Storage energy supplied (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                ) / (
                    sum(result["PV energy supplied (kWh)"][24 * day : 24 * (day + 1)])
                    + sum(result["Dumped energy (kWh)"][24 * day : 24 * (day + 1)])
                )

            utilisation_map[key] = pd.DataFrame.from_dict(
                utilisation_map[key], orient="index"
            )

        for key, result in data_map.items():
            result["grid_fraction"] = (
                result["Grid energy (kWh)"] / result["Total energy used (kWh)"]
            )

        # Plot a rolling average of the utilisations
        sns.set_palette(
            (palette := sns.cubehelix_palette(start=0.6, rot=-0.8, n_colors=12))
        )
        sns.set_context("notebook")
        sns.set_style("ticks")

        fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))

        for index, result in enumerate(utilisation_map.values()):
            plt.plot(
                (result.index / (8760 / 24)),
                result[0].rolling(window=(_window := 30)).mean(),
                # "--" if index == 5 else "-",
                color=f"C{index}",
            )
            # plt.fill_between(
            #     (result.index / (8760 / 24)),
            #     result[0].rolling(window=(_window := 30)).quantile(0.25),
            #     result[0].rolling(window=_window).quantile(0.75),
            #     alpha=0.05,
            #     color=f"C{index}",
            # )
            # plt.plot(
            #     (result:=data_map[list(data_map.keys())[index]]["grid_fraction"].dropna()).index / (8760),
            #     (result).rolling(window=(_window:=300)).mean(),
            #     # "--" if index == 5 else "-",
            #     "--",
            #     color=f"C{index}",
            # )

        plt.xlabel("Year of simulation")
        plt.ylabel("Mini-grid asset utilisation")

        plt.xlim(0, 8)
        plt.ylim(0.05, plt_ylim)

        plt.tick_params(right=True, top=True, labelrotation=0)

        norm = plt.Normalize(12.5, 24.5)
        scalar_mappable = plt.cm.ScalarMappable(
            cmap=mcolors.LinearSegmentedColormap.from_list(
                "Custom", sns.color_palette().as_hex(), 12
            ),
            norm=norm,
        )

        colorbar = axis.figure.colorbar(
            scalar_mappable, ax=ax, label="Average final grid hours per day / hours"
        )
        colorbar.set_ticks(list(range(13, 25))[::2])

        plt.savefig(
            f"utilisation_changing_grid_reliability_no_userbase_change{initial_hours_string}_{grid_prioritisation}_{INDEX}.pdf",
            dpi=DPI,
            bbox_inches="tight",
            pad_inches=0,
            transparent=True,
        )

        fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))

        for index, result in enumerate(utilisation_map.values()):
            # plt.plot(
            #     (result.index / (8760 / 24)),
            #     result[0].rolling(window=(_window := 30)).mean(),
            #     # "--" if index == 5 else "-",
            #     color=f"C{index}",
            # )
            # plt.fill_between(
            #     (result.index / (8760 / 24)),
            #     result[0].rolling(window=(_window := 30)).quantile(0.25),
            #     result[0].rolling(window=_window).quantile(0.75),
            #     alpha=0.05,
            #     color=f"C{index}",
            # )
            plt.plot(
                (
                    result := data_map[list(data_map.keys())[index]][
                        "grid_fraction"
                    ].dropna()
                ).index
                / (8760),
                (result).rolling(window=(_window := 300)).mean(),
                # "--" if index == 5 else "-",
                "--",
                color=f"C{index}",
            )

        plt.xlabel("Year of simulation")
        plt.ylabel("Fraction of electricity sourced from the national grid")

        plt.xlim(0, 8)
        plt.ylim(0.05, plt_ylim)

        plt.tick_params(right=True, top=True, labelrotation=0)

        norm = plt.Normalize(12.5, 24.5)
        scalar_mappable = plt.cm.ScalarMappable(
            cmap=mcolors.LinearSegmentedColormap.from_list(
                "Custom", sns.color_palette().as_hex(), 12
            ),
            norm=norm,
        )

        colorbar = axis.figure.colorbar(
            scalar_mappable, ax=ax, label="Average final grid hours per day / hours"
        )
        colorbar.set_ticks(list(range(13, 25))[::2])

        plt.savefig(
            f"grid_fraction_changing_grid_reliability_no_userbase_change{initial_hours_string}_{grid_prioritisation}_{INDEX}.pdf",
            dpi=DPI,
            bbox_inches="tight",
            pad_inches=0,
            transparent=True,
        )

        ##################################################################
        # Thesis Figure. Change in user base with fixed grid reliability #
        ##################################################################

        initial: int = 13
        data_map: dict[int, pd.DataFrame] = {}
        user_base_change_range = list(range(-100, 110, 10))
        for final, percent, less_more in zip(
            [13] * len(user_base_change_range),
            [abs(entry) for entry in user_base_change_range],
            ["less"] * int((len(user_base_change_range) - 1) / 2)
            + [None]
            + ["more"] * int((len(user_base_change_range) - 1) / 2),
        ):
            change_and_percent = (
                f"{percent}_percent_{less_more}_demand" if percent != 0 else "sim"
            )
            with open(
                os.path.join(
                    os.path.expanduser("~"),
                    "clover_locations",
                    "UP",
                    "outputs",
                    "simulation_outputs",
                    f"{initial}_to_{final}_{change_and_percent}{initial_hours_string}_{grid_prioritisation}",
                    "simulation_output_1.csv",
                ),
                "r",
            ) as simulation_results_file:
                data_map[(final, percent, less_more)] = pd.read_csv(
                    simulation_results_file
                )

        # Compute the utilisations as daily averages
        utilisation_map: dict[int, dict[int, float] | pd.DataFrame] = (
            collections.defaultdict(dict)
        )
        for key, result in data_map.items():
            for day in range(int(len(result) / 24)):
                utilisation_map[key][day] = (
                    sum(
                        result["Renewables energy used (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                    + sum(
                        result["Storage energy supplied (kWh)"][
                            24 * day : 24 * (day + 1)
                        ]
                    )
                ) / (
                    sum(result["PV energy supplied (kWh)"][24 * day : 24 * (day + 1)])
                    + sum(result["Dumped energy (kWh)"][24 * day : 24 * (day + 1)])
                )
            utilisation_map[key] = pd.DataFrame.from_dict(
                utilisation_map[key], orient="index"
            )

        # Plot a rolling average of the utilisations
        colourmap = mcolors.LinearSegmentedColormap.from_list(
            "", np.array(PURPLE_GREEN_COLOURS), 40
        )
        sns.set_palette(
            # sns.color_palette(
            #     colourmap(np.linspace(0, 1, len(user_base_change_range))).tolist()
            # )
            "BrBG",
            n_colors=21,
        )
        sns.set_context("notebook")
        sns.set_style("ticks")

        fig, ax = plt.subplots(figsize=(48 / 5, 32 / 5))

        for index, result in enumerate(utilisation_map.values()):
            plt.plot(
                (result.index / (8760 / 24)),
                result[0].rolling(window=(_window := 30)).mean(),
                # "--" if index == 5 else "-",
                color=f"C{index}",
            )

        plt.xlabel("Year of simulation")
        plt.ylabel("Mini-grid asset utilisation")

        plt.xlim(0, 8)
        plt.ylim(0, plt_ylim)

        plt.tick_params(right=True, top=True, labelrotation=0)

        norm = plt.Normalize(-105, 105)
        scalar_mappable = plt.cm.ScalarMappable(
            cmap=mcolors.LinearSegmentedColormap.from_list(
                "Custom", sns.color_palette().as_hex(), len(user_base_change_range)
            ),
            norm=norm,
        )

        colorbar = axis.figure.colorbar(
            scalar_mappable, ax=ax, label="Percentage user-base change / %"
        )
        colorbar.set_ticks(list(range(-100, 110, 10))[::1])

        plt.savefig(
            f"utilisation_changing_userbase_for_{initial}_hour_grid{initial_hours_string}_{grid_prioritisation}_{INDEX}.pdf",
            dpi=DPI,
            bbox_inches="tight",
            pad_inches=0,
            transparent=True,
        )

    # FIGURE #

    fig, axes = plt.subplots(1, 24)

    palette = sns.color_palette(
        [
            "#230017",
            "#560133",
            "#790149",
            "#9F0162",
            "#C7007C",
            "#EF0096",
            "#FF5AAF",
            "#FF9DC8",
            "#FFCFE2",
            "#00306F",
            "#00489E",
            "#005FCC",
            "#0079DA",
            "#009FFA",
            "#00C2F9",
            "#00E5F8",
            "#7CFFFA",
            "#004002",
            "#005A01",
            "#007702",
            "#009503",
            "#00B408",
            "#00D302",
            "#00F407",
            "#AFFF2A",
        ]
    )

    sns.set_palette((palette := sns.cubehelix_palette(rot=-0.2)))

    ################################################################################
    # Alternative Figure 2. Impact of self- vs grid-prioritisiation as a lineplot. #
    ################################################################################

    fig = plt.figure(figsize=(8, 4.5))
    sns.stripplot(
        subsidised_data[
            subsidised_data["prioritise_self_generation"] == False
        ].dropna(),
        x="grid_hours",
        y="grid_to_self_lcue_percentage_change",
        hue="service_hours_provided",
        marker="D",
        alpha=0.8,
        palette=palette,
    )
    # plt.show()

    line_plot_data: dict[int, Any] = {}
    sns.set_style("ticks")

    fig = plt.figure(figsize=(48 / 5, 32 / 5))
    sns.set_palette(palette := sns.cubehelix_palette(rot=-0.4, n_colors=24))

    for hour in reversed(range(1, 25)):
        line_plot_data[hour] = subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
            & (subsidised_data["service_hours"] == hour)
        ]
        sns.lineplot(
            data=line_plot_data[hour],
            x="grid_hours",
            y="grid_to_self_lcue_percentage_change",
            label=hour,
            color=f"C{24-hour}",
        )

    plt.ylabel("Percentage change in LCUE when\nswitching to self-prioritisation")
    plt.xlabel("Hours of national grid availability")
    plt.xlim(0, 23)
    (axis := plt.gca()).get_legend().remove()
    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)
    plt.savefig(
        f"self_prioritsation_lcue_percentage_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    ############################################################################
    # Alternative Figure 2. Emissions impact of self- vs grid-prioritisiation. #
    ############################################################################

    fig = plt.figure(figsize=(48 / 5, 32 / 5))
    sns.set_palette(palette := sns.cubehelix_palette(rot=-0.8, n_colors=24))
    for hour in reversed(range(1, 25)):
        sns.lineplot(
            data=line_plot_data[hour],
            x="grid_hours",
            y="grid_to_self_emissions_intensity_percentage_change",
            label=hour,
            color=f"C{24-hour}",
        )

    plt.ylabel(
        "Percentage change in emissions intensity\nwhen switching to self-prioritisation"
    )
    plt.xlabel("Hours of national grid availability")
    plt.xlim(0, 23)
    (axis := plt.gca()).get_legend().remove()
    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)
    plt.savefig(
        f"self_prioritsation_emissions_intensity_percentage_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    fig = plt.figure(figsize=(48 / 5, 32 / 5))
    sns.set_palette(palette := sns.cubehelix_palette(rot=-1.2, n_colors=24))
    for hour in reversed(range(1, 25)):
        sns.lineplot(
            data=line_plot_data[hour],
            x="grid_hours",
            y="grid_to_self_unmet_energy_fraction_change",
            label=hour,
            color=f"C{24-hour}",
        )

    plt.ylabel(
        "Absolute change in unmet-energy fraction\nwhen switching to self-prioritisation"
    )
    plt.xlabel("Hours of national grid availability")
    plt.xlim(0, 23)
    (axis := plt.gca()).get_legend().remove()
    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)
    plt.savefig(
        f"self_prioritsation_unmet_energy_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    ######################################################
    # Figure 2. Impact of self- vs grid-prioritisiation. #
    ######################################################

    fig, axes = plt.subplots(3, 1, figsize=(48 / 5, 64 / 5))
    sns.set_palette(
        palette := sns.cubehelix_palette(rot=-0.4, n_colors=24, reverse=True)
    )

    sns.scatterplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_lcue_percentage_change",
        hue="service_hours",
        palette=palette,
        marker="H",
        alpha=0.4,
        s=250,
        linewidth=0,
        ax=(axis := axes[0]),
    )

    sns.regplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_lcue_percentage_change",
        color="C0",
        order=5,
        scatter=False,
        ax=axis,
    )

    axis.set_ylabel("Percentage change in LCUE when\nswitching to self-prioritisation")
    # axis.set_xlabel("Hours of national grid availability")
    axis.set_xlabel("")
    axis.set_xlim(0, 23)
    axis.get_legend().remove()

    norm = plt.Normalize(0, 24)
    scalar_mappable = plt.cm.ScalarMappable(
        cmap=mcolors.LinearSegmentedColormap.from_list(
            "Custom", sns.color_palette().as_hex(), 24
        ),
        norm=norm,
    )

    colorbar = axis.figure.colorbar(
        scalar_mappable, ax=axis, label="Minimum service hours"
    )
    colorbar.set_ticks(list(range(25))[::4])

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.text(
        (_label_x_position := -0.12),
        (_label_y_position := 1.1),
        "a.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    sns.set_palette(
        palette := sns.cubehelix_palette(rot=-0.8, n_colors=24, reverse=True)
    )

    sns.scatterplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_emissions_intensity_percentage_change",
        hue="service_hours",
        palette=palette,
        marker="H",
        alpha=0.4,
        s=250,
        linewidth=0,
        ax=(axis := axes[1]),
    )

    sns.regplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_emissions_intensity_percentage_change",
        color="C0",
        order=5,
        scatter=False,
        ax=axis,
    )

    axis.set_ylabel(
        "Percentage change in emissions\nintensity when switching to self-prioritisation"
    )
    # axis.set_xlabel("Hours of national grid availability")
    axis.set_xlabel("")
    axis.set_xlim(0, 23)
    axis.get_legend().remove()

    norm = plt.Normalize(0, 24)
    scalar_mappable = plt.cm.ScalarMappable(
        cmap=mcolors.LinearSegmentedColormap.from_list(
            "Custom", sns.color_palette().as_hex(), 24
        ),
        norm=norm,
    )

    colorbar = axis.figure.colorbar(
        scalar_mappable, ax=axis, label="Minimum service hours"
    )
    colorbar.set_ticks(list(range(25))[::4])

    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.text(
        _label_x_position,
        _label_y_position,
        "b.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    sns.set_palette(palette := sns.cubehelix_palette(rot=-1.2, n_colors=24))

    sns.scatterplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_unmet_energy_fraction_change",
        hue="service_hours",
        palette=palette,
        marker="H",
        alpha=0.4,
        s=250,
        linewidth=0,
        ax=(axis := axes[2]),
    )

    sns.regplot(
        data=subsidised_data.dropna()[
            (subsidised_data["prioritise_self_generation"] == False)
        ],
        x="grid_hours",
        y="grid_to_self_unmet_energy_fraction_change",
        color="C23",
        order=5,
        scatter=False,
        ax=axis,
    )

    axis.set_ylabel(
        "Absolute change in unmet-energy\nfraction when switching to self-prioritisation"
    )
    axis.set_xlabel("Hours of national grid availability")
    axis.get_legend().remove()

    norm = plt.Normalize(0, 24)
    scalar_mappable = plt.cm.ScalarMappable(
        cmap=mcolors.LinearSegmentedColormap.from_list(
            "Custom", sns.color_palette().as_hex(), 24
        ),
        norm=norm,
    )

    colorbar = axis.figure.colorbar(
        scalar_mappable, ax=axis, label="Minimum service hours"
    )
    colorbar.set_ticks(list(range(25))[::4])

    axis.set_xlim(0, 23)
    plt.setp(axis.get_xticklabels(), visible=True)
    plt.setp(axis.get_yticklabels(), visible=True)
    plt.setp(axis.get_xticklines(), visible=True)
    plt.setp(axis.get_yticklines(), visible=True)
    plt.setp(axis.spines.values(), visible=True)

    axis.text(
        _label_x_position,
        _label_y_position,
        "c.",
        transform=axis.transAxes,
        fontsize=(_fontsize := 32),
        fontweight="bold",
        va="top",
        ha="right",
    )

    plt.savefig(
        f"self_prioritsation_comparison_{INDEX}.pdf",
        transparent=True,
        dpi=DPI,
        bbox_inches="tight",
        pad_inches=0,
    )

    pdb.set_trace(header=f"Press 'q' to exit.")


if __name__ == "__main__":
    main(sys.argv[1:])