Blame for code/recipe.py

"""The recipe that is published (docs/forecast-per-region.md, sections 1, 5 and 7), as pure functions on frames of (zone, hour, horizon h) rows.

    zones other than EE / LV / LT, fuel inputs usable   mean( LEAR-365 with fuel , LightGBM with fuel )                 `ens`
    the same zones, fuel NOT usable                      mean( LEAR-365 without fuel , LightGBM without fuel )         `ens-nofuel`
    EE, LV, LT (the Baltics: fuel does not help there)   the current LightGBM without fuel                             `lgb`

LEAR-365 = the research's `lear_r365_g` / `lear_r365`: per zone and horizon a closed-form ridge on the RESIDUAL of the naive 'level +
shape' forecast, trained on the last 365 days of settled targets (window ending the day before the last known day, exactly what was
settled when the first forecast of a block was issued), features region_models.design (110 columns: last day's shape, last week's hour,
level trend, calendar, wind/solar/temperature and their hour interactions, and - with fuel - the gas-plant cost features) plus four
price-relative lags. Refitted every day (the research refit weekly; the daily refit was 1-2 points better and costs seconds).
LightGBM = backtest.py's model: pooled over the zones, `zone` categorical, target = price - level_7d, 400 trees.

Everything here is a function of the rows only, so training rows and live rows are built and augmented by the SAME code (the leakage
rule of features.py and fuel_features.py holds for both) and the whole recipe is unit-tested without LightGBM or the network.
"""
import os
import sys

import numpy as np
import pandas as pd

DIR = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, DIR)
import family_models as FM        # noqa: E402
import fuel_features as FUEL      # noqa: E402
import region_models as RM        # noqa: E402

BALTICS = ("ee", "lv", "lt")
H_MIN, H_MAX = 2, 7
LEAR_WINDOW_DAYS = 365
LEAR_MIN_ROWS = 24 * 40           # fewer training hours than this: no linear model (the zone-horizon is not published)
EXTRA_LAGS = ["d_lag_last", "d_prev", "d_lag7", "d_lag14"]
# Hourly price lags of a LIVE row can be missing where the training rows (built in hindsight) never are: the last known day of Portugal and Ireland
# has no hour 23 (it belongs to the market day that is not published yet), a zone's price feed has a hole of an hour or a day now and then. The
# LightGBM models take a missing lag as it is; the ridge cannot, so such a feature is replaced by its training-window mean (which makes its
# contribution zero) - for at most MAX_IMPUTED of these six per row, and never for weather, calendar, level or fuel features.
IMPUTABLE = ("shape_persist", "lag7_dev", "d_lag_last", "d_prev", "d_lag7", "d_lag14")
MAX_IMPUTED = 3
RECIPES = {"ens": "mean of a 365-day rolling ridge and LightGBM, with gas and carbon inputs",
           "ens-nofuel": "mean of a 365-day rolling ridge and LightGBM, without fuel inputs",
           "lgb": "LightGBM"}
MODEL_NAME = "eu-ens-v2"


def recipe_for(zone, fuel_ok):
    """'lgb' for the Baltics, otherwise 'ens' when the fuel inputs are usable and 'ens-nofuel' when not."""
    if zone in BALTICS:
        return "lgb"
    return "ens" if fuel_ok else "ens-nofuel"


# ------------------------------------------------------------------ rows: fuel and the columns the linear model needs
def add_fuel_columns(df, table):
    """Rows with `date` (a normalised timestamp) and `h` -> the same rows with the leak-free fuel features and `mc_gap`
    (fuel_dataset.add_fuel, applied to a frame that already has `date`)."""
    out = FUEL.add_to_frame(df, table)
    out["mc_gap"] = out["ccgt_mc"] - out["level_7d"]
    return out


def augment(df, prices, tables=None):
    """Add what the linear model needs on top of features.make_features (+ fuel): ewm5 and the naive 'level + shape' baseline `b_lvlshape`
    (the point the ridge corrects), and the price-relative lags. `df` has zone (str), t, date, h; `prices` is {zone: hourly Series on the
    zone's clock}. Every added column uses only prices up to the row's last known day (target - (h - 1)); the lags k >= h - 1."""
    out = df.copy()
    out["zone"] = out["zone"].astype(str)
    tables = tables or RM.day_tables(prices)
    out = RM.attach_last_known(out, tables, ["ewm5"])
    out["shape_s"] = RM.shape(out)
    out["b_lvlshape"] = out["shape_s"] + out["level_7d"]
    out["p_prev"] = np.nan
    for (z, h), idx in out.groupby(["zone", "h"]).groups.items():
        out.loc[idx, "p_prev"] = prices[z].reindex(out.loc[idx, "t"] - pd.Timedelta(days=int(h))).to_numpy()
    out["d_lag_last"] = out["lag_last"] - out["level_7d"]
    out["d_prev"] = out["p_prev"] - out["level_7d"]
    out["d_lag7"] = out["lag7"] - out["level_7d"]
    out["d_lag14"] = out["lag14"] - out["level_7d"]
    return out


