ghi_df.describe()

ghi_df.info()

# We execpt at least 525600 datapoints if data are recorded continuously for 3 years
(
    ghi_df
    .groupby("year")
    .size()
)

groupers = [
    "year",
    "month",
    "day_week",
    "hour"
]
fig, axs = plt.subplots(2, 2, figsize=(10, 5), sharey=True)

for ax, grouper in zip(axs.flatten(), groupers):

    sns.boxplot(
        data=ghi_df,
        x=grouper,
        y="ghi",
        showfliers=False,
        ax=ax
    )
    ax.grid(alpha=0.5)

plt.tight_layout()
plt.show()

dayly_stats = (
    ghi_df
    .groupby("day_year")
    .agg(
        mean=("ghi", "mean"),
        std_err_mean=("ghi", lambda x: sem(x))
    )
    .reset_index()
    
)

plt.figure(figsize=(8, 8))
plt.scatter(
    x=dayly_stats["day_year"].values,
    y=dayly_stats["mean"].values,
    s=3,
    c=dayly_stats["mean"].values,
    cmap="magma"
)
plt.errorbar(
    x=dayly_stats["day_year"].values,
    y=dayly_stats["mean"].values,
    yerr=dayly_stats["std_err_mean"].values * 1.96,
    linewidth=0.5,
    ls='none',
    c="grey"

)
plt.ylabel("Mean GHI")
plt.xlabel("Day of Year")
plt.grid()

single_year_ghi_df = ghi_df[ghi_df["year"] == 2015]
single_month_ghi_df = ghi_df[(ghi_df["year"] == 2015) & (ghi_df["month"] == 7)]
single_day_ghi_df = ghi_df[(ghi_df["year"] == 2015) & (ghi_df["day_year"] == 175)]

fig = plt.figure(figsize=(10, 6), tight_layout=True)
gs = gridspec.GridSpec(2, 2)

ax_year = fig.add_subplot(gs[0, :])
ax_month = fig.add_subplot(gs[1, 0])
ax_day = fig.add_subplot(gs[1, 1])


axes_dfs = {
   "Year": (ax_year, single_year_ghi_df),
   "Month": (ax_month, single_month_ghi_df),
   "Day": (ax_day, single_day_ghi_df)

}

for title, (ax, df) in axes_dfs.items():

    ax = visualize_time_series(
        ax=ax,
        time_series_name="Minutely GHI",
        time_series=df.set_index("time")["ghi"],
        s=0.05,
        alpha=0.5
    )
    ax.set_title(title)
    ax.legend()

satellite_df.head(5)

satellite_df.describe()

satellite_df.info()

#| code-fold: true
#| code-summary: "Show supplementary code"
fig, axs = plt.subplots(10, 10, figsize=(10, 10), sharex=True, sharey=True)

for row, ax in enumerate(axs.flatten()):

    satellite_img = satellite_df.loc[row, :][np.arange(1, 101)].values.astype("float")
    satellite_img = satellite_img.reshape(10, 10)

    ax.imshow(satellite_img, cmap="binary")

    ax.set_yticks([])
    ax.set_xticks([])

plt.tight_layout()
plt.show()

from modules.models.baselines import static_stat_forecast

ground_truth = []
forecasts = []
early_stop_counter = 0

