notebooks/september_2025_validation.py

import marimo

__generated_with = "0.17.2"
app = marimo.App(width="medium")


@app.cell
def imports_and_setup():
    """Import libraries and set up paths."""
    import marimo as mo
    import polars as pl
    import altair as alt
    from pathlib import Path
    from datetime import datetime
    import numpy as np

    # Set up absolute paths
    project_root = Path(__file__).parent.parent
    return alt, datetime, mo, pl, project_root


@app.cell
def load_september_2025_data(datetime, pl, project_root):
    """Load September 2025 forecast results and actuals."""

    # Load actuals from HuggingFace dataset (ground truth)
    print('[INFO] Loading actuals from HuggingFace dataset...')
    from datasets import load_dataset
    import os

    dataset = load_dataset('evgueni-p/fbmc-features-24month', split='train', token=os.environ.get('HF_TOKEN'))
    df_actuals_full = pl.from_arrow(dataset.data.table)
    print(f'[INFO] HF dataset loaded: {df_actuals_full.shape}')

    # Load forecast results (full 14-day forecast with 132 borders)
    forecast_path = project_root / 'results' / 'september_2025_forecast_full_14day.parquet'

    if not forecast_path.exists():
        raise FileNotFoundError(f'Forecast file not found: {forecast_path}. Run September 2025 forecast first.')

    df_forecast_full = pl.read_parquet(forecast_path)
    print(f'[INFO] Forecast loaded: {df_forecast_full.shape}')
    print(f'[INFO] Forecast dates: {df_forecast_full["timestamp"].min()} to {df_forecast_full["timestamp"].max()}')

    # Filter actuals to September 2025 period (Aug 18 - Sept 15 for context + forecast period)
    start_date = datetime(2025, 8, 18)  # 2 weeks before forecast
    end_date = datetime(2025, 9, 16)     # Through end of forecast period

    df_actuals_filtered = df_actuals_full.filter(
        (pl.col('timestamp') >= start_date) &
        (pl.col('timestamp') < end_date)
    )

    print(f'[INFO] Actuals filtered: {df_actuals_filtered.shape[0]} hours (Aug 18 - Sept 15, 2025)')

    # Forecast period for evaluation
    forecast_start = datetime(2025, 9, 2)
    return df_actuals_filtered, df_forecast_full


@app.cell
def prepare_unified_dataframe(
    datetime,
    df_actuals_filtered,
    df_forecast_full,
    pl,
):
    """Prepare unified dataframe with forecast and actual pairs for ALL FBMC borders."""

    # Extract ALL border names from forecast columns (132 directional borders)
    # Includes both physical interconnectors and virtual trading paths
    forecast_cols_list = [col for col in df_forecast_full.columns if col.endswith('_median')]
    border_names_list = [col.replace('_median', '') for col in forecast_cols_list]

    print(f'[INFO] Processing {len(border_names_list)} FBMC borders (all directional trading paths)...')
    print(f'[INFO] Sample borders: {sorted(border_names_list)[:10]}...')

    # Start with timestamp from actuals
    df_unified_data = df_actuals_filtered.select('timestamp')

    # Add actual and forecast for each border
    for border in border_names_list:
        actual_col_source = f'target_border_{border}'
        forecast_col_source = f'{border}_median'

        # Add actuals
        if actual_col_source in df_actuals_filtered.columns:
            df_unified_data = df_unified_data.with_columns(
                df_actuals_filtered[actual_col_source].alias(f'actual_{border}')
            )
        else:
            print(f'[WARNING] Actual column missing: {actual_col_source}')
            df_unified_data = df_unified_data.with_columns(pl.lit(None).alias(f'actual_{border}'))

        # Add forecasts (join on timestamp)
        if forecast_col_source in df_forecast_full.columns:
            df_forecast_subset = df_forecast_full.select(['timestamp', forecast_col_source])
            df_unified_data = df_unified_data.join(
                df_forecast_subset,
                on='timestamp',
                how='left'
            ).rename({forecast_col_source: f'forecast_{border}'})
        else:
            print(f'[WARNING] Forecast column missing: {forecast_col_source}')
            df_unified_data = df_unified_data.with_columns(pl.lit(None).alias(f'forecast_{border}'))

    print(f'[INFO] Unified data prepared: {df_unified_data.shape}')

    # Validate no data leakage - check that forecasts don't perfectly match actuals
    sample_border = border_names_list[0]
    forecast_col_check = f'forecast_{sample_border}'
    actual_col_check = f'actual_{sample_border}'

    if forecast_col_check in df_unified_data.columns and actual_col_check in df_unified_data.columns:
        _forecast_start_check = datetime(2025, 9, 2)
        _df_forecast_check = df_unified_data.filter(pl.col('timestamp') >= _forecast_start_check)

        if len(_df_forecast_check) > 0:
            mae_check = (_df_forecast_check[forecast_col_check] - _df_forecast_check[actual_col_check]).abs().mean()
            if mae_check == 0:
                raise ValueError(f'DATA LEAKAGE DETECTED: Forecasts perfectly match actuals (MAE=0) for {sample_border}!')

    print('[INFO] Data leakage check passed - forecasts differ from actuals')
    return border_names_list, df_unified_data


