fix: filter invalid solar MVLR observations

Author: Molier

openenergyid/mvlr/__init__.py

@@ -6,14 +6,18 @@
     IndependentVariableResult,
     MultiVariableRegressionInput,
     MultiVariableRegressionResult,
+    OutlierFilteringDiagnostics,
     ValidationParameters,
 )
+from .source_data_filtering import clean_regression_frame
 
 __all__ = [
     "find_best_mvlr",
+    "clean_regression_frame",
     "IndependentVariableInput",
     "MultiVariableRegressionInput",
     "MultiVariableRegressionResult",
+    "OutlierFilteringDiagnostics",
     "ValidationParameters",
     "IndependentVariableResult",
 ]

openenergyid/mvlr/main.py

@@ -3,15 +3,19 @@
 from .helpers import resample_input_data
 from .models import MultiVariableRegressionInput, MultiVariableRegressionResult
 from .mvlr import MultiVariableLinearRegression
+from .source_data_filtering import clean_regression_frame
 
 
 def find_best_mvlr(
     data: MultiVariableRegressionInput,
 ) -> MultiVariableRegressionResult:
     """Cycle through multiple granularities and return the best model."""
     best_rsquared = 0
+    best_filtering = None
     for granularity in data.granularities:
         frame = data.data_frame()