for insample_indices, outsample_indices in tqdm(TIME_SERIES_SPLITTER.split(FULL_SERIES)):

    insample_series, outsample_series = FULL_SERIES[insample_indices], FULL_SERIES[outsample_indices]
    insample_forecast, outsample_forecast = static_stat_forecast(
        time_series=insample_series,
        steps_ahead=OUTSAMPLE_WINDOW_SIZE,
        stat_func=np.mean
    )

    forecasts.append(outsample_forecast)
    ground_truth.append(outsample_series)
    early_stop_counter += 1

    if early_stop_counter > EARLY_STOP_LIMIT:
        break

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = (
    pd.Series(np.hstack(ground_truth) - np.hstack(forecasts))
)
residual.index = outsample_ghi_df.index[:len(np.hstack(ground_truth))]

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
error_mean, error_sem = compute_rmse(
    y_true=np.array(ground_truth),
    y_forecast=np.array(forecasts),
    axis=0
)
MODELS_RESULTS["Constant Mean"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Constant Mean",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = outsample_forecast["ghi"] - outsample_forecast["ghi_forecast"]

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
y_forecast = np.array(
    np.split(
        outsample_forecast["ghi_forecast"].values,      
        int(len(outsample_forecast["ghi"].values) / 30)
    )
)
y_true = np.array(
    np.split(
        outsample_forecast["ghi"].values,      
        int(len(outsample_forecast["ghi"].values) / 30)
    )
)

error_mean, error_sem = compute_rmse(
    y_true=y_true,
    y_forecast=y_forecast,
    axis=0
)
MODELS_RESULTS["Mean Persistence"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Mean Persistence",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

from modules.models.baselines import ewm_forecast

ground_truth = []
forecasts = []
early_stop_counter = 0

for insample_indices, outsample_indices in tqdm(TIME_SERIES_SPLITTER.split(FULL_SERIES)):

    insample_series, outsample_series = FULL_SERIES[insample_indices], FULL_SERIES[outsample_indices]
    insample_forecast, outsample_forecast = ewm_forecast(
        time_series=insample_series,
        steps_ahead=30,
        alpha=0.5,
    )

    forecasts.append(outsample_forecast)
    ground_truth.append(outsample_series)
    early_stop_counter += 1

    if early_stop_counter > EARLY_STOP_LIMIT:
        break

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = (
    pd.Series(np.hstack(ground_truth) - np.hstack(forecasts))
)
residual.index = outsample_ghi_df.index[:len(np.hstack(ground_truth))]

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
error_mean, error_sem = compute_rmse(
    y_true=np.array(ground_truth),
    y_forecast=np.array(forecasts),
    axis=0
)
MODELS_RESULTS["Exponential Moving Average"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Exponential Moving Average",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=weighting_df["weight_day_to_minute"],
    time_series_name="Weighting Day-Minute",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=weighting_df["weight_day_to_minute"][-(60 * 24 * 30):],
    time_series_name="Weighting Day-Minute",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Weighting Scheme")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
fig, ax = plt.subplots(1, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=ax,
    time_series=daily_ghi_df["daily_ghi"],
    time_series_name="Daily GHI",
    s=1,
    c="k"
)
ax.plot(
    daily_ghi_df["estimated_daily_ghi"].index,
    daily_ghi_df["estimated_daily_ghi"].values,
    label="Estimated Daily GHI",
    c="r",
    linestyle="--"
)
plt.legend()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = (
    outsample_forecast["ghi"] - outsample_forecast["estimated_minutely_ghi"]
)

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
y_forecast = np.array(
    np.split(
        outsample_forecast["estimated_minutely_ghi"].values,      
        int(len(outsample_forecast["ghi"].values) / 30)
    )
)
y_true = np.array(
    np.split(
        outsample_forecast["ghi"].values,      
        int(len(outsample_forecast["ghi"].values) / 30)
    )
)

error_mean, error_sem = compute_rmse(
    y_true=y_true,
    y_forecast=y_forecast,
    axis=0
)
MODELS_RESULTS["Top Down GP"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Top Down GP",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

from torch.optim import Adam
from torch.nn import MSELoss

from modules.models.torch_backend.neural_networks.architectures import LinearModel, MLPModel
from modules.models.torch_backend.fitters import TorchModelFitter
from modules.models.torch_backend.callbacks import EarlyStopping


model = LinearModel(
    in_features=AR_WINDOW,
    out_features=OUTSAMPLE_WINDOW_SIZE,
    link_function="ReLU",
).to("cuda")
stopper = EarlyStopping(tolerance=100, min_delta=0.001)
fitter = TorchModelFitter()


model, optimizer, hisotry = fitter.fit(
    model=model,
    optimizer=Adam(model.parameters(), weight_decay=0.01),
    train_loader=train_loader,
    validation_loader=validation_loader,
    stopper=stopper,
    epochs=300,
    loss_function=MSELoss()
)

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = (
    pd.Series(np.vstack(ground_truth).flatten() - np.vstack(forecasts).flatten())
)
residual.index = outsample_ghi_df.index[:len(residual)]

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
error_mean, error_sem = compute_rmse(
    y_true=np.vstack(ground_truth),
    y_forecast=np.vstack(forecasts),
    axis=0
)
MODELS_RESULTS["AR"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Top Down GP",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

from modules.models.torch_backend.neural_networks.architectures import MLPModel

model = MLPModel(
    in_features=AR_WINDOW,
    out_features=OUTSAMPLE_WINDOW_SIZE,
    layers=(30, 30, 30),
    link_function="ReLU",
    activation_function="ReLU",
    dropout_rate=0.1,
).to("cuda")
stopper = EarlyStopping(tolerance=100, min_delta=0.001)
fitter = TorchModelFitter()


model, optimizer, hisotry = fitter.fit(
    model=model,
    optimizer=Adam(model.parameters(), weight_decay=0.01),
    train_loader=train_loader,
    validation_loader=validation_loader,
    stopper=stopper,
    epochs=300,
    loss_function=MSELoss()
)

#| code-fold: true
#| code-summary: "Show supplementary code"
residual = (
    pd.Series(np.vstack(ground_truth).flatten() - np.vstack(forecasts).flatten())
)
residual.index = outsample_ghi_df.index[:len(residual)]

fig, axs = plt.subplots(2, 1, figsize=(8, 6))

ax = visualize_time_series(
    ax=axs[0],
    time_series=residual,
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)
ax=visualize_time_series(
    axs[1],
    time_series=residual[-(60 * 24 * 30):],
    time_series_name="Residuals",
    s=0.05,
    alpha=0.5
)


for ax in axs:
    ax.axhline(
        y=0,
        linestyle=":",
        c="k"
    )

plt.suptitle("Mean Persistence Residuals")
plt.tight_layout()
plt.show()

#| code-fold: true
#| code-summary: "Show supplementary code"
error_mean, error_sem = compute_rmse(
    y_true=np.vstack(ground_truth),
    y_forecast=np.vstack(forecasts),
    axis=0
)
MODELS_RESULTS["MLP AR"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="Top Down GP",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

import numpy as np
from jax import random

from sklearn.preprocessing import MinMaxScaler

from modules.data_handling.features_generators import make_fourier_matrix


master_key = random.PRNGKey(666)

ghi_scaler = MinMaxScaler()
satellite_covariates_scaler = MinMaxScaler()

y_insample = insample_ghi_df["ghi"].values
y_outsample = outsample_ghi_df["ghi"].values

insample_satellite_covariates = insample_ghi_df[np.arange(1, 101)].values
insample_satellite_covariates = satellite_covariates_scaler.fit_transform(insample_satellite_covariates)

outsample_satellite_covariates = outsample_ghi_df[np.arange(1, 101)].values
outsample_satellite_covariates = satellite_covariates_scaler.fit_transform(outsample_satellite_covariates)

insample_time = np.arange(len(y_insample))
outsample_time = np.arange(insample_time.max(), insample_time.max() + len(y_outsample))

insample_fourier_matrix_year = make_fourier_matrix(
    time_series_indices=np.arange(len(y_insample)), 
    n_components=10, 
    period=525949.2
)
insample_fourier_matrix_day = make_fourier_matrix(
    time_series_indices=np.arange(len(y_insample)), 
    n_components=10, 
    period=1440
)
outsample_fourier_matrix_year = make_fourier_matrix(
    time_series_indices=outsample_time, 
    n_components=10, 
    period=525949.2
)
outsample_fourier_matrix_day = make_fourier_matrix(
    time_series_indices=outsample_time, 
    n_components=10, 
    period=1440
)

from modules.models.jax_backend.fitters import fit
from modules.models.jax_backend.optimizers import newton_rhapson
from modules.models.jax_backend.callbacks import EarlyStopping

params, history = fit(
    X=(
        insample_fourier_matrix_year,
        insample_fourier_matrix_day,
        insample_satellite_covariates
    ),
    y=y_insample, 
    backward=backward, 
    start_params=start_params, 
    optimizer=newton_rhapson(learning_rate=10),
    epochs=1000, 
    stopper=EarlyStopping(tol=100), 
    verbose=1, 
    batch_size=None
)

#| code-fold: true
#| code-summary: "Show supplementary code"
y_forecast = np.array(
    np.split(
        relu(sum(components_outsample.values())),      
        int(len(y_outsample) / 30)
    )
)
y_true = np.array(
    np.split(
        y_outsample,      
        int(len(y_outsample) / 30)
    )
)

error_mean, error_sem = compute_rmse(
    y_true=y_true,
    y_forecast=y_forecast,
    axis=0
)
MODELS_RESULTS["GAM - MLP"] = {"error_mean": error_mean, "error_sem": error_sem}

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax = visualize_forecast_horizon_accuracy(
    ax=ax,
    model_name="GAM MLP",
    error_name="RMSE",
    error_mean=error_mean,
    error_sem=error_sem,
    marker="o",
    markersize=5
)
plt.legend()
plt.show()

fig = visualize_time_series_components_performance(
    time_series_train=insample_ghi_df["ghi"],
    time_series_test=outsample_ghi_df["ghi"],
    components_insample=components_insample,
    components_outsample=components_outsample,
    link_function=relu,
    target_name="GHI",
    model_name="Custom Seasonal Regressor",
    guardrail_metric=compute_symmetric_absolute_percentage_error,
    loss="MAE",
    history=history,
    figsize=(15, 5),
)
plt.show()

truncated_components_insample = {
    component_name: component_value[-3000:] for component_name, component_value in components_insample.items()
}

truncated_components_outsample = {
    component_name: component_value[:3000] for component_name, component_value in components_outsample.items()
}


fig = visualize_time_series_components_performance(
    time_series_train=insample_ghi_df["ghi"][-3000:],
    time_series_test=outsample_ghi_df["ghi"][:3000],
    components_insample=truncated_components_insample,
    components_outsample=truncated_components_outsample,
    link_function=relu,
    target_name="GHI",
    model_name="Custom Seasonal Regressor",
    guardrail_metric=compute_symmetric_absolute_percentage_error,
    loss="MAE",
    history=history,
    figsize=(15, 5),
)
plt.show()

fig, ax = plt.subplots(1, 1, figsize=(20, 10))

for model_name, result in MODELS_RESULTS.items():

    ax = visualize_forecast_horizon_accuracy(
        ax=ax,
        model_name=model_name,
        error_name="RMSE",
        error_mean=result["error_mean"],
        error_sem=result["error_sem"],
        marker="o",
        markersize=5
    )

plt.legend()
plt.show()