Postprocessing¶

def get_impacts(self) -> pd.DataFrame:
    """
    Extract environmental impacts by category, process, and time.

    Returns denormalized impact values from the solved optimization model,
    organized as a pivoted DataFrame with time as rows and (category, process)
    as column MultiIndex.

    Returns
    -------
    pd.DataFrame
        Pivoted DataFrame with 'Time' as index and MultiIndex columns for
        (Category, Process) combinations. Values represent environmental
        impacts in the units of the characterization method.
    """
    if hasattr(self, "df_impacts"):
        return self.df_impacts

    impacts = {}
    cat_scales = getattr(self.m, "scales", {}).get("characterization", 1.0)
    fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0)
    impacts = {
        (c, p, t): pyo.value(self.m.specific_impact[c, p, t])
        * cat_scales[c]
        * fg_scale  # Unscale foreground impacts
        for c in self.m.CATEGORY
        for p in self.m.PROCESS
        for t in self.m.SYSTEM_TIME
    }
    df = pd.DataFrame.from_dict(impacts, orient="index", columns=["Value"])
    df.index = pd.MultiIndex.from_tuples(
        df.index, names=["Category", "Process", "Time"]
    )
    df = df.reset_index()
    df_pivot = df.pivot(
        index="Time", columns=["Category", "Process"], values="Value"
    )
    self.df_impacts = df_pivot
    return self.df_impacts

def get_dynamic_inventory(self, biosphere_database: str = "ecoinvent-3.12-biosphere") -> pd.DataFrame:
    """
    Extract the dynamic inventory from the solved model.

    Returns a DataFrame with elementary flows over time, formatted for use
    with dynamic_characterization.

    Parameters
    ----------
    biosphere_database : str, default="ecoinvent-3.12-biosphere"
        Name of the biosphere database to look up flow IDs.

    Returns
    -------
    pd.DataFrame
        DataFrame with columns: activity, flow, date, amount.
        - activity: process code (str)
        - flow: biosphere flow ID (int)
        - date: datetime of emission
        - amount: flow amount (float)
    """
    fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0)

    # The inventory is assembled from the solved flow values and the background
    # inventory data rather than by evaluating the `scaled_inventory` expression
    # once per (process, elementary flow, year): the latter walks a deep
    # expression tree hundreds of thousands of times.
    mapped_inventory = getattr(self.m, "_mapped_inventory_by_flow", None)
    inventory = {}
    if mapped_inventory is not None:
        for (p, i, t) in self.m.scaled_technosphere_flow:
            amount = pyo.value(self.m.scaled_technosphere_flow[p, i, t])
            if not amount:
                continue
            for e, factor in mapped_inventory.get((i, t), ()):
                key = (p, e, t)
                inventory[key] = inventory.get(key, 0.0) + amount * factor
        for (p, e, t) in self.m.scaled_biosphere_flow:
            amount = pyo.value(self.m.scaled_biosphere_flow[p, e, t])
            if amount:
                key = (p, e, t)
                inventory[key] = inventory.get(key, 0.0) + amount
    else:  # pragma: no cover - models built by earlier versions
        inventory = {
            (p, e, t): pyo.value(self.m.scaled_inventory[p, e, t])
            for p in self.m.PROCESS
            for e in self.m.ELEMENTARY_FLOW
            for t in self.m.SYSTEM_TIME
        }

    df = pd.DataFrame.from_records(
        [
            (p, e, t, value * fg_scale)
            for (p, e, t), value in inventory.items()
            if value
        ],
        columns=["activity", "flow", "date", "amount"]
    ).astype({
        "activity": "str",
        "flow": "str",
        "amount": "float64"
    })

    # Convert year integers to datetime
    df["date"] = pd.to_datetime(df["date"].astype(int), format="%Y")

    # Convert flow codes to database IDs, from a single pass over the biosphere
    # database instead of one query per row.
    biosphere_ids = {
        node["code"]: node["id"] for node in bd.Database(biosphere_database)
    }
    df["flow"] = df["flow"].map(biosphere_ids)

    self.df_dynamic_inventory = df
    return self.df_dynamic_inventory

