Previsione Energie Rinnovabili con ML: Python LSTM per Solare e Eolico

import numpy as np
import pandas as pd
from pvlib import solarposition, irradiance, atmosphere
from sklearn.preprocessing import MinMaxScaler

def create_solar_features(df: pd.DataFrame,
                           lat: float,
                           lon: float,
                           altitude: float = 0,
                           panel_tilt: float = 30,
                           panel_azimuth: float = 180) -> pd.DataFrame:
    """
    Costruisce feature complete per forecasting solare.

    Args:
        df: DataFrame con colonne [ghi, dni, dhi, temp_air, wind_speed, cloud_cover]
            e indice DatetimeTzIndex
        lat, lon: coordinate impianto
        altitude: quota s.l.m. [m]
        panel_tilt: inclinazione pannelli [gradi]
        panel_azimuth: azimut pannelli [gradi, 180=sud]

    Returns:
        DataFrame con tutte le feature engineered
    """
    feat = df.copy()
    times = feat.index

    # --- Geometria solare ---
    solpos = solarposition.get_solarposition(times, lat, lon, altitude)
    feat['solar_zenith'] = solpos['apparent_zenith']
    feat['solar_azimuth'] = solpos['azimuth']
    feat['solar_elevation'] = solpos['apparent_elevation']
    feat['solar_elevation_rad'] = np.radians(feat['solar_elevation'])

    # Angolo di incidenza pannello (AOI)
    aoi = irradiance.aoi(panel_tilt, panel_azimuth,
                         solpos['apparent_zenith'],
                         solpos['azimuth'])
    feat['aoi'] = aoi
    feat['cos_aoi'] = np.cos(np.radians(aoi)).clip(0, 1)

    # --- Irradianza sul piano inclinato (POA) ---
    airmass = atmosphere.get_relative_airmass(solpos['apparent_zenith'])
    poa = irradiance.get_total_irradiance(
        panel_tilt, panel_azimuth,
        solpos['apparent_zenith'], solpos['azimuth'],
        feat['dni'], feat['ghi'], feat['dhi'],
        airmass=airmass
    )
    feat['poa_global'] = poa['poa_global']
    feat['poa_direct'] = poa['poa_direct']
    feat['poa_diffuse'] = poa['poa_diffuse']

    # --- Temperatura cella (modello NOCT) ---
    NOCT = 45  # grado Celsius tipico modulo standard
    feat['temp_cell'] = (feat['temp_air'] +
                         (NOCT - 20) * feat['poa_global'] / 800)

    # Derating per temperatura (coefficiente -0.0045 /degC)
    feat['temp_derating'] = 1 + (-0.0045) * (feat['temp_cell'] - 25)

    # --- Clearness Index (kt) ---
    # Extra-terrestrial radiation per calcolo clearness
    et_rad = irradiance.get_extra_radiation(times)
    feat['clearness_index'] = (feat['ghi'] /
                                (et_rad * np.cos(np.radians(
                                    solpos['apparent_zenith']))
                                ).clip(lower=1e-3)).clip(0, 1.2)

    # --- Cloud cover encoding ---
    feat['cloud_cover_norm'] = feat['cloud_cover'] / 100.0
    feat['clear_sky_fraction'] = 1.0 - feat['cloud_cover_norm']

    # --- Encoding ciclico temporale ---
    feat['hour_sin'] = np.sin(2 * np.pi * times.hour / 24)
    feat['hour_cos'] = np.cos(2 * np.pi * times.hour / 24)
    feat['doy_sin'] = np.sin(2 * np.pi * times.dayofyear / 365.25)
    feat['doy_cos'] = np.cos(2 * np.pi * times.dayofyear / 365.25)
    feat['month_sin'] = np.sin(2 * np.pi * times.month / 12)
    feat['month_cos'] = np.cos(2 * np.pi * times.month / 12)
    feat['dow_sin'] = np.sin(2 * np.pi * times.dayofweek / 7)
    feat['dow_cos'] = np.cos(2 * np.pi * times.dayofweek / 7)

    # --- Feature lag (produzione storica) ---
    for lag in [1, 2, 3, 6, 12, 24, 48]:
        feat[f'power_lag_{lag}h'] = feat['power'].shift(lag)

    # Rolling stats (finestra 6h, 24h)
    feat['power_roll6h_mean'] = feat['power'].rolling(6, min_periods=1).mean()
    feat['power_roll24h_mean'] = feat['power'].rolling(24, min_periods=1).mean()
    feat['ghi_roll3h_mean'] = feat['ghi'].rolling(3, min_periods=1).mean()

    # --- Ramping feature ---
    feat['power_delta_1h'] = feat['power'].diff(1)
    feat['ghi_delta_1h'] = feat['ghi'].diff(1)

    # Maschera ore notturne (sun elevation < 5 gradi)
    feat['is_daytime'] = (feat['solar_elevation'] > 5).astype(int)

    return feat.dropna()


def prepare_sequences(features: np.ndarray,
                       targets: np.ndarray,
                       lookback: int = 48,
                       horizon: int = 24) -> tuple:
    """
    Crea sequenze input/output per LSTM.

    Returns:
        X: (N, lookback, n_features)
        y: (N, horizon)
    """
    X, y = [], []
    for i in range(lookback, len(features) - horizon + 1):
        X.append(features[i - lookback:i])
        y.append(targets[i:i + horizon])
    return np.array(X, dtype=np.float32), np.array(y, dtype=np.float32)

import numpy as np
import pandas as pd
from scipy.stats import weibull_min