+        frame, filtering = clean_regression_frame(frame, data.dependent_variable)
+        best_filtering = filtering
         frame = resample_input_data(data=frame, granularity=granularity)
         mvlr = MultiVariableLinearRegression(
             data=frame,
@@ -27,8 +31,17 @@ def find_best_mvlr(
             max_f_pvalue=data.validation_parameters.f_pvalue,
             max_pvalues=data.validation_parameters.pvalues,
         ):
-            return MultiVariableRegressionResult.from_mvlr(mvlr)
+            result = MultiVariableRegressionResult.from_mvlr(mvlr)
+            result.outlier_filtering = filtering
+            return result
         best_rsquared = max(best_rsquared, mvlr.fit.rsquared_adj)
-    raise ValueError(
-        f"No valid model found. Best R²: {best_rsquared:.3f} (need ≥{data.validation_parameters.rsquared})"
+    detail = (
+        f"No valid model found. Best R²: {best_rsquared:.3f} "
+        f"(need ≥{data.validation_parameters.rsquared})"
     )
+    if best_filtering and best_filtering.applied:
+        detail += (
+            f"; outlier filtering removed {best_filtering.removed_observation_count}/"
+            f"{best_filtering.original_observation_count} observations"
+        )
+    raise ValueError(detail)

openenergyid/mvlr/models.py

@@ -200,6 +200,23 @@ def from_fit(cls, fit: fm.ols, name: str) -> "IndependentVariableResult":
         )
 
 
+class OutlierFilteringDiagnostics(BaseModel):
+    """Diagnostics about source observations removed before regression fitting."""
+
+    enabled: bool = True
+    original_observation_count: int = Field(alias="originalObservationCount")
+    retained_observation_count: int = Field(alias="retainedObservationCount")
+    removed_observation_count: int = Field(alias="removedObservationCount")
+    removed_non_finite_count: int = Field(alias="removedNonFiniteCount", default=0)
+    removed_negative_count: int = Field(alias="removedNegativeCount", default=0)
+    removed_zero_with_solar_count: int = Field(alias="removedZeroWithSolarCount", default=0)
+    removed_ratio_outlier_count: int = Field(alias="removedRatioOutlierCount", default=0)
+    applied: bool
+    reason: str | None = None
+
+    model_config = ConfigDict(populate_by_name=True)
+
+
 class MultiVariableRegressionResult(BaseModel):
     """Result of a multivariable regression model."""
 
@@ -212,6 +229,10 @@ class MultiVariableRegressionResult(BaseModel):
     intercept: IndependentVariableResult
     granularity: Granularity
     frame: TimeDataFrame
+    outlier_filtering: OutlierFilteringDiagnostics | None = Field(
+        default=None,
+        alias="outlierFiltering",
+    )
 
     model_config = ConfigDict(populate_by_name=True)
 

openenergyid/mvlr/source_data_filtering.py

@@ -0,0 +1,129 @@
+"""Source-data filtering helpers for multi-variable regression."""
+
+import numpy as np
+import pandas as pd
+
+from .models import OutlierFilteringDiagnostics
+
+MINIMUM_RETAINED_FRACTION = 0.50
+MINIMUM_RETAINED_ROWS = 30
+SOLAR_REFERENCE_NAMES = ("solarPowerGeneration", "solarRadiation")
+
+
+def _solar_reference_column(frame: pd.DataFrame) -> str | None:
+    for name in SOLAR_REFERENCE_NAMES:
+        if name in frame.columns:
+            return name
+
+    for column in frame.columns:
+        lower = column.lower()
+        if "solar" in lower and ("generation" in lower or "radiation" in lower):
+            return column
+
+    return None
+
+
+def _is_solar_production_model(dependent_variable: str, frame: pd.DataFrame) -> bool:
+    dependent = dependent_variable.lower()
+    if "solarphotovoltaic" in dependent:
+        return True
+    if "solar" in dependent and "production" in dependent:
+        return True
+    return "production" in dependent and _solar_reference_column(frame) is not None
+
+
+def _positive_reference_threshold(series: pd.Series) -> float:
+    positive = series[series > 0]
+    if positive.empty:
+        return 0.0
+    return max(float(positive.median()) * 0.10, 0.05)
+
+
+def _robust_ratio_outlier_mask(ratio: pd.Series) -> pd.Series:
+    if len(ratio) < MINIMUM_RETAINED_ROWS:
+        return pd.Series(False, index=ratio.index)
+
+    median = float(ratio.median())
+    mad = float((ratio - median).abs().median())
+    if not np.isfinite(mad) or mad <= 0:
+        q1 = float(ratio.quantile(0.25))
+        q3 = float(ratio.quantile(0.75))
+        iqr = q3 - q1
+        if not np.isfinite(iqr) or iqr <= 0:
+            return pd.Series(False, index=ratio.index)
+        return (ratio < q1 - 3.0 * iqr) | (ratio > q3 + 3.0 * iqr)
+
+    robust_z = 0.6745 * (ratio - median).abs() / mad
+    return robust_z > 4.5
+
+
+def clean_regression_frame(
+    frame: pd.DataFrame,
+    dependent_variable: str,
+) -> tuple[pd.DataFrame, OutlierFilteringDiagnostics]:
+    """Remove obvious bad source observations before fitting a regression model."""
+
+    original_count = len(frame)
+    diagnostics = OutlierFilteringDiagnostics(
+        originalObservationCount=original_count,
+        retainedObservationCount=original_count,
+        removedObservationCount=0,
+        applied=False,
+    )
+
+    if original_count == 0 or dependent_variable not in frame.columns:
+        diagnostics.reason = "empty frame or missing dependent variable"
+        return frame, diagnostics
+
+    numeric_frame = frame.apply(pd.to_numeric, errors="coerce")
+    keep = pd.Series(True, index=numeric_frame.index)
+
+    finite_mask = np.isfinite(numeric_frame).all(axis=1)
+    diagnostics.removed_non_finite_count = int((keep & ~finite_mask).sum())
+    keep &= finite_mask
+
+    if not _is_solar_production_model(dependent_variable, numeric_frame):
+        cleaned = numeric_frame.loc[keep].copy()
+        diagnostics.retained_observation_count = len(cleaned)
+        diagnostics.removed_observation_count = original_count - len(cleaned)
+        diagnostics.applied = diagnostics.removed_observation_count > 0
+        diagnostics.reason = "generic non-finite filtering only"
+        return cleaned, diagnostics
+
+    y = numeric_frame[dependent_variable]
+    negative_mask = y < 0
+    diagnostics.removed_negative_count = int((keep & negative_mask).sum())
+    keep &= ~negative_mask
+
+    solar_column = _solar_reference_column(numeric_frame)
+    if solar_column is not None:
+        solar_reference = numeric_frame[solar_column]
+        solar_threshold = _positive_reference_threshold(solar_reference[keep])
+
+        zero_with_solar_mask = (y <= 0) & (solar_reference > solar_threshold)
+        diagnostics.removed_zero_with_solar_count = int((keep & zero_with_solar_mask).sum())
+        keep &= ~zero_with_solar_mask
+
+        ratio_candidates = keep & (y > 0) & (solar_reference > solar_threshold)
+        ratios = y[ratio_candidates] / solar_reference[ratio_candidates]
+        ratio_outliers = _robust_ratio_outlier_mask(ratios)
+        diagnostics.removed_ratio_outlier_count = int(ratio_outliers.sum())
+        keep.loc[ratio_outliers[ratio_outliers].index] = False
+
+    cleaned = numeric_frame.loc[keep].copy()
+    retained_count = len(cleaned)
+    removed_count = original_count - retained_count
+
+    if retained_count < MINIMUM_RETAINED_ROWS:
+        diagnostics.reason = "too few observations retained after filtering"
+        return numeric_frame, diagnostics
+
+    if retained_count / original_count < MINIMUM_RETAINED_FRACTION:
+        diagnostics.reason = "too much source data would be removed"
+        return numeric_frame, diagnostics
+
+    diagnostics.retained_observation_count = retained_count
+    diagnostics.removed_observation_count = removed_count
+    diagnostics.applied = removed_count > 0
+    diagnostics.reason = "solar production source-data filtering"
+    return cleaned, diagnostics

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "openenergyid"
-version = "0.1.40"
+version = "0.1.41"
 description = "Open Source Python library for energy analytics and simulations"
 authors = [
     { name = "Jan Pecinovsky", email = "jan@energieid.be" },

tests/mvlr/test_source_data_filtering.py

@@ -0,0 +1,104 @@
+"""Tests for MVLR source-data filtering."""
+
+import math
+
+import pandas as pd
+
+from openenergyid.models import TimeDataFrame
+from openenergyid.mvlr import (
+    MultiVariableRegressionInput,
+    clean_regression_frame,
+    find_best_mvlr,
+)
+
+DEPENDENT = "energyProduction/solarPhotovoltaic"
+SOLAR_REFERENCE = "solarPowerGeneration"
+
+
+def _solar_regression_input(
+    *,
+    zero_slice: slice = slice(20, 30),
+    spikes: dict[int, float] | None = None,
+) -> MultiVariableRegressionInput:
+    index = pd.date_range("2025-04-01", periods=90, freq="D", tz="Europe/Brussels")
+    solar_reference = pd.Series(
+        [1.0 + (i % 35) / 8.0 for i in range(len(index))],
+        index=index,
+        dtype=float,
+    )
+    production = 4.35 * solar_reference + 0.15
+
+    if zero_slice:
+        production.iloc[zero_slice] = 0.0
+    for idx, value in (spikes or {45: 70.0, 60: 75.0, 75: 80.0}).items():
+        production.iloc[idx] = value
+
+    frame = pd.DataFrame(
+        {
+            DEPENDENT: production,
+            SOLAR_REFERENCE: solar_reference,
+        },
+        index=index,
+    )
+
+    return MultiVariableRegressionInput.model_validate(
+        {
+            "timeZone": "Europe/Brussels",
+            "independentVariables": [
+                {
+                    "name": SOLAR_REFERENCE,
+                    "allowNegativeCoefficient": False,
+                },
+            ],
+            "dependentVariable": DEPENDENT,
+            "frame": TimeDataFrame.from_pandas(frame).model_dump(),
+            "granularities": ["P1D"],
+            "allowNegativePredictions": False,
+            "validationParameters": {
+                "rsquared": 0.95,
+                "f_pvalue": 0.05,
+                "pvalues": 0.05,
+            },
+        },
+    )
+
+
+def test_clean_regression_frame_removes_solar_source_outliers() -> None:
+    """Solar production cleaning should drop zero-line and ratio outliers."""
+    data = _solar_regression_input()
+    frame = data.data_frame()
+
+    cleaned, diagnostics = clean_regression_frame(frame, DEPENDENT)
+
+    assert diagnostics.applied
+    assert diagnostics.original_observation_count == 90
+    assert diagnostics.removed_zero_with_solar_count == 10
+    assert diagnostics.removed_ratio_outlier_count == 3
+    assert diagnostics.removed_observation_count == 13
+    assert len(cleaned) == 77
+    assert (cleaned[DEPENDENT] > 0).all()
+
+
+def test_find_best_mvlr_returns_filtering_diagnostics() -> None:
+    """A model should fit after bad source observations are excluded."""
+    data = _solar_regression_input()
+
+    result = find_best_mvlr(data)
+
+    assert result.r2 > 0.99
+    assert result.outlier_filtering is not None
+    assert result.outlier_filtering.applied
+    assert result.outlier_filtering.removed_observation_count == 13
+    assert math.isclose(result.independent_variables[0].coef, 4.35, rel_tol=0.02)
+
+
+def test_clean_regression_frame_keeps_original_data_when_filtering_too_much() -> None:
+    """Filtering should not apply when too little source data would remain."""
+    data = _solar_regression_input(zero_slice=slice(0, 55), spikes={})
+    frame = data.data_frame()
+
+    cleaned, diagnostics = clean_regression_frame(frame, DEPENDENT)
+
+    assert not diagnostics.applied
+    assert diagnostics.reason == "too much source data would be removed"
+    assert len(cleaned) == 90