def get_characterized_dynamic_inventory(
    self,
    base_lcia_method: tuple,
    metric: str = "radiative_forcing",
    time_horizon: int = 100,
    fixed_time_horizon: bool = True,
    biosphere_database: str = "ecoinvent-3.12-biosphere",
    df_inventory: pd.DataFrame = None,
) -> pd.DataFrame:
    """
    Characterize the dynamic inventory using dynamic_characterization.

    Parameters
    ----------
    base_lcia_method : tuple
        The LCIA method tuple for characterization (e.g., ('IPCC', 'GWP100')).
    metric : str, default="radiative_forcing"
        Characterization metric. Options: "radiative_forcing", "GWP".
    time_horizon : int, default=100
        Time horizon for characterization in years.
    fixed_time_horizon : bool, default=True
        If True, use fixed time horizon; if False, use dynamic time horizon.
    biosphere_database : str, default="ecoinvent-3.12-biosphere"
        Name of the biosphere database (used if df_inventory not provided).
    df_inventory : pd.DataFrame, optional
        Pre-computed inventory DataFrame. If not provided, calls get_dynamic_inventory().

    Returns
    -------
    pd.DataFrame
        Characterized inventory DataFrame with columns: date, amount.
    """
    from dynamic_characterization import characterize

    if df_inventory is None:
        df_inventory = self.get_dynamic_inventory(biosphere_database=biosphere_database)

    df_characterized = characterize(
        df_inventory,
        metric=metric,
        base_lcia_method=base_lcia_method,
        time_horizon=time_horizon,
        fixed_time_horizon=fixed_time_horizon,
    )

    self.df_characterized_inventory = df_characterized
    return self.df_characterized_inventory

def plot_characterized_dynamic_inventory(
    self,
    base_lcia_method: tuple = None,
    metric: str = "radiative_forcing",
    time_horizon: int = 100,
    fixed_time_horizon: bool = True,
    biosphere_database: str = "ecoinvent-3.12-biosphere",
    df_characterized: pd.DataFrame = None,
):
    """
    Plot the characterized dynamic inventory aggregated by year.

    Parameters
    ----------
    base_lcia_method : tuple, optional
        The LCIA method tuple for characterization. Required if df_characterized
        is not provided.
    metric : str, default="radiative_forcing"
        Characterization metric (used if df_characterized not provided).
    time_horizon : int, default=100
        Time horizon for characterization (used if df_characterized not provided).
    fixed_time_horizon : bool, default=True
        If True, use fixed time horizon (used if df_characterized not provided).
    biosphere_database : str, default="ecoinvent-3.12-biosphere"
        Name of the biosphere database (used if df_characterized not provided).
    df_characterized : pd.DataFrame, optional
        Pre-computed characterized inventory. If not provided, calls
        get_characterized_dynamic_inventory().
    """
    if df_characterized is None:
        if base_lcia_method is None:
            raise ValueError("base_lcia_method is required when df_characterized is not provided")
        df_characterized = self.get_characterized_dynamic_inventory(
            base_lcia_method=base_lcia_method,
            metric=metric,
            time_horizon=time_horizon,
            fixed_time_horizon=fixed_time_horizon,
            biosphere_database=biosphere_database,
        )

    # Ensure date column is datetime
    df_plot = df_characterized.copy()
    df_plot["date"] = pd.to_datetime(df_plot["date"])

    # Round to nearest year and aggregate
    df_grouped = (
        df_plot
        .assign(date_rounded=(df_plot["date"] + pd.offsets.MonthBegin(6)).dt.to_period("Y").dt.to_timestamp())
        .groupby("date_rounded")["amount"]
        .sum()
        .reset_index()
    )

    # Create plot
    fig, axes = self._create_clean_axes()
    ax = axes[0]

    ax.plot(
        df_grouped["date_rounded"],
        df_grouped["amount"],
        marker=self._plot_config["line_marker"],
        linewidth=self._plot_config["line_width"],
        color=self._plot_config["line_color"],
        label=metric.replace("_", " ").title(),
    )

    ax.set_ylabel(f"{metric.replace('_', ' ').title()}", fontsize=self._plot_config["label_fontsize"])
    ax.set_title(f"Dynamic {metric.replace('_', ' ').title()}", fontsize=self._plot_config["title_fontsize"])
    ax.set_axisbelow(True)
    ax.grid(
        axis="both",
        linestyle=self._plot_config["grid_linestyle"],
        alpha=self._plot_config["grid_alpha"],
    )

    self._add_legend(ax, position="bottom", bbox_to_anchor=(0.5, -0.35))

    fig.tight_layout()
    plt.show()