def create_wind_features(df: pd.DataFrame,
                          hub_height: float = 100,
                          ref_height: float = 10,
                          rated_power: float = 2000,  # kW
                          cut_in: float = 3.0,        # m/s
                          rated_speed: float = 12.0,  # m/s
                          cut_out: float = 25.0       # m/s
                          ) -> pd.DataFrame:
    """
    Feature engineering per forecasting eolico.

    df deve contenere:
      - wind_speed_10m: velocità vento a 10m [m/s]
      - wind_dir_10m: direzione vento [gradi 0-360]
      - temp_air: temperatura aria [degC]
      - pressure: pressione atmosferica [Pa]
      - roughness_length: lunghezza rugosita terreno [m]
    """
    feat = df.copy()

    # --- Wind shear: extrapolazione all'altezza hub ---
    # Profilo logaritmico (più accurato del power law per terreni complessi)
    z0 = feat.get('roughness_length', 0.03)  # default terreno agricolo
    feat['wind_speed_hub'] = (feat['wind_speed_10m'] *
                               np.log(hub_height / z0) /
                               np.log(ref_height / z0))

    # Power law come alternativa
    shear_exp = 0.14  # terreno aperto standard
    feat['wind_speed_hub_powerlaw'] = (feat['wind_speed_10m'] *
                                        (hub_height / ref_height) ** shear_exp)

    # --- Densita aria (funzione di T e P) ---
    R = 287.05  # costante gas reale aria secca [J/kg/K]
    feat['air_density'] = feat['pressure'] / (R * (feat['temp_air'] + 273.15))

    # Velocita normalizzata per densita (v_equiv = v * (rho/rho_ref)^(1/3))
    rho_ref = 1.225  # kg/m3 a 15degC, 1013.25 hPa
    feat['wind_speed_density_corrected'] = (feat['wind_speed_hub'] *
        (feat['air_density'] / rho_ref) ** (1/3))

    # --- Curva di potenza teorica ---
    v = feat['wind_speed_density_corrected'].values

    def power_curve(v: np.ndarray) -> np.ndarray:
        """Curva di potenza cubica con saturazione."""
        p = np.zeros_like(v)
        # Zona cubica: cut_in <= v < rated_speed
        mask_cubic = (v >= cut_in) & (v < rated_speed)
        p[mask_cubic] = rated_power * ((v[mask_cubic] - cut_in) /
                                        (rated_speed - cut_in)) ** 3
        # Zona nominale: rated_speed <= v < cut_out
        mask_rated = (v >= rated_speed) & (v < cut_out)
        p[mask_rated] = rated_power
        # Zona spenta: v < cut_in o v >= cut_out
        return p

    feat['theoretical_power'] = power_curve(v)
    feat['power_ratio'] = feat['theoretical_power'] / rated_power

    # --- Componenti vettoriali del vento ---
    wind_dir_rad = np.radians(feat['wind_dir_10m'])
    feat['wind_u'] = -feat['wind_speed_hub'] * np.sin(wind_dir_rad)
    feat['wind_v'] = -feat['wind_speed_hub'] * np.cos(wind_dir_rad)

    # --- Turbolenza (Turbulence Intensity) ---
    feat['turbulence_intensity'] = (feat['wind_speed_hub'].rolling(10).std() /
                                     feat['wind_speed_hub'].rolling(10).mean().clip(lower=0.1))

    # --- Weibull: parametri stimati su finestra mobile ---
    def estimate_weibull_shape(series: pd.Series, min_count: int = 20) -> float:
        """Stima parametro di forma k della distribuzione Weibull."""
        vals = series.dropna().values
        if len(vals) < min_count:
            return np.nan
        try:
            params = weibull_min.fit(vals[vals > 0], floc=0)
            return params[0]  # k (shape parameter)
        except Exception:
            return np.nan

    feat['weibull_k_72h'] = (feat['wind_speed_hub']
                              .rolling(72, min_periods=20)
                              .apply(estimate_weibull_shape, raw=False))

    # --- Encoding temporale ciclico ---
    times = feat.index
    feat['hour_sin'] = np.sin(2 * np.pi * times.hour / 24)
    feat['hour_cos'] = np.cos(2 * np.pi * times.hour / 24)
    feat['doy_sin'] = np.sin(2 * np.pi * times.dayofyear / 365.25)
    feat['doy_cos'] = np.cos(2 * np.pi * times.dayofyear / 365.25)
    feat['month_sin'] = np.sin(2 * np.pi * times.month / 12)
    feat['month_cos'] = np.cos(2 * np.pi * times.month / 12)

    # --- Lag e rolling per potenza storica ---
    for lag in [1, 2, 3, 6, 12, 24, 48, 168]:  # 168h = 1 settimana
        feat[f'power_lag_{lag}h'] = feat['power'].shift(lag)

    feat['power_roll6h_std'] = feat['power'].rolling(6).std()
    feat['power_roll24h_mean'] = feat['power'].rolling(24).mean()
    feat['wind_roll3h_mean'] = feat['wind_speed_hub'].rolling(3).mean()

    # --- Wake effect proxy (vento in direzione impianto vicino) ---
    # Semplificato: penalizzazione direzione principale
    feat['wake_factor'] = np.where(
        (feat['wind_dir_10m'] > 170) & (feat['wind_dir_10m'] < 220),
        0.85,  # direzione dominante = wake effect stimato 15%
        1.0
    )

    return feat.dropna()

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import numpy as np
from typing import Optional