@app.cell
def create_border_selector(border_names_list, mo):
    """Create interactive border selection dropdown."""

    border_selector_widget = mo.ui.dropdown(
        options={border: border for border in sorted(border_names_list)},
        value='CZ_PL',  # Default to Polish border to showcase fix
        label='Select Border:'
    )
    return (border_selector_widget,)


@app.cell
def display_border_selector(border_selector_widget, mo):
    """Display the border selector UI."""
    mo.md(f"""
    ## Forecast Validation: September 2025 (All FBMC Borders)

    **Select a border to view:**
    {border_selector_widget}

    Chart shows:
    - **2 weeks historical** (Aug 18 - Sept 1, 2025): Actual flows only
    - **2 weeks forecast** (Sept 2-15, 2025): Forecast vs Actual comparison
    - **Context**: 336 hours forecast period (14 days)
    - **Borders shown**: All 132 FBMC directional borders (66 country pairs x 2 directions)
    - **Note**: Includes both physical interconnectors and virtual trading paths
    """)
    return


@app.cell
def filter_data_for_selected_border(
    border_selector_widget,
    datetime,
    df_unified_data,
    pl,
):
    """Filter data for the selected border."""

    selected_border_name = border_selector_widget.value

    # Extract columns for selected border
    actual_col_name = f'actual_{selected_border_name}'
    forecast_col_name = f'forecast_{selected_border_name}'

    # Check if columns exist
    if actual_col_name not in df_unified_data.columns:
        df_selected_border = None
        print(f'[ERROR] Actual column {actual_col_name} not found')
    else:
        df_selected_border = df_unified_data.select([
            'timestamp',
            pl.col(actual_col_name).alias('actual'),
            pl.col(forecast_col_name).alias('forecast') if forecast_col_name in df_unified_data.columns else pl.lit(None).alias('forecast')
        ])

        # Add period marker (historical vs forecast)
        forecast_start_time = datetime(2025, 9, 2)
        df_selected_border = df_selected_border.with_columns(
            pl.when(pl.col('timestamp') >= forecast_start_time)
            .then(pl.lit('Forecast Period'))
            .otherwise(pl.lit('Historical'))
            .alias('period')
        )
    return df_selected_border, forecast_start_time, selected_border_name


@app.cell
def create_time_series_chart(
    alt,
    df_selected_border,
    forecast_start_time,
    selected_border_name,
):
    """Create Altair time series visualization."""

    if df_selected_border is None:
        chart_time_series = alt.Chart().mark_text(text='No data available', size=20)
    else:
        # Convert to pandas for Altair (CLAUDE.md Rule #37)
        df_plot = df_selected_border.to_pandas()

        # Create base chart
        base = alt.Chart(df_plot).encode(
            x=alt.X('timestamp:T', title='Date', axis=alt.Axis(format='%b %d'))
        )

        # Actual line (blue, solid)
        line_actual = base.mark_line(color='blue', strokeWidth=2).encode(
            y=alt.Y('actual:Q', title='Flow (MW)', scale=alt.Scale(zero=False)),
            tooltip=[
                alt.Tooltip('timestamp:T', title='Time', format='%Y-%m-%d %H:%M'),
                alt.Tooltip('actual:Q', title='Actual (MW)', format='.0f')
            ]
        )

        # Forecast line (orange, dashed) - only for forecast period
        df_plot_forecast = df_plot[df_plot['period'] == 'Forecast Period']

        if len(df_plot_forecast) > 0 and df_plot_forecast['forecast'].notna().any():
            line_forecast = alt.Chart(df_plot_forecast).mark_line(
                color='orange',
                strokeWidth=2,
                strokeDash=[5, 5]
            ).encode(
                x=alt.X('timestamp:T'),
                y=alt.Y('forecast:Q'),
                tooltip=[
                    alt.Tooltip('timestamp:T', title='Time', format='%Y-%m-%d %H:%M'),
                    alt.Tooltip('forecast:Q', title='Forecast (MW)', format='.0f'),
                    alt.Tooltip('actual:Q', title='Actual (MW)', format='.0f')
                ]
            )
        else:
            line_forecast = alt.Chart().mark_point()  # Empty chart

        # Vertical line at forecast start
        rule_forecast_start = alt.Chart(
            alt.Data(values=[{'x': forecast_start_time}])
        ).mark_rule(color='red', strokeDash=[3, 3], strokeWidth=1).encode(
            x='x:T'
        )

        # Combine layers
        chart_time_series = (line_actual + line_forecast + rule_forecast_start).properties(
            width=800,
            height=400,
            title=f'Border: {selected_border_name} | Hourly Flows (Aug 18 - Sept 15, 2025)'
        ).configure_axis(
            labelFontSize=12,
            titleFontSize=14
        ).configure_title(
            fontSize=16
        )
    return (chart_time_series,)