def get_installation(self) -> pd.DataFrame:
    """
    Extracts the installation data from the model and returns it as a DataFrame.

    Values are the decision variable `var_installation[p, v]` itself, unchanged:
    the number of process UNITS built in vintage year v. One unit delivers its full
    lifetime production (the sum of its production temporal distribution) spread
    over its operation window.

    Two things follow, and both matter when reading plots:

    1. The index is the year of INSTALLATION, not a year of production. A unit
       installed in 2030 with an operation window of tau 1-20 produces in 2031-2050.
    2. The values are a lifetime quantity, so they are not an annual capacity and
       must not be compared with the per-year values from `get_production()` or
       `get_demand()`. Use `get_production_capacity()` for the annual capacity that
       installed units and existing stock make available in each year.

    Returns
    -------
    pd.DataFrame
        Time (vintage year) as index, Process as columns, units as values.
    """
    # var_installation is already in real units, no scaling needed
    installation_matrix = {
        (t, p): pyo.value(self.m.var_installation[p, t])
        for p in self.m.PROCESS
        for t in self.m.SYSTEM_TIME
    }
    df = pd.DataFrame.from_dict(
        installation_matrix, orient="index", columns=["Value"]
    )
    df.index = pd.MultiIndex.from_tuples(df.index, names=["Time", "Process"])
    df = df.reset_index()
    df_pivot = df.pivot(index="Time", columns="Process", values="Value")
    self.df_installation = df_pivot
    return self.df_installation

def get_operation(self, aggregate_vintages: bool = True) -> pd.DataFrame:
    """
    Extracts the operation data from the model and returns it as a DataFrame.

    Values are the decision variable `var_operation[p, v, t]` itself: the number of
    UNITS of vintage v that run in year t, summed over vintages by default. Unlike
    `get_installation()`, the index is the year of OPERATION.

    Units running are not a production volume: multiply by the output per unit and
    year (the production entry at the vintage's lifecycle stage) to get production,
    or simply use `get_production()`. Units running can be compared directly with
    `get_installation()` only per vintage, since operation of a vintage is bounded
    by the units installed in that vintage.

    Parameters
    ----------
    aggregate_vintages : bool, default=True
        If True (default), sum operation across vintages for each (process, time)
        to provide backward-compatible 2D output.
        If False, return full 3D data with (Process, Vintage) as MultiIndex columns.

    Returns
    -------
    pd.DataFrame
        If aggregate_vintages=True: DataFrame with Time as index, Process as columns.
        If aggregate_vintages=False: DataFrame with Time as index,
            (Process, Vintage) MultiIndex columns.
        Values are counts of running units in both cases.

    Note: var_operation is not scaled because when both demand and
    foreground_production are scaled by the same factor, the scaling
    cancels out in the constraint: demand = production * operation.
    """
    if aggregate_vintages:
        # Aggregate across vintages (backward compatible)
        operation_matrix = {}
        for p in self.m.PROCESS:
            for t in self.m.SYSTEM_TIME:
                total_op = sum(
                    pyo.value(self.m.var_operation[proc, v, time])
                    for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME
                    if proc == p and time == t
                )
                operation_matrix[(t, p)] = total_op
        df = pd.DataFrame.from_dict(operation_matrix, orient="index", columns=["Value"])
        df.index = pd.MultiIndex.from_tuples(df.index, names=["Time", "Process"])
        df = df.reset_index()
        df_pivot = df.pivot(index="Time", columns="Process", values="Value")
        self.df_operation = df_pivot
        return self.df_operation
    else:
        # Return full 3D data with (Process, Vintage) columns
        operation_matrix = {
            (t, p, v): pyo.value(self.m.var_operation[p, v, t])
            for (p, v, t) in self.m.ACTIVE_VINTAGE_TIME
        }
        df = pd.DataFrame.from_dict(operation_matrix, orient="index", columns=["Value"])
        df.index = pd.MultiIndex.from_tuples(df.index, names=["Time", "Process", "Vintage"])
        df = df.reset_index()
        df_pivot = df.pivot(index="Time", columns=["Process", "Vintage"], values="Value")
        return df_pivot