class SolarEncoder(nn.Module):
    """Encoder LSTM per sequenze storiche di input."""

    def __init__(self, input_size: int, hidden_size: int, num_layers: int,
                 dropout: float = 0.2):
        super().__init__()
        self.hidden_size = hidden_size
        self.num_layers = num_layers

        self.lstm = nn.LSTM(
            input_size=input_size,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=dropout if num_layers > 1 else 0.0,
            bidirectional=False
        )
        self.layer_norm = nn.LayerNorm(hidden_size)

    def forward(self, x: torch.Tensor):
        # x: (batch, seq_len, input_size)
        outputs, (hidden, cell) = self.lstm(x)
        # outputs: (batch, seq_len, hidden_size) - tutti gli hidden state
        outputs = self.layer_norm(outputs)
        return outputs, hidden, cell


class BahdanauAttention(nn.Module):
    """Meccanismo di attention Bahdanau (additive attention)."""

    def __init__(self, hidden_size: int):
        super().__init__()
        self.W_q = nn.Linear(hidden_size, hidden_size, bias=False)
        self.W_k = nn.Linear(hidden_size, hidden_size, bias=False)
        self.v = nn.Linear(hidden_size, 1, bias=False)

    def forward(self, query: torch.Tensor,
                keys: torch.Tensor) -> tuple:
        # query: (batch, 1, hidden_size) - decoder hidden state
        # keys: (batch, seq_len, hidden_size) - encoder outputs
        q = self.W_q(query)        # (batch, 1, hidden_size)
        k = self.W_k(keys)         # (batch, seq_len, hidden_size)
        scores = self.v(torch.tanh(q + k))  # (batch, seq_len, 1)
        weights = torch.softmax(scores, dim=1)  # (batch, seq_len, 1)
        context = (weights * keys).sum(dim=1, keepdim=True)  # (batch, 1, hidden)
        return context, weights.squeeze(-1)