def lear_design(df, fuel):
    """The ridge's feature matrix (float32 frame): region_models.design + the four price-relative lags."""
    return RM.design(df, fuel).join(df[EXTRA_LAGS].astype(np.float32))


# ------------------------------------------------------------------ LEAR-365
def lear_predict(history, live, fuel, window=LEAR_WINDOW_DAYS, min_rows=LEAR_MIN_ROWS, impute=True):
    """Predict the `live` rows with one ridge per (zone, h).

    history  settled training rows (zone, h, date, y, b_lvlshape + features), already augmented
    live     rows to predict (same columns, y absent), already augmented
    For a (zone, h) group of live rows whose first target day is T, the training rows are those of the same zone and h whose target day
    lies in (T - h - window, T - h]: T - h is the day before the last known day (T - h + 1), the newest target already settled when the
    forecast is issued. Returns a Series aligned to live.index; NaN where there are too few training hours or a feature is missing -
    never a patched value. The one exception is `impute` (see IMPUTABLE): a missing hourly price lag takes the training mean. The number of
    imputed values is in the result's `.attrs["imputed"]`."""
    out = pd.Series(np.nan, index=live.index)
    out.attrs["imputed"] = 0
    if live.empty:
        return out
    Xl = lear_design(live, fuel)
    imp = [i for i, c in enumerate(Xl.columns) if c in IMPUTABLE] if impute else []
    hist = history[history["h"].isin(live["h"].unique())]
    for (z, h), g in live.groupby(["zone", "h"]):
        first = g["date"].min()
        end = first - pd.Timedelta(days=int(h))
        tr = hist[(hist["zone"] == z) & (hist["h"] == h) & (hist["date"] <= end) & (hist["date"] > end - pd.Timedelta(days=window))]
        if len(tr) < min_rows:
            continue
        Xt = lear_design(tr, fuel).to_numpy(np.float64)
        yt = (tr["y"] - tr["b_lvlshape"]).to_numpy(np.float64)
        ok = np.isfinite(Xt).all(axis=1) & np.isfinite(yt)
        if ok.sum() < min_rows:
            continue
        model = FM.ridge_fit(Xt[ok], yt[ok], np.ones(int(ok.sum())))
        Xg = Xl.loc[g.index].to_numpy(np.float64)
        if imp:
            nan = ~np.isfinite(Xg)
            fill = nan[:, imp]
            other = nan.copy()
            other[:, imp] = False
            usable = ~other.any(axis=1) & (fill.sum(axis=1) <= MAX_IMPUTED)
            for k, j in enumerate(imp):
                m = fill[:, k] & usable
                Xg[m, j] = model[0][j]
                out.attrs["imputed"] += int(m.sum())
        good = np.isfinite(Xg).all(axis=1) & g["b_lvlshape"].notna().to_numpy()
        if good.any():
            out.loc[g.index[good]] = g.loc[g.index[good], "b_lvlshape"].to_numpy() + FM.linear_predict(model, Xg[good])
    return out


# ------------------------------------------------------------------ the ensemble
MAX_FALLBACK_HOURS = 6


def point_forecast(zone, parts, fuel_ok, day_ids=None, max_fallback_hours=MAX_FALLBACK_HOURS):
    """The published central estimate of one zone from the model columns `parts` (dict name -> array/Series aligned to the zone's rows):
    lear_fuel, lgb_fuel, lear_nofuel, lgb_base. NaN wherever a needed component is NaN (the day is then dropped by the caller)...

    ... with ONE exception, for continuity: when only the ridge is missing in an hour (a gap in the price feed that imputation could not bridge) and
    `day_ids` (one id per row: the target day) is given, the LightGBM alone stands in for that hour, if at most `max_fallback_hours` hours of the same day
    need it. The API shows forecast days only without a hole after the settled days, so one lost day would hide every day after it. Each model alone was
    about as much better than level+shape as the mean at D+2 (docs/forecast-per-region.md, section 1); a day that needs it for more than a quarter of
    its hours (a zone with no ridge at all) is dropped, not published on a model the gate did not measure."""
    r = recipe_for(zone, fuel_ok)
    if r == "lgb":
        return np.asarray(parts["lgb_base"], dtype=float)
    a, b = ("lear_fuel", "lgb_fuel") if r == "ens" else ("lear_nofuel", "lgb_base")
    a, b = np.asarray(parts[a], dtype=float), np.asarray(parts[b], dtype=float)
    out = (a + b) / 2
    if day_ids is not None:
        only_tree = ~np.isfinite(a) & np.isfinite(b)
        if only_tree.any():
            n = pd.Series(only_tree).groupby(np.asarray(day_ids)).transform("sum").to_numpy()
            out = np.where(only_tree & (n <= max_fallback_hours), b, out)
    return out