def get_production(self) -> pd.DataFrame:
    """
    Extracts the production data from the model and returns it as a DataFrame.
    The DataFrame will have a MultiIndex with 'Process', 'Product', and
    'Time'. The values are the total production for each process and product
    at each time step.

    With 3D var_operation[p, v, t], production is summed across all active
    vintages at each time step.
    """
    production_tensor = {}
    fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0)

    # Get production overrides data from model (if exists)
    production_overrides = getattr(self.m, "_production_vintage_overrides", {})
    production_overrides_index = getattr(self.m, "_production_overrides_index", frozenset())

    def get_production_value(p, r, tau, vintage):
        """Get production value, checking sparse overrides first."""
        key = (p, r, tau, vintage)
        if key in production_overrides:
            return production_overrides[key]
        return pyo.value(self.m.foreground_production[p, r, tau])

    def has_production_overrides(p, r):
        """Check if any vintage overrides exist for this process/product."""
        return (p, r) in production_overrides_index

    for p in self.m.PROCESS:
        for f in self.m.PRODUCT:
            for t in self.m.SYSTEM_TIME:
                # Sum production across all active vintages at time t
                total_production = 0
                for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME:
                    if proc != p or time != t:
                        continue

                    tau = t - v  # lifecycle stage of this vintage in year t
                    if tau not in self.m.PROCESS_TIME:
                        continue

                    # Annual output per running unit at this lifecycle stage
                    if has_production_overrides(p, f):
                        production_rate = pyo.value(get_production_value(p, f, tau, v))
                    else:
                        production_rate = pyo.value(
                            self.m.foreground_production[p, f, tau]
                        )

                    # Production from this vintage
                    total_production += production_rate * pyo.value(self.m.var_operation[p, v, t])

                production_tensor[(p, f, t)] = total_production * fg_scale

    df = pd.DataFrame.from_dict(
        production_tensor, orient="index", columns=["Value"]
    )
    df.index = pd.MultiIndex.from_tuples(
        df.index, names=["Process", "Product", "Time"]
    )
    df = df.reset_index()
    df_pivot = df.pivot(
        index="Time", columns=["Process", "Product"], values="Value"
    )
    self.df_production = df_pivot
    return self.df_production

def get_demand(self) -> pd.DataFrame:
    """
    Extracts the demand data from the model and returns it as a DataFrame.
    The DataFrame will have a MultiIndex with 'Product' and 'Time'.
    The values are the demand for each Product at each time step.
    """
    fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0)
    demand_matrix = {
        (f, t): self.m.demand[f, t] * fg_scale
        for f in self.m.PRODUCT
        for t in self.m.SYSTEM_TIME
    }
    df = pd.DataFrame.from_dict(demand_matrix, orient="index", columns=["Value"])
    df.index = pd.MultiIndex.from_tuples(
        df.index, names=["Product", "Time"]
    )
    df = df.reset_index()
    df_pivot = df.pivot(index="Time", columns="Product", values="Value")
    self.df_demand = df_pivot
    return self.df_demand

def plot_impacts(self, df_impacts=None, annotated=True):
    """
    Plot a stacked bar chart for impacts by category and process over time.

    Creates a figure with one subplot per impact category, showing process
    contributions as stacked bars. Automatically denormalizes scaled values
    and optionally displays human-readable process names.

    Parameters
    ----------
    df_impacts : DataFrame, optional
        DataFrame with Time as index, Categories and Processes as columns.
        Columns must be a MultiIndex: (Category, Process). If not provided,
        automatically extracted via get_impacts().
    annotated : bool, default=True
        If True, show human-readable names from Brightway database instead
        of process codes.
    """
    if df_impacts is None:
        if getattr(self, "df_impacts", None) is not None:
            df_impacts = self.df_impacts
        else:
            df_impacts = self.get_impacts()

    categories = df_impacts.columns.get_level_values(0).unique()
    n_categories = len(categories)
    ncols = min(self._plot_config["subplot_ncols"], n_categories)
    nrows = math.ceil(n_categories / ncols)

    base_w, base_h = self._plot_config["figsize"]
    fig_w = base_w * ncols
    fig_h = base_h * nrows

    fig, axes = self._create_clean_axes(nrows=nrows, ncols=ncols, figsize=(fig_w, fig_h))

    all_handles = []
    all_labels = []
    for i, category in enumerate(categories):
        ax = axes[i]
        sub_df = df_impacts[category]
        # Filter out processes with all zero values in this category
        sub_df = sub_df.loc[:, (sub_df != 0).any(axis=0)]
        # Get colors BEFORE annotation (using codes)
        colors = self._get_colors_for_dataframe(sub_df)
        # Annotate if requested
        sub_df = self._annotate_dataframe(sub_df, annotated)
        self._apply_bar_styles(sub_df, ax, colors, title=category)
        ax.set_xlabel("")
        ax.set_ylabel("Impact", fontsize=self._plot_config["label_fontsize"])

        # Collect handles/labels for shared legend
        h, l = ax.get_legend_handles_labels()
        for handle, label in zip(h, l):
            if label not in all_labels:
                all_handles.append(handle)
                all_labels.append(label)

    # Hide unused axes
    for j in range(i + 1, len(axes)):
        fig.delaxes(axes[j])

    # Shared legend to the right of figure
    if all_handles:
        all_handles, all_labels = self._reorder_legend_row_first(all_handles, all_labels, 2)
        fig.legend(
            handles=all_handles,
            labels=all_labels,
            loc="upper center",
            bbox_to_anchor=(0.5, -0.02),
            ncol=2,
            fontsize=self._plot_config["legend_fontsize"],
            frameon=False,
        )

    fig.tight_layout()
    plt.show()