class SolarDecoder(nn.Module):
    """Decoder LSTM con attention per forecast multi-step."""

    def __init__(self, input_size: int, hidden_size: int, num_layers: int,
                 dropout: float = 0.2):
        super().__init__()
        self.hidden_size = hidden_size
        self.num_layers = num_layers

        self.attention = BahdanauAttention(hidden_size)
        self.lstm = nn.LSTM(
            input_size=input_size + hidden_size,  # input + context
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=dropout if num_layers > 1 else 0.0
        )
        self.output_proj = nn.Sequential(
            nn.Linear(hidden_size, hidden_size // 2),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(hidden_size // 2, 1)
        )
        self.layer_norm = nn.LayerNorm(hidden_size)

    def forward(self, x: torch.Tensor,
                encoder_outputs: torch.Tensor,
                hidden: torch.Tensor,
                cell: torch.Tensor) -> tuple:
        # x: (batch, 1, input_size) - input corrente del decoder
        context, attn_weights = self.attention(
            hidden[-1:].permute(1, 0, 2),  # (batch, 1, hidden)
            encoder_outputs
        )
        lstm_input = torch.cat([x, context], dim=-1)
        output, (hidden, cell) = self.lstm(lstm_input, (hidden, cell))
        output = self.layer_norm(output)
        prediction = self.output_proj(output)  # (batch, 1, 1)
        return prediction.squeeze(-1), hidden, cell, attn_weights


class SolarForecastModel(nn.Module):
    """Modello encoder-decoder completo per previsione solare."""

    def __init__(self,
                 encoder_features: int,    # numero feature input storiche
                 decoder_features: int,    # numero feature future note (NWP)
                 hidden_size: int = 256,
                 num_layers: int = 3,
                 horizon: int = 24,
                 dropout: float = 0.2):
        super().__init__()
        self.horizon = horizon
        self.hidden_size = hidden_size

        self.encoder = SolarEncoder(encoder_features, hidden_size,
                                     num_layers, dropout)
        self.decoder = SolarDecoder(decoder_features, hidden_size,
                                     num_layers, dropout)

        # Proiezione per adattare hidden encoder a decoder (se layers != 1)
        self.hidden_proj = nn.Linear(hidden_size, hidden_size)

    def forward(self, encoder_input: torch.Tensor,
                decoder_input: torch.Tensor,
                teacher_forcing_ratio: float = 0.5) -> torch.Tensor:
        """
        Args:
            encoder_input: (batch, lookback, encoder_features) - storia
            decoder_input: (batch, horizon, decoder_features) - NWP futuro
            teacher_forcing_ratio: prob. di usare ground truth nel training

        Returns:
            predictions: (batch, horizon)
        """
        batch_size = encoder_input.size(0)

        # Encoding della storia
        encoder_outputs, hidden, cell = self.encoder(encoder_input)

        # Proiezione hidden state
        hidden = self.hidden_proj(hidden)

        predictions = torch.zeros(batch_size, self.horizon,
                                   device=encoder_input.device)

        # Decoding step-by-step
        decoder_in = decoder_input[:, 0:1, :]  # primo passo

        for t in range(self.horizon):
            pred, hidden, cell, _ = self.decoder(
                decoder_in, encoder_outputs, hidden, cell
            )
            predictions[:, t] = pred.squeeze(-1)

            # Teacher forcing: usa ground truth o previsione precedente
            if t < self.horizon - 1:
                decoder_in = decoder_input[:, t+1:t+2, :]

        return predictions


def train_solar_model(
        model: SolarForecastModel,
        train_loader: DataLoader,
        val_loader: DataLoader,
        epochs: int = 100,
        lr: float = 1e-3,
        patience: int = 15,
        device: str = 'cuda' if torch.cuda.is_available() else 'cpu'
) -> dict:
    """
    Training loop completo con early stopping e LR scheduler.

    Returns:
        dict con history di loss e best_epoch
    """
    model = model.to(device)
    optimizer = optim.AdamW(model.parameters(), lr=lr, weight_decay=1e-4)
    criterion = nn.HuberLoss(delta=1.0)  # più robusto a outlier vs MSE

    # Cosine Annealing con Warm Restarts
    scheduler = optim.lr_scheduler.CosineAnnealingWarmRestarts(
        optimizer, T_0=20, T_mult=2, eta_min=1e-6
    )

    history = {'train_loss': [], 'val_loss': [], 'lr': []}
    best_val_loss = float('inf')
    patience_counter = 0
    best_state = None

    for epoch in range(epochs):
        # Training
        model.train()
        train_loss = 0.0
        for batch in train_loader:
            enc_in, dec_in, targets = [b.to(device) for b in batch]
            optimizer.zero_grad()

            # Teacher forcing decay: da 0.8 a 0.1 nel corso del training
            tf_ratio = max(0.1, 0.8 - epoch * 0.007)
            preds = model(enc_in, dec_in, tf_ratio)

            loss = criterion(preds, targets)
            loss.backward()
            # Gradient clipping per stabilità
            torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
            optimizer.step()
            train_loss += loss.item()

        # Validation
        model.eval()
        val_loss = 0.0
        with torch.no_grad():
            for batch in val_loader:
                enc_in, dec_in, targets = [b.to(device) for b in batch]
                preds = model(enc_in, dec_in, teacher_forcing_ratio=0.0)
                val_loss += criterion(preds, targets).item()

        train_loss /= len(train_loader)
        val_loss /= len(val_loader)
        scheduler.step()

        current_lr = optimizer.param_groups[0]['lr']
        history['train_loss'].append(train_loss)
        history['val_loss'].append(val_loss)
        history['lr'].append(current_lr)

        print(f"Epoch {epoch+1:3d} | Train: {train_loss:.4f} | "
              f"Val: {val_loss:.4f} | LR: {current_lr:.2e}")

        # Early stopping
        if val_loss < best_val_loss - 1e-4:
            best_val_loss = val_loss
            patience_counter = 0
            best_state = {k: v.clone() for k, v in model.state_dict().items()}
        else:
            patience_counter += 1
            if patience_counter >= patience:
                print(f"Early stopping at epoch {epoch+1}")
                break

    # Ripristina il miglior modello
    if best_state:
        model.load_state_dict(best_state)
    history['best_epoch'] = epoch + 1 - patience_counter

    return history

import torch
import torch.nn as nn
import torch.nn.functional as F
from typing import List


class CausalConv1d(nn.Module):
    """Convoluzione causale: usa solo input passati, mai futuri."""

    def __init__(self, in_channels: int, out_channels: int,
                 kernel_size: int, dilation: int = 1):
        super().__init__()
        self.dilation = dilation
        self.kernel_size = kernel_size
        # Padding causale: solo a sinistra
        self.padding = (kernel_size - 1) * dilation
        self.conv = nn.Conv1d(in_channels, out_channels,
                               kernel_size=kernel_size,
                               padding=self.padding,
                               dilation=dilation)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (batch, channels, seq_len)
        out = self.conv(x)
        # Rimuovi padding aggiunto a destra per causality
        if self.padding > 0:
            out = out[:, :, :-self.padding]
        return out


class TCNResidualBlock(nn.Module):
    """Blocco residuale TCN con dilated causal convolutions."""

    def __init__(self, in_channels: int, out_channels: int,
                 kernel_size: int, dilation: int, dropout: float = 0.1):
        super().__init__()
        self.conv1 = CausalConv1d(in_channels, out_channels,
                                   kernel_size, dilation)
        self.conv2 = CausalConv1d(out_channels, out_channels,
                                   kernel_size, dilation)
        self.bn1 = nn.BatchNorm1d(out_channels)
        self.bn2 = nn.BatchNorm1d(out_channels)
        self.dropout = nn.Dropout(dropout)
        self.activation = nn.GELU()

        # Skip connection con 1x1 conv se canali diversi
        self.skip = (nn.Conv1d(in_channels, out_channels, 1)
                     if in_channels != out_channels else nn.Identity())

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        residual = self.skip(x)
        out = self.activation(self.bn1(self.conv1(x)))
        out = self.dropout(out)
        out = self.activation(self.bn2(self.conv2(out)))
        out = self.dropout(out)
        return self.activation(out + residual)


class TCNWindForecast(nn.Module):
    """
    Temporal Convolutional Network per previsione eolica multi-step.

    Architettura: Input -> TCN blocks (dilatazione esponenziale)
                          -> Global Average Pool -> FC -> Output
    """

    def __init__(self,
                 input_size: int,
                 num_channels: List[int],
                 kernel_size: int = 3,
                 horizon: int = 24,
                 dropout: float = 0.1):
        """
        Args:
            input_size: numero di feature in input
            num_channels: canali per ogni livello TCN, es. [64, 128, 256, 256]
            kernel_size: dimensione kernel (3 o 5 tipici)
            horizon: passi di previsione futura
            dropout: tasso di dropout
        """
        super().__init__()
        self.horizon = horizon

        # Input projection
        self.input_proj = nn.Conv1d(input_size, num_channels[0], kernel_size=1)

        # TCN blocks con dilatazione esponenziale
        layers = []
        for i, (in_ch, out_ch) in enumerate(
                zip(num_channels[:-1], num_channels[1:])):
            dilation = 2 ** i  # 1, 2, 4, 8, 16, ...
            layers.append(TCNResidualBlock(in_ch, out_ch, kernel_size,
                                            dilation, dropout))
        self.tcn = nn.Sequential(*layers)

        # Output head
        final_channels = num_channels[-1]
        self.output_head = nn.Sequential(
            nn.AdaptiveAvgPool1d(1),   # global average pooling temporale
            nn.Flatten(),
            nn.Linear(final_channels, final_channels // 2),
            nn.GELU(),
            nn.Dropout(dropout),
            nn.Linear(final_channels // 2, horizon)
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # x: (batch, seq_len, input_size) -> porta in (batch, input_size, seq_len)
        x = x.permute(0, 2, 1)
        x = self.input_proj(x)    # (batch, ch[0], seq_len)
        x = self.tcn(x)           # (batch, ch[-1], seq_len)
        out = self.output_head(x) # (batch, horizon)
        return out


# Esempio utilizzo
def build_wind_model(n_features: int = 35, horizon: int = 24) -> TCNWindForecast:
    return TCNWindForecast(
        input_size=n_features,
        num_channels=[64, 128, 256, 256, 128],
        kernel_size=3,
        horizon=horizon,
        dropout=0.15
    )


# Calcolo del campo ricettivo TCN
def receptive_field(num_levels: int, kernel_size: int) -> int:
    """Campo ricettivo del TCN (quante ore passate 'vede')."""
    return 1 + 2 * (kernel_size - 1) * (2 ** num_levels - 1)

# Con 5 livelli e kernel 3: RF = 1 + 2*2*(32-1) = 125 ore ~5 giorni
print(f"Receptive field: {receptive_field(5, 3)} timesteps")

import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from typing import Dict, List
import lightgbm as lgb


class WeightedEnsembleForecaster:
    """
    Ensemble adattivo che combina previsioni da più modelli base.

    Strategia: weighted average con pesi ottimizzati via Ridge regression
    sul validation set, condizionati sul contesto meteorologico.
    """

    def __init__(self, model_names: List[str], n_quantiles: int = 4):
        self.model_names = model_names
        self.n_quantiles = n_quantiles
        # Meta-learner per ogni slot temporale (es. 24 ore)
        self.meta_learners: Dict = {}
        self.scalers: Dict = {}

    def fit(self, base_predictions: Dict[str, np.ndarray],
            targets: np.ndarray,
            context_features: np.ndarray) -> 'WeightedEnsembleForecaster':
        """
        Allena i meta-learner sull'insieme di validazione.

        Args:
            base_predictions: {model_name: array (N, horizon)}
            targets: array (N, horizon)
            context_features: array (N, n_context) - GHI, ora, stagione, ecc.
        """
        horizon = targets.shape[1]
        n_models = len(self.model_names)

        for h in range(horizon):
            # Stack previsioni base per passo h
            X = np.column_stack([
                base_predictions[m][:, h] for m in self.model_names
            ])  # (N, n_models)

            # Aggiungi feature di contesto
            X_ctx = np.hstack([X, context_features])  # (N, n_models + n_ctx)

            y = targets[:, h]

            # Scaler
            scaler = StandardScaler()
            X_scaled = scaler.fit_transform(X_ctx)
            self.scalers[h] = scaler

            # Ridge regression con coefficienti non negativi (weight positivi)
            # Usiamo Ridge standard + clip post-hoc
            meta = Ridge(alpha=1.0, fit_intercept=True)
            meta.fit(X_scaled, y)
            self.meta_learners[h] = meta

        return self

    def predict(self, base_predictions: Dict[str, np.ndarray],
                context_features: np.ndarray) -> np.ndarray:
        """
        Genera previsione ensemble.

        Returns:
            array (N, horizon)
        """
        horizon = len(self.meta_learners)
        N = list(base_predictions.values())[0].shape[0]
        ensemble_pred = np.zeros((N, horizon))

        for h in range(horizon):
            X = np.column_stack([
                base_predictions[m][:, h] for m in self.model_names
            ])
            X_ctx = np.hstack([X, context_features])
            X_scaled = self.scalers[h].transform(X_ctx)
            ensemble_pred[:, h] = self.meta_learners[h].predict(X_scaled)

        return ensemble_pred

    def get_model_weights(self, hour_idx: int) -> Dict[str, float]:
        """Restituisce i pesi normalizzati dei modelli base per un'ora."""
        meta = self.meta_learners[hour_idx]
        # Solo i pesi dei modelli base (esclude context features)
        raw_weights = meta.coef_[:len(self.model_names)]
        raw_weights = np.clip(raw_weights, 0, None)  # non negativi
        total = raw_weights.sum()
        if total > 0:
            norm_weights = raw_weights / total
        else:
            norm_weights = np.ones(len(self.model_names)) / len(self.model_names)
        return dict(zip(self.model_names, norm_weights))


# LightGBM come modello base alternativo a LSTM per contesti tabulari
def train_lgbm_solar(X_train: np.ndarray, y_train: np.ndarray,
                      X_val: np.ndarray, y_val: np.ndarray,
                      horizon: int = 24) -> List[lgb.Booster]:
    """Allena un modello LightGBM per ogni passo del forecast."""
    models = []
    for h in range(horizon):
        params = {
            'objective': 'huber',
            'metric': 'huber',
            'num_leaves': 63,
            'learning_rate': 0.05,
            'feature_fraction': 0.8,
            'bagging_fraction': 0.8,
            'bagging_freq': 5,
            'verbose': -1,
            'n_jobs': -1
        }
        train_data = lgb.Dataset(X_train, label=y_train[:, h])
        val_data = lgb.Dataset(X_val, label=y_val[:, h], reference=train_data)

        model = lgb.train(
            params, train_data,
            num_boost_round=500,
            valid_sets=[val_data],
            callbacks=[
                lgb.early_stopping(stopping_rounds=30),
                lgb.log_evaluation(period=100)
            ]
        )
        models.append(model)
    return models

import torch
import torch.nn as nn
import numpy as np
from typing import List


class PinballLoss(nn.Module):
    """
    Pinball loss (quantile loss) per quantile regression.
    Supporta multipli quantili simultaneamente.
    """

    def __init__(self, quantiles: List[float]):
        super().__init__()
        self.quantiles = torch.FloatTensor(quantiles)

    def forward(self, preds: torch.Tensor,
                targets: torch.Tensor) -> torch.Tensor:
        """
        Args:
            preds: (batch, horizon, n_quantiles)
            targets: (batch, horizon)

        Returns:
            loss scalare
        """
        q = self.quantiles.to(preds.device)  # (n_quantiles,)
        # Espandi targets: (batch, horizon, 1)
        y = targets.unsqueeze(-1)
        # errors: (batch, horizon, n_quantiles)
        errors = y - preds
        # Pinball: max(q*e, (q-1)*e)
        loss = torch.max(q * errors, (q - 1) * errors)
        return loss.mean()


class QuantileLSTM(nn.Module):
    """
    LSTM che produce simultaneamente previsioni per N quantili.
    Più efficiente che allenare N modelli separati.
    """

    def __init__(self,
                 input_size: int,
                 hidden_size: int,
                 num_layers: int,
                 quantiles: List[float],
                 horizon: int,
                 dropout: float = 0.2):
        super().__init__()
        self.quantiles = quantiles
        self.n_q = len(quantiles)
        self.horizon = horizon

        self.lstm = nn.LSTM(
            input_size=input_size,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=dropout if num_layers > 1 else 0.0
        )
        self.dropout = nn.Dropout(dropout)

        # Head per ogni quantile con shared trunk + separate heads
        self.shared_head = nn.Sequential(
            nn.Linear(hidden_size, hidden_size // 2),
            nn.LayerNorm(hidden_size // 2),
            nn.GELU()
        )

        # Output separato per ogni quantile (garantisce ordering)
        self.quantile_heads = nn.ModuleList([
            nn.Linear(hidden_size // 2, horizon)
            for _ in quantiles
        ])

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        Returns:
            (batch, horizon, n_quantiles) - quantili ORDINATI
        """
        lstm_out, (hidden, _) = self.lstm(x)
        # Usa l'ultimo hidden state
        last_hidden = self.dropout(hidden[-1])  # (batch, hidden_size)
        shared = self.shared_head(last_hidden)  # (batch, hidden//2)

        # Produce raw quantili
        raw_quantiles = torch.stack([
            head(shared) for head in self.quantile_heads
        ], dim=-1)  # (batch, horizon, n_quantiles)

        # Quantile crossing fix: ordina i quantili per monotonia
        # Tecnica: softplus cumsum per garantire q1 <= q2 <= ... <= qN
        ordered = torch.cumsum(
            F.softplus(raw_quantiles), dim=-1
        ) - F.softplus(raw_quantiles[..., :1])
        # Ancora al quantile mediano della prima head
        median_base = self.quantile_heads[len(self.quantiles)//2](shared)
        ordered = ordered + median_base.unsqueeze(-1)

        return ordered


def compute_coverage_and_sharpness(
        y_true: np.ndarray,
        q_low: np.ndarray,
        q_high: np.ndarray,
        nominal_level: float = 0.8) -> dict:
    """
    Valuta calibrazione (coverage) e sharpness degli intervalli predittivi.

    Args:
        y_true: valori reali (N,)
        q_low: quantile inferiore (N,) - es. q10
        q_high: quantile superiore (N,) - es. q90
        nominal_level: livello nominale (es. 0.80 per 10%-90%)

    Returns:
        {'coverage': float, 'sharpness': float, 'winkler_score': float}
    """
    coverage = np.mean((y_true >= q_low) & (y_true <= q_high))
    sharpness = np.mean(q_high - q_low)  # larghezza media intervallo

    # Winkler Score: penalizza sia scarsa coverage che larghezza eccessiva
    alpha = 1 - nominal_level  # 0.20 per 80% PI
    winkler = sharpness.copy() if hasattr(sharpness, 'copy') else sharpness
    undercover = y_true < q_low
    overcover = y_true > q_high
    winkler_arr = (q_high - q_low).copy()
    winkler_arr[undercover] += (2 / alpha) * (q_low[undercover] - y_true[undercover])
    winkler_arr[overcover] += (2 / alpha) * (y_true[overcover] - q_high[overcover])

    return {
        'coverage': float(coverage),
        'nominal_level': nominal_level,
        'sharpness_mean_kw': float(sharpness),
        'winkler_score': float(np.mean(winkler_arr)),
        'calibration_error': float(abs(coverage - nominal_level))
    }

"""
Pipeline di produzione per forecasting rinnovabili.
Stack: MLflow (tracking + registry) + FastAPI (serving) + APScheduler (scheduling)
"""
import mlflow
import mlflow.pytorch
import torch
import numpy as np
import pandas as pd
from fastapi import FastAPI, HTTPException, BackgroundTasks
from pydantic import BaseModel, Field
from apscheduler.schedulers.background import BackgroundScheduler
from datetime import datetime, timedelta
import logging
from typing import Optional, List
import httpx

logger = logging.getLogger(__name__)

# ---- MLflow Experiment Setup ----

MLFLOW_TRACKING_URI = "http://mlflow-server:5000"
EXPERIMENT_NAME = "solar_forecast_puglia_5mw"

mlflow.set_tracking_uri(MLFLOW_TRACKING_URI)
mlflow.set_experiment(EXPERIMENT_NAME)


def log_training_run(model, metrics: dict, params: dict,
                     tags: dict) -> str:
    """Logga training run su MLflow e registra il modello."""
    with mlflow.start_run() as run:
        # Log parametri
        mlflow.log_params(params)
        # Log metriche
        mlflow.log_metrics(metrics)
        # Log tags
        mlflow.set_tags(tags)

        # Salva modello PyTorch
        mlflow.pytorch.log_model(
            pytorch_model=model,
            artifact_path="model",
            registered_model_name="solar_forecast_encoder_decoder"
        )

        run_id = run.info.run_id
        logger.info(f"MLflow run {run_id} logged successfully")
        return run_id


def load_production_model(model_name: str,
                           stage: str = "Production") -> tuple:
    """Carica modello dalla MLflow Model Registry."""
    client = mlflow.tracking.MlflowClient()
    versions = client.get_latest_versions(model_name, stages=[stage])

    if not versions:
        raise ValueError(f"No model in stage {stage} for {model_name}")

    model_version = versions[0]
    model_uri = f"models:/{model_name}/{stage}"
    model = mlflow.pytorch.load_model(model_uri)
    model.eval()

    return model, model_version


# ---- FastAPI Serving ----

app = FastAPI(
    title="Renewable Energy Forecast API",
    description="Day-ahead and intraday power forecast for solar/wind plants",
    version="1.0.0"
)


class ForecastRequest(BaseModel):
    plant_id: str = Field(..., description="ID impianto (es. 'solar_puglia_5mw')")
    forecast_date: str = Field(..., description="Data previsione YYYY-MM-DD")
    horizon_hours: int = Field(default=24, ge=1, le=72,
                                description="Orizzonte forecast in ore")
    include_quantiles: bool = Field(default=True,
                                     description="Includi intervalli probabilistici")


class ForecastResponse(BaseModel):
    plant_id: str
    forecast_date: str
    generated_at: str
    model_version: str
    forecasts: List[dict]  # {timestamp, power_kw, q10, q25, q75, q90}
    metrics: Optional[dict] = None


# Cache modello in memoria
_model_cache: dict = {}


def get_model(plant_id: str):
    """Lazy loading del modello con cache."""
    if plant_id not in _model_cache:
        model_name = f"forecast_{plant_id}"
        model, version = load_production_model(model_name)
        _model_cache[plant_id] = {'model': model, 'version': version}
    return _model_cache[plant_id]


async def fetch_nwp_data(lat: float, lon: float,
                          start_date: str, hours: int) -> pd.DataFrame:
    """Scarica previsioni NWP da Open-Meteo API."""
    params = {
        "latitude": lat,
        "longitude": lon,
        "hourly": "shortwave_radiation,direct_radiation,diffuse_radiation,"
                  "temperature_2m,wind_speed_10m,wind_direction_10m,"
                  "cloud_cover,precipitation",
        "start_date": start_date,
        "end_date": start_date,
        "forecast_days": max(1, hours // 24 + 1),
        "models": "best_match"
    }
    async with httpx.AsyncClient() as client:
        resp = await client.get("https://api.open-meteo.com/v1/forecast",
                                 params=params, timeout=30)
        resp.raise_for_status()
        data = resp.json()

    # Converti in DataFrame
    hourly = data['hourly']
    df = pd.DataFrame(hourly)
    df['timestamp'] = pd.to_datetime(df['time'])
    df = df.set_index('timestamp').drop(columns=['time'])
    return df


@app.post("/forecast", response_model=ForecastResponse)
async def generate_forecast(request: ForecastRequest):
    """Endpoint principale per generare previsioni."""
    try:
        # Configurazione impianto (in produzione: da database)
        plant_config = {
            'solar_puglia_5mw': {'lat': 40.5, 'lon': 17.2,
                                   'rated_power': 5000, 'type': 'solar'},
            'wind_campania_12mw': {'lat': 40.8, 'lon': 15.1,
                                    'rated_power': 12000, 'type': 'wind'}
        }

        if request.plant_id not in plant_config:
            raise HTTPException(404, f"Plant {request.plant_id} not found")

        config = plant_config[request.plant_id]
        cached = get_model(request.plant_id)
        model = cached['model']
        model_version = str(cached['version'].version)

        # Fetch dati NWP
        nwp_data = await fetch_nwp_data(
            config['lat'], config['lon'],
            request.forecast_date, request.horizon_hours
        )

        # Preprocessing e feature engineering (semplificato)
        # In produzione: pipeline standardizzata e validata
        features = preprocess_for_inference(nwp_data, config)

        # Inference
        with torch.no_grad():
            enc_input = torch.FloatTensor(features['encoder']).unsqueeze(0)
            dec_input = torch.FloatTensor(features['decoder']).unsqueeze(0)

            if request.include_quantiles:
                predictions = model.predict_quantiles(enc_input, dec_input)
                # predictions: (1, horizon, 5) per q10,q25,q50,q75,q90
                preds_np = predictions.squeeze(0).numpy()
            else:
                predictions = model(enc_input, dec_input)
                preds_np = predictions.squeeze(0).numpy()

        # Costruisci risposta
        start_ts = pd.Timestamp(request.forecast_date)
        forecasts = []
        for h in range(min(request.horizon_hours, len(preds_np))):
            ts = start_ts + pd.Timedelta(hours=h)
            entry = {
                'timestamp': ts.isoformat(),
                'power_kw': max(0, float(preds_np[h, 2] if request.include_quantiles
                                         else preds_np[h]))
            }
            if request.include_quantiles:
                entry.update({
                    'q10': max(0, float(preds_np[h, 0])),
                    'q25': max(0, float(preds_np[h, 1])),
                    'q75': max(0, float(preds_np[h, 3])),
                    'q90': max(0, float(preds_np[h, 4]))
                })
            forecasts.append(entry)

        return ForecastResponse(
            plant_id=request.plant_id,
            forecast_date=request.forecast_date,
            generated_at=datetime.utcnow().isoformat(),
            model_version=model_version,
            forecasts=forecasts
        )

    except Exception as e:
        logger.error(f"Forecast error for {request.plant_id}: {e}",
                     exc_info=True)
        raise HTTPException(500, f"Forecast generation failed: {str(e)}")


def preprocess_for_inference(nwp_data: pd.DataFrame,
                              config: dict) -> dict:
    """Preprocessing NWP per inference (placeholder)."""
    # In produzione: replicare esattamente la pipeline di training
    n_enc = 48  # lookback 48 ore
    n_dec = 24  # horizon 24 ore
    n_enc_feat = 25
    n_dec_feat = 10

    return {
        'encoder': np.random.randn(n_enc, n_enc_feat).astype(np.float32),
        'decoder': np.random.randn(n_dec, n_dec_feat).astype(np.float32)
    }


# ---- Scheduler per forecast automatici ----

scheduler = BackgroundScheduler()


@app.on_event("startup")
async def startup_event():
    """Avvia scheduler al boot dell'applicazione."""
    # Genera forecast day-ahead ogni giorno alle 11:00 (dopo asta MGP)
    scheduler.add_job(
        trigger_daily_forecast,
        'cron', hour=11, minute=0,
        id='day_ahead_forecast',
        replace_existing=True
    )
    # Aggiorna forecast intraday ogni 4 ore
    scheduler.add_job(
        trigger_intraday_update,
        'cron', hour='*/4',
        id='intraday_update',
        replace_existing=True
    )
    scheduler.start()
    logger.info("Forecast scheduler started")


def trigger_daily_forecast():
    """Genera forecast day-ahead per tutti gli impianti."""
    tomorrow = (datetime.now() + timedelta(days=1)).strftime('%Y-%m-%d')
    plants = ['solar_puglia_5mw', 'wind_campania_12mw']
    for plant_id in plants:
        try:
            logger.info(f"Generating day-ahead forecast for {plant_id} on {tomorrow}")
            # In produzione: chiamata asincrona all'endpoint /forecast
        except Exception as e:
            logger.error(f"Failed forecast for {plant_id}: {e}")

from evidently import ColumnMapping
from evidently.report import Report
from evidently.metric_preset import RegressionPreset, DataDriftPreset
import pandas as pd
import numpy as np
from datetime import datetime, timedelta


class ForecastMonitor:
    """
    Monitoring continuo delle performance di forecasting.
    Usa Evidently per drift detection e report automatici.
    """

    def __init__(self, baseline_window_days: int = 30,
                 alert_rmse_threshold: float = 0.15,  # nRMSE > 15% = alert
                 alert_drift_threshold: float = 0.3):
        self.baseline_window_days = baseline_window_days
        self.alert_rmse_threshold = alert_rmse_threshold
        self.alert_drift_threshold = alert_drift_threshold
        self.performance_history: list = []

    def evaluate_forecast(self,
                           y_true: np.ndarray,
                           y_pred: np.ndarray,
                           rated_power: float,
                           timestamp: datetime) -> dict:
        """Calcola metriche e aggiunge alla storia."""
        rmse = np.sqrt(np.mean((y_true - y_pred) ** 2))
        mae = np.mean(np.abs(y_true - y_pred))
        nrmse = rmse / rated_power

        # Skill score vs persistence
        y_persist = np.roll(y_true, 24)  # previsione persistence 24h
        rmse_persist = np.sqrt(np.mean((y_true[24:] - y_persist[24:]) ** 2))
        skill = 1 - rmse / rmse_persist if rmse_persist > 0 else 0

        metrics = {
            'timestamp': timestamp.isoformat(),
            'rmse': float(rmse),
            'mae': float(mae),
            'nrmse': float(nrmse),
            'skill_score': float(skill),
            'alert': nrmse > self.alert_rmse_threshold
        }
        self.performance_history.append(metrics)

        if metrics['alert']:
            logger.warning(
                f"PERFORMANCE ALERT at {timestamp}: "
                f"nRMSE={nrmse:.3f} > threshold {self.alert_rmse_threshold}"
            )

        return metrics

    def generate_drift_report(self,
                               reference_features: pd.DataFrame,
                               current_features: pd.DataFrame,
                               output_path: str = "drift_report.html") -> Report:
        """Genera report Evidently per data drift."""
        column_mapping = ColumnMapping(
            target="power",
            prediction="forecast",
            numerical_features=reference_features.select_dtypes(
                include=[np.number]).columns.tolist()
        )
        report = Report(metrics=[
            DataDriftPreset(),
            RegressionPreset()
        ])
        report.run(
            reference_data=reference_features,
            current_data=current_features,
            column_mapping=column_mapping
        )
        report.save_html(output_path)
        logger.info(f"Drift report saved to {output_path}")
        return report

    def should_retrain(self) -> bool:
        """Determina se e necessario un retraining."""
        if len(self.performance_history) < 7:
            return False
        recent = self.performance_history[-7:]  # ultima settimana
        avg_nrmse = np.mean([m['nrmse'] for m in recent])
        n_alerts = sum(1 for m in recent if m['alert'])

        return avg_nrmse > self.alert_rmse_threshold * 1.2 or n_alerts >= 3