@app.cell
def calculate_summary_statistics(
    df_selected_border,
    forecast_start_time,
    pl,
    selected_border_name,
):
    """Calculate comprehensive evaluation metrics for the selected border."""

    if df_selected_border is None:
        stats_summary_text = 'No data available'
    else:
        # Filter to forecast period only
        df_forecast_period = df_selected_border.filter(
            pl.col('timestamp') >= forecast_start_time
        )

        if len(df_forecast_period) == 0 or df_forecast_period['forecast'].is_null().all():
            stats_summary_text = 'No forecast data available for this period'
        else:
            # Calculate accuracy metrics
            forecast_vals = df_forecast_period['forecast'].drop_nulls()
            actual_vals = df_forecast_period['actual'].drop_nulls()

            # Align forecast and actual (remove any nulls)
            df_eval = df_forecast_period.filter(
                pl.col('forecast').is_not_null() & pl.col('actual').is_not_null()
            )

            if len(df_eval) == 0:
                stats_summary_text = 'No overlapping forecast and actual data'
            else:
                # Error metrics
                errors = (df_eval['forecast'] - df_eval['actual'])
                abs_errors = errors.abs()

                mae_value = abs_errors.mean()
                rmse_value = (errors.pow(2).mean() ** 0.5)
                mape_value = (abs_errors / df_eval['actual'].abs()).mean() * 100

                # Bias metrics
                mean_error = errors.mean()

                # Forecast quality metrics
                unique_count = forecast_vals.n_unique()
                std_forecast = forecast_vals.std()
                std_actual = actual_vals.std()

                # Range metrics
                forecast_range = forecast_vals.max() - forecast_vals.min()
                actual_range = actual_vals.max() - actual_vals.min()

                stats_summary_text = f"""
    ### Forecast Quality Metrics

    **Border**: {selected_border_name}
    **Period**: September 2-15, 2025 (336 hours)
    **Evaluation Points**: {len(df_eval)} hours

    #### Accuracy Metrics
    - **MAE** (Mean Absolute Error): {mae_value:.0f} MW
    - **RMSE** (Root Mean Squared Error): {rmse_value:.0f} MW
    - **MAPE** (Mean Absolute Percentage Error): {mape_value:.1f}%
    - **Bias** (Mean Error): {mean_error:+.0f} MW

    #### Forecast Variation
    - **Unique Values**: {unique_count} / {len(df_eval)} ({unique_count/len(df_eval)*100:.0f}%)
    - **Forecast StdDev**: {std_forecast:.0f} MW
    - **Actual StdDev**: {std_actual:.0f} MW
    - **Forecast Range**: {forecast_range:.0f} MW
    - **Actual Range**: {actual_range:.0f} MW

    #### Interpretation
    - **MAE < 150 MW**: ✓ Excellent (zero-shot baseline target)
    - **MAE 150-300 MW**: Good
    - **MAE > 300 MW**: Needs improvement
    - **Variation**: {unique_count} unique values indicates {'VALID time-varying forecast' if unique_count > 50 else 'LOW VARIATION - may be constant'}
    - **Bias**: {'Overforecasting' if mean_error > 50 else 'Underforecasting' if mean_error < -50 else 'Balanced'}
    """
    return (stats_summary_text,)


@app.cell
def display_chart_and_stats(chart_time_series, mo, stats_summary_text):
    """Display the chart and statistics."""
    mo.vstack([
        chart_time_series,
        mo.md(stats_summary_text)
    ])
    return


if __name__ == "__main__":
    app.run()