def plot_installation(self, df_installation=None, annotated=True):
    """
    Plot a stacked bar chart for installation data.

    Bars show `var_installation`: the number of process units built in each vintage
    year, a LIFETIME quantity plotted at the year of installation. This is not an
    annual capacity and does not line up with production in the same year - use
    `plot_capacity_balance()` to compare production against the annual capacity
    those units make available.

    Parameters
    ----------
    df_installation : DataFrame, optional
        DataFrame with Time as index, Processes as columns
    annotated : bool, default=True
        If True, show human-readable names instead of codes
    """
    if df_installation is None:
        df_installation = self.get_installation()

    # Filter out columns with all zero values
    df_installation = df_installation.loc[:, (df_installation != 0).any(axis=0)]

    # Get colors BEFORE annotation (using codes)
    colors = self._get_colors_for_dataframe(df_installation)

    # Annotate if requested
    df_installation = self._annotate_dataframe(df_installation, annotated)

    fig, axes = self._create_clean_axes()
    ax = axes[0]
    self._apply_bar_styles(
        df_installation, ax, colors, title="Installed Units by Vintage"
    )
    ax.set_ylabel(
        "Installed units (lifetime)", fontsize=self._plot_config["label_fontsize"]
    )


    # Legend at bottom
    h, l = ax.get_legend_handles_labels()
    if h:
        h, l = self._reorder_legend_row_first(h, l, 2)
        fig.legend(
            handles=h,
            labels=l,
            loc="upper center",
            bbox_to_anchor=(0.5, -0.02),
            ncol=2,
            fontsize=self._plot_config["legend_fontsize"],
            frameon=False,
        )

    fig.tight_layout()
    plt.show()

def plot_operation(self, df_operation=None, annotated=True):
    """
    Plot a stacked bar chart for operation data.

    Bars show `var_operation` summed over vintages: how many units run in each year.
    Note the different meaning of the x-axis compared to `plot_installation()`,
    which is indexed by the year of installation.

    Parameters
    ----------
    df_operation : DataFrame, optional
        DataFrame with Time as index, Processes as columns
    annotated : bool, default=True
        If True, show human-readable names instead of codes
    """
    if df_operation is None:
        df_operation = self.get_operation()

    # Filter out columns with all zero values
    df_operation = df_operation.loc[:, (df_operation != 0).any(axis=0)]

    # Get colors BEFORE annotation (using codes)
    colors = self._get_colors_for_dataframe(df_operation)

    # Annotate if requested
    df_operation = self._annotate_dataframe(df_operation, annotated)

    fig, axes = self._create_clean_axes()
    ax = axes[0]
    self._apply_bar_styles(
        df_operation, ax, colors, title="Units Running"
    )
    ax.set_ylabel("Units running", fontsize=self._plot_config["label_fontsize"])


    # Legend at bottom
    h, l = ax.get_legend_handles_labels()
    if h:
        h, l = self._reorder_legend_row_first(h, l, 2)
        fig.legend(
            handles=h,
            labels=l,
            loc="upper center",
            bbox_to_anchor=(0.5, -0.02),
            ncol=2,
            fontsize=self._plot_config["legend_fontsize"],
            frameon=False,
        )

    fig.tight_layout()
    plt.show()

def get_existing_capacity(self) -> pd.DataFrame:
    """
    Extract existing (brownfield) capacity data from the model.

    Returns a DataFrame showing which processes have existing capacity,
    when they were installed, and their operational status at each time step.

    Returns
    -------
    pd.DataFrame
        DataFrame with Time as index and (Process, Type) as MultiIndex columns.
        Type can be 'existing_capacity' (total existing) or 'existing_operating'
        (existing capacity in operation phase at that time).
    """
    existing_cap_dict = getattr(self.m, "_existing_capacity_dict", {})

    if not existing_cap_dict:
        # Return empty DataFrame if no existing capacity
        return pd.DataFrame()

    data = {}
    for t in self.m.SYSTEM_TIME:
        for p in self.m.PROCESS:
            # Calculate existing capacity in operation at time t
            op_start = pyo.value(self.m.process_operation_start[p])
            op_end = pyo.value(self.m.process_operation_end[p])

            existing_operating = 0
            existing_total = 0

            for (proc, inst_year), capacity in existing_cap_dict.items():
                if proc == p:
                    existing_total += capacity
                    tau_existing = t - inst_year
                    if op_start <= tau_existing <= op_end:
                        existing_operating += capacity

            if existing_total > 0:
                data[(t, p, "existing_capacity")] = existing_total
                data[(t, p, "existing_operating")] = existing_operating

    if not data:
        return pd.DataFrame()

    df = pd.DataFrame.from_dict(data, orient="index", columns=["Value"])
    df.index = pd.MultiIndex.from_tuples(df.index, names=["Time", "Process", "Type"])
    df = df.reset_index()
    df_pivot = df.pivot(index="Time", columns=["Process", "Type"], values="Value")

    return df_pivot

def get_production_capacity(self) -> pd.DataFrame:
    """
    Calculate maximum available ANNUAL production capacity for each product at
    each time step.

    Installed units (`get_installation()`) are counted in process units, and one
    unit yields its full lifetime production over the whole operation window. This
    method converts those units into the output they can deliver *in a given year*:
    for every vintage active at time t, its unit count is multiplied by the
    production coefficient at its current lifecycle stage tau = t - v. This is the
    quantity to compare against actual production (`get_production()`), which is
    also per year. Includes both new installations (from var_installation) and
    existing (brownfield) capacity.

    Returns
    -------
    pd.DataFrame
        DataFrame with Time as index and Products as columns.
        Values represent maximum annual production capacity (not actual production).
    """
    capacity_tensor = {}
    fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0)
    existing_cap_dict = getattr(self.m, "_existing_capacity_dict", {})

    # Get production overrides data from model (if exists)
    production_overrides = getattr(self.m, "_production_vintage_overrides", {})
    production_overrides_index = getattr(self.m, "_production_overrides_index", frozenset())

    def get_production_value(p, r, tau, vintage):
        """Get production value, checking sparse overrides first."""
        key = (p, r, tau, vintage)
        if key in production_overrides:
            return production_overrides[key]
        return pyo.value(self.m.foreground_production[p, r, tau])

    def has_production_overrides(p, r):
        """Check if any vintage overrides exist for this process/product."""
        return (p, r) in production_overrides_index

    def units_available(p, v):
        """Units of vintage v that can run: greenfield installs or brownfield stock."""
        if v in self.m.SYSTEM_TIME:
            return pyo.value(self.m.var_installation[p, v])
        return existing_cap_dict.get((p, v), 0)

    for f in self.m.PRODUCT:
        for t in self.m.SYSTEM_TIME:
            # Sum annual capacity over all vintages active at time t
            total_capacity = 0

            for (p, v, time) in self.m.ACTIVE_VINTAGE_TIME:
                if time != t:
                    continue

                tau = t - v  # lifecycle stage of this vintage in year t
                if tau not in self.m.PROCESS_TIME:
                    continue

                if has_production_overrides(p, f):
                    annual_production_per_unit = get_production_value(p, f, tau, v)
                else:
                    annual_production_per_unit = pyo.value(
                        self.m.foreground_production[p, f, tau]
                    )

                total_capacity += annual_production_per_unit * units_available(p, v)

            # Store denormalized capacity
            capacity_tensor[(f, t)] = total_capacity * fg_scale

    # Convert to DataFrame
    df = pd.DataFrame.from_dict(capacity_tensor, orient="index", columns=["Value"])
    df.index = pd.MultiIndex.from_tuples(df.index, names=["Product", "Time"])
    df = df.reset_index()
    df_pivot = df.pivot(index="Time", columns="Product", values="Value")

    return df_pivot