#%%
import pandas as pd
import Models as models
import pre_processing as pp
import seaborn as sns
import Main_functions_DP as MfD
from importlib import reload
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.font_manager as fm

import matplotlib.gridspec as gridspec
import matplotlib as mpl
from sklearn.metrics import mean_absolute_error, r2_score, mean_squared_error
import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import locale
import seaborn as sns
import matplotlib.font_manager as font_manager
import pandas as pd
from collections import defaultdict
import statsmodels.api as sm
from sklearn.preprocessing import MinMaxScaler
from scipy.stats import gaussian_kde
from matplotlib import cm
# reload the models module
MfD = reload(MfD)
season_colors = {'Winter': '#8bb3e6', 'Spring': '#90d17f', 'Summer': '#ff9a87', 'Fall': '#ffa04c'}
model_colors_pred = {'ANN': '#56d660', 'LGBM': '#6056d6', 'LR': '#d66056'}
model_colors_real = {'ANN': '#060e00', 'LGBM': '#060e00', 'LR': '#060e00'}
binary_colours = {"blue": "#ff0000", "orange": "#9500ff"}
basic_colours = {"blue": "#3070CF", "orange": "#CF8F30"}

pred_W_trnas = pd.DataFrame()
df = pp.general_file_preprocessing(pd.read_excel('data_dp.xlsx'))

df["mst_TG_calc"] = df["total_steam"] - (df["m_st_bo"]+df["de_air"]+df["external_cons"]+df["lin_bp"])
df_filled = MfD.fill_missing_Q_req(df, 'Q_req.xlsx')


# Set the global font to be Times New Roman
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']
# Define the font properties
font = fm.FontProperties(family='Times New Roman', size=16)

def plot_quantile_kde(df, x_col, y_col, x_quantiles=(0.05, 0.95), y_quantiles=(0.05, 0.95), x_label=None, y_label=None, save_name='figure.svg'):

    # Set 'Times New Roman' as the default font.
    mpl.rcParams['font.family'] = 'Times New Roman'
    mpl.rcParams['font.size'] = 16
    # Calculate quantiles for x_col and y_col
    x_low, x_high = df[x_col].quantile(x_quantiles)
    y_low, y_high = df[y_col].quantile(y_quantiles)
    # Filter data based on quantiles
    filtered_df = df[(df[x_col].between(x_low, x_high)) & (df[y_col].between(y_low, y_high))]
    # Create jointplot with KDE and 'coolwarm' colormap
    g = sns.jointplot(data=filtered_df, x=x_col, y=y_col, kind='kde', fill=False, cmap='coolwarm')
    # Set equal axis scales for comparability
    lim = max(x_high, y_high)
    g.ax_joint.set_xlim(0, lim)
    g.ax_joint.set_ylim(0, lim)
    # Add dashed diagonal line
    g.ax_joint.plot([0, lim], [0, lim], color='blue', linestyle='--', alpha=0.8)
    g.ax_joint.plot([-25, lim-25], [0, lim], color='pink', linestyle='--', alpha=0.5)
    # Set the font properties
    for label in (g.ax_joint.get_xticklabels() + g.ax_joint.get_yticklabels()):
        label.set_fontname('Times New Roman')
        label.set_fontsize(16)
    # Set axis labels using \mathrm{} for TeX text

    
    
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.savefig(save_name, format='svg', bbox_inches='tight')
    plt.show()

def remove_outliers(df, col_a, col_b, tolerance):
    # Make a copy of the dataframe
    df_copy = df.copy()
    # Calculate the absolute difference between col_a and col_b
    diff = abs(df_copy[col_a] - df_copy[col_b])
    # Create a mask where the difference is less than or equal to the tolerance
    mask = diff <= tolerance
    # Number of points to be removed
    num_removed = len(df_copy) - mask.sum()
    # Remove rows where the difference is greater than the tolerance
    df_clean = df_copy.loc[mask]
    # Print the count of removed points
    print(f'Removed data: {num_removed} ({num_removed/len(df)*100:.2f}%)')
    return df_clean
# Use the function
df_ready = df_filled


# set up the arguments
classes_names = ["m_st_total","m_st_ec","ms_st_bo","m_st_da","m_st_sc_ex"]
main_columns = ["total_steam","external_cons","m_st_bo","de_air","mst_sf_demi"]
general_factors = ["boilers_active","Fall","Spring","Summer","Winter","month","day","hour",
                 "t_avg","h_avg","rain_avg","T_amb"]
general_factors_TDW = general_factors.copy()
zscore_threshold = 3
classes = []
# call the function
for col in main_columns:
    pp.plot_column_distribution(df_ready, col, bins=80, zscore_threshold=2, lower_quantile=0.10, upper_quantile=0.90)

classes, df_ready = MfD.create_steam_stream_classes(df_ready, classes_names, main_columns, general_factors, zscore_threshold)

df_ready, classes, general_factors, general_factors_TDW, main_columns = MfD.preprocess_W_gen_and_W_trans(df_ready, classes, general_factors, general_factors_TDW, main_columns)


def plot_removed_values(original_df, cleaned_df, season_colors):
    # Find the removed rows by index
    removed_df = original_df[~original_df.index.isin(cleaned_df.index)]

    # Count removed values by year and season
    counts = defaultdict(lambda: defaultdict(int))
    for _, row in removed_df.iterrows():
        year = row['year']
        for season in ['Winter', 'Spring', 'Summer', 'Fall']:
            if row[season] == 1:
                counts[year][season] += 1

    # Convert counts to dataframe
    count_df = pd.DataFrame(counts).T.fillna(0)

    # Sort by year
    count_df.sort_index(inplace=True)

    # Create the plot
    plt.figure(figsize=(10, 6))
    bottom = [0] * len(count_df)
    for season in ['Winter', 'Spring', 'Summer', 'Fall']:
        plt.bar(count_df.index, count_df[season], bottom=bottom, color=season_colors[season], label=season)
        bottom += count_df[season]

    plt.xlabel('Year', fontname='Times New Roman', fontsize=24)
    plt.ylabel('Removed Values Count ', fontname='Times New Roman', fontsize=24)
    plt.title('Number of Removed Values by Year and Season', fontname='Times New Roman', fontsize=24)
    plt.legend(prop={'family': 'Times New Roman', 'size': 12 * 1.5})

    # Remove gridlines
    plt.grid(False)

    # Set the font for axis labels
    font_path = fm.findfont(fm.FontProperties(family='Times New Roman'))
    font_prop = fm.FontProperties(fname=font_path, size=14)
    plt.gca().set_xticklabels(plt.gca().get_xticks(), fontproperties=font_prop)
    plt.gca().set_yticklabels(plt.gca().get_yticks(), fontproperties=font_prop)

    plt.gca().yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
    plt.gca().xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
    plt.savefig("removal distirbution.svg", format='svg', bbox_inches='tight')
    
    plt.show()


season_colors = {'Winter': '#8bb3e6', 'Spring': '#90d17f', 'Summer': '#ff9a87', 'Fall': '#ffa04c'}

plot_removed_values(df_filled[:"2016"], df_ready[:"2016"], season_colors)
df_ready = remove_outliers(df_ready, 'mst_TG', 'mst_TG_calc', 25)
plot_removed_values(df, df_ready, season_colors)
MfD = reload(MfD)
df_ready = df_ready.drop(columns=["year"])
MfD.model_creation(df_ready, classes)

df_ready["Q_req_og"] = df_ready["Q_req"]
df_ready["Q_req"] = df_ready["Q_req_og"]  * 0.75
df_test = df_ready["2016"][:].copy()




#%%

sns.set(style="whitegrid", font="Times New Roman", rc={"grid.linewidth": 0, "axes.grid": False})
def calculate_cascade(df, main_columns=["total_steam", "external_cons", "de_air", "m_st_bo", "mst_sf_demi"],
                      general_factors=["boilers_active", "Fall", "Spring", "Summer", "Winter", "month", "day", "hour",
                                       "t_avg", "h_avg", "rain_avg", "T_amb"], window_size=24, days=2):
    # Combine main columns and general factors
    all_columns = main_columns + general_factors
    # Slice the dataframe to include only required columns
    df = df[all_columns]
    # Scale the dataframe values between 0 and 1
    scaler = MinMaxScaler()
    df_scaled = pd.DataFrame(scaler.fit_transform(df), columns=df.columns, index=df.index)  # Scales all the columns between 0 and 1
    # Apply moving average
    df_ma = df_scaled.rolling(window_size).mean().dropna()  # Applies a moving average of window_size on the data
    # Calculate the correlation matrix for the main columns
    correlation_matrix = df_ma.corr(method='pearson')  # Calculate the correlation matrix
    # Print correlation matrix
    print(correlation_matrix.loc[main_columns, main_columns])
    # Calculate the sum of correlations for each column
    correlation_sum = correlation_matrix.sum()
    # Sort by the sum of correlations
    colors = ['#568CD6',"#ffffff","#f1b460", '#ff4040']  # Red and Blue
    # Create a custom colormap using the defined colors
    print(correlation_sum[main_columns])
    cascade_order = correlation_sum[main_columns].sort_values().index.tolist()
    # Plot correlation matrix
    fig, ax = plt.subplots(figsize=(10, 5))
    sns.heatmap(correlation_matrix.loc[main_columns, main_columns], annot=True, fmt=".2f", ax=ax, cmap="coolwarm",
                cbar_kws={'label': 'Correlation coefficient'})  # Plotting the correlation matrix
    plt.title('Correlation Matrix and Order of Prediction', fontsize=18)
    ax.set_xticklabels([r'$\mathregular{\dot{m}_{st,total}}$', r'$\mathregular{\dot{m}_{external}}$', r'$\mathregular{\dot{m}_{de-air}}$', 
                        r'$\mathregular{\dot{m}_{st,blow-off}}$', r'$\mathregular{\dot{m}_{st,demi}}$'], fontsize=18) 
    ax.set_yticklabels([r'$\mathregular{\dot{m}_{st,total}}$', r'$\mathregular{\dot{m}_{external}}$', r'$\mathregular{\dot{m}_{de-air}}$', 
                        r'$\mathregular{\dot{m}_{st,blow-off}}$', r'$\mathregular{\dot{m}_{st,demi}}$'], fontsize=18)  
    plt.savefig("TEST4.svg", format='svg', bbox_inches='tight')
    plt.show()
    # Print order of prediction and their sums
    print("\nOrder of prediction and their sums:")
    for col in cascade_order:
        print(f"{col}: {correlation_sum[col]:.2f}")

    # Plot X days slice of data
    fig, ax = plt.subplots(figsize=(8, 5))
    df_ma[main_columns].iloc[24 * (days + 1500):24 * (days + 1510)].plot(ax=ax)  # Plotting X days slice of data
    plt.title('Last X days slice of data', fontsize=18)
    plt.ylabel('Scaled steam flow rate')
    plt.legend(loc='lower right', fontsize=14)

    plt.show()
    return cascade_order

cascade_order = calculate_cascade(df_ready)
print(cascade_order)


#%%
def predict_general(model,scaler, real_data_df):
    cols_atm_aviable = scaler.get_feature_names_out().tolist()
    X = real_data_df[cols_atm_aviable].reset_index(drop=True)
    X = scaler.transform(X)
    # Make predictions
    y_pred = model.predict(X)
    return y_pred

df = df_ready
df_test = df_ready["2016"][:].copy()

#df_test = df_ready[df_ready["year"]>2015].copy()
#compute_predictions ------- (df, classes, main_columns,general_factors, start, end):
cls_steam = classes[0]

real_df = df_test.copy()
real_data_X, real_data_Y, pred_df = MfD.create_lagged_data(real_df, main_columns,general_factors, start=24, end=48)
real_data_X, pred_df = MfD.add_pred_to_df(cls_steam, real_data_X, pred_df)

model_used = ["_LGBM","_LGBM","_LR","_LR","_LGBM"]
real_data_X, real_data_Y, pred_df = MfD.predict_classes(real_data_X, real_data_Y, pred_df, classes[:-2],model_used)

# Iterate over each column pair and compute the error and residual
cls = classes[-2]


w_gen_cols = ["total_steam","external_cons","de_air","m_st_bo",'t_avg',"h_avg"]
real_df_W_max = df_ready[w_gen_cols].loc[real_data_X.index].copy().dropna()
Q_req = df_ready['Q_req'].loc[real_df_W_max.index].values
real_df_W_max["mst_bp"] = 0

real_df_W_max['mst_TG_calc'] = MfD.calculate_mst_TG(real_df_W_max)
#hell
import numpy as np
mst_TG =  real_df_W_max['mst_TG_calc'].values
mst_ex_sc = df['mst_sf_demi'].loc[real_df_W_max.index].values
mst_ex_cal = np.zeros_like(Q_req)
mst_bp = np.zeros_like(Q_req)
T_min = df['T_min'].loc[pred_df.index].values
for i in range(len(Q_req)):
    progress_pct = (i+1)/len(Q_req)*100
    print(f"\rCalculating MST Ex original... {progress_pct:.2f}% complete", end="")
    mst_ex_cal[i],mst_bp[i] = pp.calculate_mst_ex(Q_req[i], mst_TG[i],  mst_ex_sc[i],T_min =T_min[i])
    print("\n")


real_df_W_max['mst_ex']  = mst_ex_cal
real_df_W_max['mst_bp'] = mst_bp
real_df_W_max['mst_TG_calc'] = MfD.calculate_mst_TG(real_df_W_max)
model_name = "_ANN"
cols_atm_aviable = cls.scaler_dict[cls.column_name + model_name].get_feature_names_out().tolist()


model_dict = classes[-2].model_dict
scaler_dict = classes[-2].scaler_dict

A_model = model_dict['W_gen_LGBM']
A_scaler = scaler_dict['W_gen_LGBM']
A_pred = predict_general(A_model,A_scaler, real_df_W_max)
B_model = model_dict['W_gen_ANN']
B_scaler = scaler_dict['W_gen_ANN']
B_pred = predict_general(B_model,B_scaler, real_df_W_max)
C_model = model_dict['W_gen_LR']
C_scaler = scaler_dict['W_gen_LR']
C_pred = predict_general(C_model,C_scaler, real_df_W_max)

real_df_W_max["W_gen_0_LGBM"] = A_pred
real_df_W_max["W_gen_0_ANN"] = B_pred
real_df_W_max["W_gen_0_LR"] = C_pred


df_real_vals = df[w_gen_cols].loc[real_data_X.index].copy().dropna()
Q_req = df_ready['Q_req'].loc[real_df_W_max.index].values
df_real_vals["mst_bp"] = df["lin_bp"].loc[real_data_X.index].copy().dropna()

df_real_vals['mst_TG_calc'] = MfD.calculate_mst_TG(df_real_vals)
#hell
import numpy as np
mst_TG =  df_real_vals['mst_TG_calc'].values
mst_ex_sc = df['mst_sf_demi'].loc[df_real_vals.index].values
mst_ex_cal = np.zeros_like(Q_req)
mst_bp = np.zeros_like(Q_req)
T_min = df['T_min'].loc[pred_df.index].values
for i in range(len(Q_req)):
    progress_pct = (i+1)/len(Q_req)*100
    print(f"\rCalculating MST Ex original... {progress_pct:.2f}% complete", end="")
    mst_ex_cal[i],mst_bp[i] = pp.calculate_mst_ex(Q_req[i], mst_TG[i],  mst_ex_sc[i],T_min =T_min[i])
    print("\n")


df_real_vals['mst_ex']  = mst_ex_cal
df_real_vals['mst_bp'] = df_real_vals["mst_bp"] + mst_bp
df_real_vals['mst_TG_calc'] = MfD.calculate_mst_TG(df_real_vals)

df_real_vals["W_gen_bp_LGBM"] = predict_general(A_model,A_scaler, df_real_vals)
df_real_vals["W_gen_bp_ANN"] = predict_general(B_model,B_scaler, df_real_vals)
df_real_vals["W_gen_bp_LR"] = predict_general(C_model,C_scaler, df_real_vals)





pred_df['mst_TG_calc'] = pred_df['total_steam_pred'] \
                                    - pred_df['external_cons_pred'] \
                                    - pred_df['de_air_pred'] \
                                    - pred_df['m_st_bo_pred']
pred_df['mst_TG'] = real_df_W_max['mst_TG_calc']

#hell
Q_req = df_ready['Q_req'].loc[real_data_X.index].values
mst_TG =  pred_df['mst_TG_calc'].values
mst_ex_sc = pred_df['mst_sf_demi_pred'].values
mst_ex_cal = np.zeros_like(Q_req)
mst_bp = np.zeros_like(Q_req)
for i in range(len(Q_req)):
    progress_pct = (i+1)/len(Q_req)*100
    print(f"\rCalculating MST Ex prediction... {progress_pct:.2f}% complete", end="")
    mst_ex_cal[i],mst_bp[i] = pp.calculate_mst_ex(Q_req[i], mst_TG[i],  mst_ex_sc[i],T_min =T_min[i])
    print("\n")

pred_df['t_avg'] = df["t_avg"].loc[real_data_X.index]
pred_df['mst_ex'] = mst_ex_cal
pred_df["mst_bp_pred"]= mst_bp
pred_df['h_avg'] = df["h_avg"].loc[real_data_X.index]

pred_df['mst_TG_calc'] = pred_df['total_steam_pred'] \
                                - pred_df['external_cons_pred'] \
                                - pred_df['de_air_pred'] \
                                - pred_df['m_st_bo_pred'] - pred_df["mst_bp_pred"]


pred_df["W_gen_0_LGBM_pred"] = predict_general(A_model,A_scaler, pred_df)
pred_df["W_gen_0_ANN_pred"] = predict_general(B_model,B_scaler, pred_df)
pred_df["W_gen_0_LR_pred"] = predict_general(C_model,C_scaler, pred_df)




pred_df["W_gen_bp_LGBM"] = df_real_vals["W_gen_bp_LGBM"] 
pred_df["W_gen_bp_ANN"] = df_real_vals["W_gen_bp_ANN"]
pred_df["W_gen_bp_LR"] = df_real_vals["W_gen_bp_LR"] 

pred_df["W_gen_0_LGBM"] = real_df_W_max["W_gen_0_LGBM"] 
pred_df["W_gen_0_ANN"] = real_df_W_max["W_gen_0_ANN"] 
pred_df["W_gen_0_LR"] = real_df_W_max["W_gen_0_LR"] 






pred_df = pred_df.rename(columns={'mst_TG_calc': 'mst_TG_pred'})

pred_df['mst_TG_calc'] = pred_df['total_steam_pred_max'] \
                                    - pred_df['external_cons_pred'] \
                                    - pred_df['de_air_pred'] \
                                    - pred_df['m_st_bo_pred'] - pred_df["mst_bp_pred"]

pred_df["W_gen_MAX_LGBM"] = predict_general(A_model,A_scaler, pred_df)
pred_df["W_gen_MAX_ANN"] = predict_general(B_model,B_scaler, pred_df)
pred_df["W_gen_MAX_LR"] = predict_general(C_model,C_scaler, pred_df)


pred_df = pred_df.rename(columns={'mst_TG_calc': 'mst_TG_calc_max'})

pred_df['mst_TG_calc'] = pred_df['total_steam_pred_min'] \
                                    - pred_df['external_cons_pred'] \
                                    - pred_df['de_air_pred'] \
                                    - pred_df['m_st_bo_pred'] - pred_df["mst_bp_pred"]


pred_df["W_gen_min_LGBM"] = predict_general(A_model,A_scaler, pred_df)
pred_df["W_gen_min_ANN"] = predict_general(B_model,B_scaler, pred_df)
pred_df["W_gen_min_LR"] = predict_general(C_model,C_scaler, pred_df)

pred_df = pred_df.rename(columns={'mst_TG_calc': 'mst_TG_calc_min'})


pred_df['mst_ex_pred'] = pred_df['mst_ex']  
pred_df["mst_ex"] = real_df_W_max['mst_ex'] 
pred_df["mst_bp"] = real_df_W_max['mst_bp'] 

pred_df['mst_ex_real'] = df_real_vals['mst_ex']  
pred_df["mst_bp_real"] = df_real_vals['mst_bp'] 


pred_df["W_plan"] = 0
# Define a lambda function to apply the logic to each row
def calculate_total_steam_plan(row,cls):
    if abs(row["W_gen_0_ANN_pred"] - row["W_gen_min_ANN"]) > 2:
        return row["W_gen_min_ANN"] + (1/4) * (row["W_gen_0_ANN_pred"] - row["W_gen_min_ANN"])
    else:
        return row["W_gen_min_ANN"] + (1/2) * (row["W_gen_0_ANN_pred"] - row["W_gen_min_ANN"])
    
pred_df['W_plan'] = pred_df.apply(lambda row: calculate_total_steam_plan(row,cls), axis=1)
pred_df[["W_gen_bp_ANN","W_gen_min_ANN","W_plan"]].sum()


pred_df['W_plan'] = pred_df["W_gen_min_ANN"]

pred_df_safe = pred_df.copy() 


from sklearn.metrics import mean_absolute_error

cls = classes[-1]
cols_atm_aviable = cls.scaler_dict[cls.column_name + "_LR"].get_feature_names_out().tolist()

real_data_X["W_gen"] = pred_df["W_plan"]
real_data_X[cols_atm_aviable]
model_trans_A = cls.model_dict["W_trans_LR"]
scaler_trans_A = cls.scaler_dict["W_trans_LR"]
model_trans_B = cls.model_dict["W_trans_ANN"]
scaler_trans_B = cls.scaler_dict["W_trans_ANN"]

model_trans_C = cls.model_dict["W_trans_LGBM"]
scaler_trans_C = cls.scaler_dict["W_trans_LGBM"]



cur_list = scaler_trans_A.get_feature_names_out().tolist()
model_name = "W_trans_LR"
target_col = cls.column_name

if scaler_trans_A:
    X_test_scaled = scaler_trans_A.transform(real_data_X[cur_list])
    if 'LR' in model_name and target_col != "W_gen":
        X_test_scaled_df = pd.DataFrame(X_test_scaled, columns=cur_list)
        X_test_scaled_df["const"] = 1
        X_test_scaled = X_test_scaled_df.values
    y_pred = model_trans_A.predict(X_test_scaled)

real_data_X = real_data_X.rename(columns={'W_gen': 'W_plan'})

pred_df["W_trans_pred"] = y_pred


cur_list = scaler_trans_A.get_feature_names_out().tolist()
model_name = "W_trans_LR"
real_data_X["W_gen"] = pred_df["W_gen_0_ANN"]
if scaler_trans_A:
    X_test_scaled = scaler_trans_A.transform(real_data_X[cur_list])
    if 'LR' in model_name and target_col != "W_gen":
        X_test_scaled_df = pd.DataFrame(X_test_scaled, columns=cur_list)
        X_test_scaled_df["const"] = 1
        X_test_scaled = X_test_scaled_df.values
    y_pred = model_trans_A.predict(X_test_scaled)
    pred_df["W_trans_0_LR"] = y_pred

    cur_list = scaler_trans_B.get_feature_names_out().tolist()
    X_test_scaled = scaler_trans_B.transform(real_data_X[cur_list])
    y_pred = model_trans_B.predict(X_test_scaled)
    pred_df["W_trans_0_ANN"] = y_pred

    cur_list = scaler_trans_C.get_feature_names_out().tolist()
    X_test_scaled = scaler_trans_C.transform(real_data_X[cur_list])
    y_pred = model_trans_C.predict(X_test_scaled)
    pred_df["W_trans_0_LGBM"] = y_pred


real_data_X = real_data_X.rename(columns={'W_gen': 'W_gen_0_Ann'})

preds_df = pd.DataFrame()
# Calculate the differences between ANN, LGBM, and LR predictions
preds_df["W_trans_0_ANN"] = pred_df["W_trans_0_ANN"]
preds_df["W_trans_0_LGBM"] = pred_df["W_trans_0_LGBM"]
preds_df["W_trans_0_LR"] = pred_df["W_trans_0_LR"]


preds_df['ANN_diff_LR'] = pred_df['W_trans_0_ANN'] - pred_df['W_trans_0_LR']
preds_df['LGBM_diff_LR'] = pred_df['W_trans_0_LGBM'] - pred_df['W_trans_0_LR']

# Calculate the quantiles
quantiles = [0.33, 0.66, 0.99]
quantile_names = ['Q33', 'Q66', 'Q99']

results = []

for q, q_name in zip(quantiles, quantile_names):
    quantile_value = preds_df['W_trans_0_LR'].quantile(q)
    mask = preds_df['W_trans_0_LR'] >= quantile_value

    mae_ANN = mean_absolute_error(preds_df.loc[mask, 'W_trans_0_LR'], preds_df.loc[mask, 'W_trans_0_ANN'])
    mae_LGBM = mean_absolute_error(preds_df.loc[mask, 'W_trans_0_LR'], preds_df.loc[mask, 'W_trans_0_LGBM'])

    results.append({'quantile': q_name, 'model': 'W_trans_0_ANN', 'mae': mae_ANN})
    results.append({'quantile': q_name, 'model': 'W_trans_0_LGBM', 'mae': mae_LGBM})

mae_quantiles_df = pd.DataFrame(results)
print(mae_quantiles_df)



sns.set(style="whitegrid")
plt.figure(figsize=(10, 6))

# Create a bar plot
mae_quantiles_plot = sns.barplot(x='quantile', y='mae', hue='model', data=mae_quantiles_df)

# Set plot labels and title
plt.xlabel('Quantile', fontsize=14)
plt.ylabel('Mean Absolute Error (MAE)', fontsize=14)
plt.title('Comparison of ANN and LGBM Model Performance at Different Quantiles of LR Model', fontsize=16)
# Show the plot
plt.show()


def plot_predictions(df, ann_col, lr_col, lgbm_col, title):
    fig, ax = plt.subplots(figsize=(8, 4))
    ann_color = model_colors_pred["ANN"]
    lr_color = model_colors_pred["LR"]
    lgbm_color = model_colors_pred["LGBM"]
    # Plot the predictions for each model
    sns.lineplot(data=df, x=df.index, y=df[ann_col], label='ANN', color=ann_color, alpha=0.7)
    sns.lineplot(data=df, x=df.index, y=df[lr_col], label='LR', color=lr_color, linestyle='--', linewidth=2)
    sns.lineplot(data=df, x=df.index, y=df[lgbm_col], label='LGBM', color=lgbm_color, alpha=0.7)

    # Set plot labels and title
    ax.set_xlabel('Index', fontsize=14)
    ax.set_ylabel('Predicted Value', fontsize=14)
    ax.set_title(title, fontsize=16)
    #    models.save_figure(fig, title)
    # Add legend and show the plot
    ax.legend(fontsize=12)
    plt.show()



# Example usage:
# Replace 'pred_df' with your DataFrame containing the predictions, and
# adjust the column names accordingly
plot_predictions(pred_df[1000:1100], 'W_trans_0_ANN', 'W_trans_0_LR', 'W_trans_0_LGBM', 'W_trans Predictions Comparison')
def plot_heatmap(x_col, y_col, pred_df):
    # Plot x_col vs. y_col as a heatmap
    fig, ax = plt.subplots(figsize=(8, 6))
    hb = ax.hexbin(pred_df[x_col], pred_df[y_col], gridsize=50, cmap='binary')
    ax.set_xlabel(x_col)
    ax.set_ylabel(y_col)
    ax.set_title(f'{y_col} vs. {x_col}')
    cb = fig.colorbar(hb, ax=ax)
    cb.set_label('count')
#we stick with LR!
plot_heatmap('W_trans_0_LR', 'W_trans_pred', pred_df)


def plot_heatmap(df, x_col, y_col, gridsize=33, cmap='binary'):
    fig, ax = plt.subplots(figsize=(8, 6))
    hb = ax.hexbin(df[x_col], df[y_col], gridsize=gridsize, cmap=cmap, mincnt=1)
    ax.set_xlabel(x_col)
    ax.set_ylabel(y_col)
    ax.set_title(f'{y_col} vs. {x_col}')
    cb = fig.colorbar(hb, ax=ax)
    cb.set_label('count')

    # Add boundary lines
    ax.plot([min(df[x_col]), max(df[x_col])], [min(df[x_col])+0.5, max(df[x_col])+0.5], 'r--')
    ax.plot([min(df[x_col]), max(df[x_col])], [min(df[x_col])-0.5, max(df[x_col])-0.5], 'r--')
    
    # Calculate percentage of points above, within, and below boundaries
    count_above = (df[y_col] - df[x_col]) > 0.5
    count_below = (df[y_col] - df[x_col]) < -0.5
    count_within = (df[y_col] - df[x_col]).between(-0.5, 0.5)
    
    percent_above = count_above.sum() / len(df) * 100
    percent_below = count_below.sum() / len(df) * 100
    percent_within = count_within.sum() / len(df) * 100
    
    ax.set_title(f'{y_col} vs. {x_col}\n{percent_above:.1f}% above, {percent_within:.1f}% within, {percent_below:.1f}% below')

    plt.show()

plot_heatmap(pred_df, 'W_trans_0_LR', 'W_trans_pred')


cur_list = scaler_trans_A.get_feature_names_out().tolist()
model_name = "W_trans_LR"
real_data_X["W_gen"] = pred_df["W_gen_bp_ANN"]

if scaler_trans_A:
    X_test_scaled = scaler_trans_A.transform(real_data_X[cur_list])
    if 'LR' in model_name and target_col != "W_gen":
        X_test_scaled_df = pd.DataFrame(X_test_scaled, columns=cur_list)
        X_test_scaled_df["const"] = 1
        X_test_scaled = X_test_scaled_df.values
    y_pred = model_trans_A.predict(X_test_scaled)
real_data_X = real_data_X.rename(columns={'W_gen': 'W_gen_bp_ANN'})


# Apply seaborn styling
sns.set(style="whitegrid")

# Prepare data
df_final = pd.DataFrame()
df_final["W_trans_0"] = pred_df["W_trans_0_ANN"]
df_final["W_trans_pred"] = pred_df["W_trans_pred"]
#df_final["W_trans_pred_min"] = pred_df["W_trans_pred"]
df_final["W_trans_bp"] = y_pred

# Create 'season' column
month_to_season = {1: 'Winter', 2: 'Winter', 3: 'Spring', 4: 'Spring', 5: 'Spring', 6: 'Summer',
                   7: 'Summer', 8: 'Summer', 9: 'Fall', 10: 'Fall', 11: 'Fall', 12: 'Winter'}
df_final['month'] = df_final.index.month
df_final['season'] = df_final['month'].apply(lambda x: month_to_season[x])

# Calculate success rates
df_final['success'] = (df_final['W_trans_pred'] - df_final['W_trans_0']) <= 0.5
success_rate = (df_final['success'].sum() / len(df_final)) * 100
monthly_success_rate = df_final.groupby('month')['success'].mean() * 100

# Define season colors


# Map month to season colors
monthly_season_colors = df_final[['month', 'season']].drop_duplicates().set_index('month')['season'].map(season_colors)

# Print success rates
print(f"Success rate within ±0.5: {success_rate:.2f}%")
print("\nSuccess rate by month:")
print(monthly_success_rate)

# Bar plot for success rate by month
fig, ax = plt.subplots(figsize=(8, 4))
monthly_success_rate.plot(kind='bar', color=monthly_season_colors, ax=ax)
ax.set_title('Success Rate by Month', fontname='Times New Roman')
ax.set_ylabel('Success Rate (%)', fontname='Times New Roman')
ax.set_xticklabels(['Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'], rotation=0, fontname='Times New Roman')
ax.grid(False)
ax.axhline(95, color='gray', linestyle='--')  # Add line at 95%
ax.yaxis.set_ticks([0, 75, 95, 100])  # Add tick at 95%
ax.axhline(75, color='gray', linestyle='--')  #
# Add season legend with 4 columns
seasons = ['Winter', 'Spring', 'Summer', 'Fall']
handles = [plt.Rectangle((0, 0), 1, 1, color=season_colors[s]) for s in seasons]
plt.xlabel("") 
ax.legend(handles, seasons, loc='lower center', ncol=4, bbox_to_anchor=(0.5, -0.21), prop={'family': 'Times New Roman'})
plt.savefig("TEST3.svg", format='svg', bbox_inches='tight')
plt.show()




plt.show()


#%%

def plot_residual_kde(df, main_col, pred_col, tolerance=0.5, x_label=None, y_label=None, y_limit=None, save_name='figure.svg'):

    # Set 'Times New Roman' as the default font.
    plt.rcParams['font.family'] = 'Times New Roman'
    plt.rcParams['font.size'] = 12

    # Calculate residuals
    df['residual'] = df[main_col] - df[pred_col]

    # Create jointplot with KDE (only lines)
    g = sns.jointplot(data=df, x=main_col, y='residual', kind='kde', fill=False, cmap='coolwarm', height=6, ratio=6)

    # Add zero line and tolerance lines using colors from the coolwarm colormap
    g.ax_joint.axhline(0, color=cm.coolwarm(0.5), linestyle='--', alpha=0.8)
    g.ax_joint.axhline(-tolerance, color=cm.coolwarm(0), linestyle='--', alpha=0.5)
    g.ax_joint.axhline(tolerance, color=cm.coolwarm(1), linestyle='--', alpha=0.5)

    # Calculate the number of instances within the tolerance range
    within_tolerance = df[(df['residual'] >= -tolerance) & (df['residual'] <= tolerance)].shape[0]
    above_tolerance = df[df['residual'] > tolerance].shape[0]
    below_tolerance = df[df['residual'] < -tolerance].shape[0]
    total_count = df.shape[0]

    # Add text about the relative count of instances within, above, and below tolerance
    text_str = f"Within tolerance: {within_tolerance} ({100 * within_tolerance / total_count:.1f}%)\n"
    text_str += f"Above tolerance: {above_tolerance} ({100 * above_tolerance / total_count:.1f}%)\n"
    text_str += f"Below tolerance: {below_tolerance} ({100 * below_tolerance / total_count:.1f}%)"
    props = dict(boxstyle='round', facecolor='lightskyblue', alpha=0.2)
    g.ax_joint.text(0.05, 0.95, text_str, transform=g.ax_joint.transAxes, fontsize=14,
        verticalalignment='top', bbox=props)

    # Set the font propertiesrlů
    for label in (g.ax_joint.get_xticklabels() + g.ax_joint.get_yticklabels()):
        label.set_fontname('Times New Roman')
        label.set_fontsize(12)

    # Set axis labels 'W_{transfered}^{\overline{bp}\approx0}'
    plt.xlabel(x_label if x_label else r'Power transfered  ($\mathregular{W_{transfered}^{\overline{bp}\approx0}}$) [MWh]', fontsize=12, fontname='Times New Roman')
    plt.ylabel(y_label if y_label else 'Residuals [MWh]', fontsize=12, fontname='Times New Roman')

    # Set y-axis limit
    if y_limit is not None:
        plt.ylim(y_limit)
        plt.xlim((0, 12))
    # Remove grid
    g.ax_joint.grid(False)
    g.ax_marg_x.grid(False)  # Remove marginal gridlines
    g.ax_marg_y.grid(False)
    
    # Save and display the figure
    plt.savefig(save_name, format='svg', bbox_inches='tight')
    plt.show()

# Usage:
ylbl = r'Residua $\mathregular{(W_{transfered}^{\overline{bp}\approx0}}$ - $\mathregular{\hat{W}_{transfered})}$ [MWh]'
plot_residual_kde(df_final,  "W_trans_0","W_trans_pred", tolerance=0.5, y_limit=(-2, 4), save_name='figure.svg',y_label=ylbl)



ylbl = r'Residua $\mathregular{(W_{transfered}^{\overline{bp}\approx0}}$ - $\mathregular{W_{transfered}^{\overline{bp}>0})}$ [MWh]'
plot_residual_kde(df_final,  "W_trans_0","W_trans_bp", tolerance=0.5, x_label=None, y_label=ylbl, save_name='figureBP.svg', y_limit=(-1, 6))


ylbl = r'Residua $\mathregular{(W_{transfered}^{\overline{bp}\approx0}}$ - $\mathregular{W_{tranQfered}^{\overline{bp}>0})}$ [MWh]'
plot_residual_kde(df_final,  "W_trans_0","W_trans_bp", tolerance=0.5, x_label=None, y_label=ylbl, save_name='figureBP.svg', y_limit=(-1, 6))




x_lbl = r'Total steam production ($\mathregular{\dot{m}_{st,total}}$) [t/h]'
ylbl = r'Residua $\mathregular{\dot{m}_{st,total}}$ - $\mathregular{\hat{\dot{m}}_{st,total,q=0.5}}$ [t/h]'
plot_residual_kde(pred_df,  "total_steam","total_steam_pred", tolerance=4,save_name='Total_steam_5.svg',x_label=x_lbl, y_label=ylbl)

ylbl = r'Residua $\mathregular{\dot{m}_{st,total}}$ - $\mathregular{\hat{\dot{m}}_{st,total,q=0.05}}$ [t/h]'
plot_residual_kde(pred_df,  "total_steam","total_steam_pred_min", tolerance=4,save_name='Total_steam_min.svg',x_label=x_lbl, y_label=ylbl)





def plot_with_highlights(df, date_ranges, column1, column2, column3, interval=0.5, highlight_color_below="#f1b460", highlight_color_above="#D656CC", season_colors=None):
    if season_colors is None:
        season_colors = {'Winter': '#568CD6', 'Spring': '#76D7C4', 'Summer': '#F1C40F', 'Fall': '#D35400'}
    plt.rcParams['font.family'] = 'Times New Roman'
    plt.rcParams['font.size'] = 12
    fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 14))

    for ax, (start, end) in zip(axes, date_ranges):
        current_season = df[start:end]['season'].iloc[0]

        expression_0 = r'$\mathregular{\hat{W}_{transfered}}$'
        expression_1 = r'$\mathregular{W_{transfered}^{\overline{bp}\approx0}}$'
        expression_2 = r'$\mathregular{W_{transfered}^{\overline{bp}>0}}$'

        highlight_mask_below = (df[start:end][column2] - df[start:end][column1]) < -interval
        highlight_mask_above = (df[start:end][column2] - df[start:end][column1]) > interval
        highlight_mask_within = ~highlight_mask_below & ~highlight_mask_above

        within_tolerance = highlight_mask_within.sum()
        above_tolerance = highlight_mask_above.sum()
        below_tolerance = highlight_mask_below.sum()
        total_count = len(df[start:end])
        
        percent_within = (within_tolerance / total_count) * 100
        percent_below = (below_tolerance / total_count) * 100
        percent_above = (above_tolerance / total_count) * 100

        ax.plot(df[start:end][column1], label=f"Power tranfered ({expression_1})", color='black', linestyle='--')
        ax.fill_between(df[start:end].index, df[start:end][column1] - interval, df[start:end][column1] + interval, color='lightgray', alpha=0.5)
        ax.plot(df[start:end][column2], label=f"Power tranfered ({expression_0}) - prediction", color=season_colors[current_season], linewidth=2)
        ax.plot(df[start:end][column3], label=f"Power tranfered ({expression_2}) - original", color='black', linestyle=':')
        ax.scatter(df[start:end][highlight_mask_below].index, df[start:end][highlight_mask_below][column2], color=highlight_color_below)
        ax.scatter(df[start:end][highlight_mask_above].index, df[start:end][highlight_mask_above][column2], color=highlight_color_above)

        ax.set_ylabel('Power transfer [MWh]')
        ax.set_title(current_season) 
        ax.grid(False)

        text_str = f"Within tolerance: {within_tolerance} ({percent_within:.1f}%)\n"
        text_str += f"Above tolerance: {above_tolerance} ({percent_above:.1f}%)\n"
        text_str += f"Below tolerance: {below_tolerance} ({percent_below:.1f}%)"
        props = dict(boxstyle='round', facecolor=season_colors[current_season], alpha=0.2)
        ax.text(0.99, 0.95, text_str, transform=ax.transAxes, ha='right', va='top', bbox=props)

    # Only show the legend under the last subplot
    handles, labels = ax.get_legend_handles_labels()
    legend = fig.legend(handles, labels, loc='lower center', ncol=2, bbox_to_anchor=(0.5, 0.05))
    plt.setp(legend.get_texts(), fontproperties=font_manager.FontProperties(family='Times New Roman'))

    fig.tight_layout()
    plt.subplots_adjust(bottom=0.15)  # Add some space at the bottom for the legend

    plt.savefig("TEST2.svg", format='svg', bbox_inches='tight')
    plt.show()

# Select a 7-day period in the middle of each season
winter_start = pd.to_datetime('2016-02-10')
spring_start = pd.to_datetime('2016-05-15')
summer_start = pd.to_datetime('2016-08-15')
fall_start = pd.to_datetime('2016-10-15')

# Define the date ranges for each season
date_ranges = [
    (winter_start, winter_start + pd.DateOffset(days=6)),
    (spring_start, spring_start + pd.DateOffset(days=6)),
    (summer_start, summer_start + pd.DateOffset(days=6)),
    (fall_start, fall_start + pd.DateOffset(days=6))
]

# Plot the data for each season
plot_with_highlights(df_final, date_ranges, 'W_trans_0', 'W_trans_pred', 'W_trans_bp', interval=0.5)


#%%
df_cash = pd.read_excel('Prices.xlsx')

# Drop the 'Day' and 'Hour' columns
df_cash = df_cash.drop('Hour', axis=1)

# Set the index of the prices dataframe to the datetime column
df_cash.set_index('DateTime', inplace=True)
# Add a column to df_final for the difference between W_trans_0 - W_trans_bp
df_final['Diff_W_trans_0_W_trans_bp'] = df_final['W_trans_0'] - df_final['W_trans_bp']

# Add a column to df_final for the difference between W_trans_0 - W_trans_pred
df_final['Diff_W_trans_0_W_trans_pred'] = df_final['W_trans_0'] - df_final['W_trans_pred']

# Merge the prices dataframe with df_final using the datetime index
df_final = df_final.merge(df_cash, how='left', left_index=True, right_index=True)

# find the 10th percentile of the 'average_price (Kč/MWh)' column
q10 = df_final['average_price (Kč/MWh)'].quantile(0.02)

# replace values below the 10th percentile with the mean value of the column
mean_price = df_final['average_price (Kč/MWh)'].mean()
df_final.loc[df_final['average_price (Kč/MWh)'] < q10, 'average_price (Kč/MWh)'] = mean_price

df_final['average_price (Kč/MWh)'] = df_final['average_price (Kč/MWh)'].fillna(mean_price)


# Multiply the price by the two difference columns to get the potential gain columns
df_final['Potential_Gain_BP'] = df_final['Diff_W_trans_0_W_trans_bp'] * df_final['average_price (Kč/MWh)']
df_final['Potential_Gain_Pred'] = df_final['Diff_W_trans_0_W_trans_pred'] * df_final['average_price (Kč/MWh)']



def calculate_fines_and_profit(df, w_pred):
    # Set locale to Czech
    locale.setlocale(locale.LC_ALL, 'cs_CZ.utf8')

    # Constants
    C_under = 0.15
    C_over = 3.0
    tolerance = 0.5
    w_trans_pred = df[w_pred]
    # Mask for values within the tolerance range
    interval_mask = (w_trans_pred >= df['W_trans_0'] - tolerance) & (w_trans_pred <= df['W_trans_0'] + tolerance)

    # Masks for underestimation and overestimation
    under_mask = (w_trans_pred - df['W_trans_0']) < -tolerance
    over_mask = (w_trans_pred - df['W_trans_0']) > tolerance

    # Calculate fines for underestimation and overestimation
    df['fines'] = 0
    df.loc[under_mask, 'fines'] = abs(w_trans_pred[under_mask] - df[under_mask]['W_trans_0']) * df[under_mask]['average_price (Kč/MWh)'] * C_under
    df.loc[over_mask, 'fines'] = abs(w_trans_pred[over_mask] - df[over_mask]['W_trans_0']) * df[over_mask]['average_price (Kč/MWh)'] * C_over

    # Calculate total fines
    total_fines = df['fines'].sum()

    # Format the total fines as financial value in CZK with separation of thousands and millions
    formatted_total_fines = locale.format_string('%.0f Kč', total_fines, grouping=True)

    # Calculate profit
    profit = (df[interval_mask]['W_trans_0']*df[interval_mask]['average_price (Kč/MWh)']).sum() - total_fines
    values = [df.loc[under_mask, 'fines'].sum(), df.loc[over_mask, 'fines'].sum(), profit]
    # Format the profit as financial value in CZK with separation of thousands and millions
    formatted_profit = locale.format_string('%.0f Kč', profit, grouping=True)
    return formatted_total_fines, formatted_profit,values

formatted_total_fines, formatted_profit,values_bar1 = calculate_fines_and_profit(df_final, "W_trans_pred")
print(formatted_total_fines, formatted_profit,values_bar1)

formatted_total_fines, formatted_profit,values_bar2 = calculate_fines_and_profit(df_final, "W_trans_bp")
print(formatted_total_fines, formatted_profit,values_bar2)



values1 = values_bar1
values2 = values_bar2
labels = ['Underestimation Fines', 'Overestimation Fines', 'Profit']

fig, axes = plt.subplots(1, 3, figsize=(10, 4))

for i, ax in enumerate(axes):
    ax.bar(0, values1[i], width=0.7, label='List 1', color='#568CD6')
    ax.bar(1, values2[i], width=0.7, label='List 2', color='#f1b460')
    ax.set_xticks([0, 1])
    ax.set_xticklabels(['Novel approach', 'Current approach'])
    ax.set_ylabel('Fined value [mil CZK]')
    ax.set_title(labels[i])
    ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x/1_000_000:.1f}'))
    ax.grid(False)
    # Add bar values
        # Adjust y-axis limits
    y_margin = 0.15 * max(values1[i], values2[i])
    ax.set_ylim(0, max(values1[i], values2[i]) + y_margin)
    ax.text(0, values1[i], f'{values1[i]/1_000_000:,.1f} mil CZK', ha='center', va='bottom', color='black',fontsize=10)
    ax.text(1, values2[i], f'{values2[i]/1_000_000:,.1f} mil CZK', ha='center', va='bottom', color='black',fontsize=10)

ax.set_ylabel('Profit for year 2016 [mil CZK]')
plt.tight_layout()
plt.savefig("TEST1.svg", format='svg', bbox_inches='tight')
plt.show()




def calculate_quarterly_profit(df, w_pred):
    # Set locale to Czech
    locale.setlocale(locale.LC_ALL, 'cs_CZ.utf8')

    # Constants
    C_under = 0.15
    C_over = 3.0
    tolerance = 0.5
    w_trans_pred = df[w_pred]

    # Mask for values within the tolerance range
    interval_mask = (w_trans_pred >= df['W_trans_0'] - tolerance) & (w_trans_pred <= df['W_trans_0'] + tolerance)

    # Masks for underestimation and overestimation
    under_mask = (w_trans_pred - df['W_trans_0']) < -tolerance
    over_mask = (w_trans_pred - df['W_trans_0']) > tolerance

    # Calculate fines for underestimation and overestimation
    df['under_fines'] = 0
    df['over_fines'] = 0
    df.loc[under_mask, 'under_fines'] = abs(w_trans_pred[under_mask] - df.loc[under_mask, 'W_trans_0']) * df.loc[under_mask, 'average_price (Kč/MWh)'] * C_under
    df.loc[over_mask, 'over_fines'] = abs(w_trans_pred[over_mask] - df.loc[over_mask, 'W_trans_0']) * df.loc[over_mask, 'average_price (Kč/MWh)'] * C_over

    # Calculate profit
    df['profit'] = (df.loc[interval_mask, 'W_trans_0']*df.loc[interval_mask, 'average_price (Kč/MWh)']) - (df['under_fines'] + df['over_fines'])

    # Group by quarter and calculate total profit for each quarter
    quarterly_df = df.resample('Q').sum()[['profit']]

    return quarterly_df

def compare_predictions(df, w_pred1, w_pred2):
    # Get the quarterly profit for each prediction
    quarterly_df1 = calculate_quarterly_profit(df, w_pred1)
    quarterly_df2 = calculate_quarterly_profit(df, w_pred2)
    
    # Rename columns to distinguish between the two predictions
    quarterly_df1.columns = [f'{col}_{w_pred1}' for col in quarterly_df1.columns]
    quarterly_df2.columns = [f'{col}_{w_pred2}' for col in quarterly_df2.columns]

    # Join the two DataFrames
    comparison_df = quarterly_df1.join(quarterly_df2)

    return comparison_df



def plot_comparisons(comparison_df, w_pred1, w_pred2):
    # Plot the comparison of profit for the two predictions
    fig, ax = plt.subplots(figsize=(8, 4))
    
    bar_width = 0.4  # Adjust the width of the bars
    bar_gap = 0.1  # Adjust the gap between the bars
    
    # Plot the bars with adjusted positions
    bar_positions = np.arange(len(comparison_df.index))
    bar_positions2 = bar_positions + bar_width + bar_gap
    
    bars1 = ax.bar(bar_positions, comparison_df[f'profit_{w_pred1}'], width=bar_width, color='#568CD6')
    bars2 = ax.bar(bar_positions2, comparison_df[f'profit_{w_pred2}'], width=bar_width, color='#f1b460')
    
    plt.title('Profit Comparison by quarters', fontname='Times New Roman')
    plt.ylabel('Profit [CZK]', fontname='Times New Roman')
    ax.set_xticks(bar_positions + (bar_width + bar_gap) / 2)  # Set X tick positions at the center of the bars
    ax.set_xticklabels(['Q 1', 'Q 2', 'Q 3', 'Q 4'], rotation=0, ha='center')  # Align tick labels at the center
    
    ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x/1000000:.1f}'))  # Format y-axis labels
    ax.grid(False)  # Remove grid lines
    ax.legend((bars1[0], bars2[0]), ['Novel approach', 'Current approach'])
    ax.set_ylabel('Profit for year 2016 [mil. CZK]')
    
    # Add bar values
    for rect in bars1 + bars2:
        height = rect.get_height()
        ax.text(rect.get_x() + rect.get_width() / 2, height, f'{height/1_000_000:,.1f} mil. CZK', ha='center', va='bottom',fontsize=10)

        # Adjust Y-axis limits
    y_min, y_max = ax.get_ylim()
    ax.set_ylim(y_min, y_max * 1.05)  # Increase Y-axis limit by 10%
    plt.tight_layout()
    plt.savefig("sdsd.svg", format='svg', bbox_inches='tight')
    plt.show()



# Usage:
# comparison_df = compare_predictions(df, 'w_pred1', 'w_pred2')
# plot_comparisons(comparison_df, 'w_pred1', 'w_pred2')

comparison_df = compare_predictions(df_final, 'W_trans_pred', 'W_trans_bp')
plot_comparisons(comparison_df, 'W_trans_pred', 'W_trans_bp')

def calculate_ratio(dataframe, column1, column2):
    sum1 = dataframe[column1].sum()
    sum2 = dataframe[column2].sum()
    
    ratio = sum1 / sum2
    
    print(f"The ratio of the sum of {column1} to the sum of {column2} is: {ratio}")
calculate_ratio(df_final, 'W_trans_pred', 'W_trans_bp')
#%%





def financial_analysis(df, interval=0.5, fine_val=3):
    # Calculate masks and create new columns
    df['within_interval_mask'] = ((df['W_trans_pred'] - df['W_trans_0']).abs() <= interval)
    df['below_zero_mask'] = (df['W_trans_pred'] - df['W_trans_0']) < -interval
    df['above_zero_mask'] = (df['W_trans_pred'] - df['W_trans_0']) > interval
    df['total_gains_pred'] = df['W_trans_pred'] * df['average_price (Kč/MWh)']
    df['total_gains_bp'] = df['W_trans_bp'] * df['average_price (Kč/MWh)']
    df['gains_difference'] = df['total_gains_pred'] - df['total_gains_bp']
    
    # Calculate gains_within_interval, fines, and adjusted_income_diff
    df['gains_within_interval'] = df['gains_difference'] * df['within_interval_mask']
    df['fines'] = (df['above_zero_mask'] * (df['W_trans_pred'] - df['W_trans_0']) * df['average_price (Kč/MWh)'] * fine_val)
    df['adjusted_income_diff'] = df['gains_difference'] * df['below_zero_mask'] * 0.85

    # Calculate cash flow
    df['cash_flow'] = df['gains_within_interval'] - df['fines'] + df['adjusted_income_diff']

    # Group data by month
    df['month'] = df.index.month
    monthly_cash_flow = df.groupby('month')[['gains_within_interval', 'fines', 'adjusted_income_diff', 'cash_flow', 'total_gains_pred', 'total_gains_bp']].sum()

    # Calculate cumulative gain
    monthly_cash_flow['cumulative_gain'] = np.cumsum(monthly_cash_flow['gains_within_interval'] + monthly_cash_flow['adjusted_income_diff'] - monthly_cash_flow['fines'] )

    # Add quartile information
    monthly_cash_flow['quartile'] = pd.cut(monthly_cash_flow.index, bins=[0, 3, 6, 9, 12], labels=['Q1', 'Q2', 'Q3', 'Q4'], include_lowest=True)

    # Plot total monthly profit comparison (bar plot)
    # Plot total monthly profit comparison (bar plot)
    fig, ax = plt.subplots(figsize=(10, 6))
    monthly_cash_flow.pivot_table(index='quartile', values=['total_gains_pred', 'total_gains_bp'], aggfunc='sum').plot.bar(ax=ax)
    ax.set_title('Total Monthly Profit Comparison')
    ax.set_xlabel('Quartile')
    ax.set_ylabel('Amount')
    ax.legend(loc='upper left')
    plt.show()
    models.save_figure(fig, 'TTotal Monthly Profit Comparison')
    # Plot cumulative gain (line plot)
    fig, ax = plt.subplots(figsize=(10, 6))
    monthly_cash_flow['cumulative_gain'].plot(ax=ax)
    ax.set_title('Cumulative Gain')
    ax.set_xlabel('Month')
    ax.set_ylabel('Cumulative Gain')
    plt.show()
    models.save_figure(fig, 'Cumulative Gain')
    return monthly_cash_flow

cash_flow_analysis = financial_analysis(df_final, interval=0.5, fine_val=3)



def summarize_financial_analysis(monthly_cash_flow):
    summary = monthly_cash_flow.pivot_table(index='quartile', values=['total_gains_pred', 'total_gains_bp', 'gains_within_interval', 'fines', 'adjusted_income_diff', 'cash_flow'], aggfunc='sum')
    yearly_summary = monthly_cash_flow[['total_gains_pred', 'total_gains_bp', 'gains_within_interval', 'fines', 'adjusted_income_diff', 'cash_flow']].sum().to_frame().transpose()
    yearly_summary.index = ['Yearly']
    
    summary_table = pd.concat([summary, yearly_summary])
    
    # Convert values to millions
    summary_table = summary_table / 1_000_000
    
    return summary_table

summary_table = summarize_financial_analysis(cash_flow_analysis)
display(summary_table)



#%%
import pandas as pd

# Assuming you have a DataFrame called 'df'
df = real_data_X

# Get the column names from the DataFrame
real_data = df.columns.tolist()

# Cycle through the classes and save the models, scalers, and factor columns
ANN_model = {}
LR_model = {}
LGBM_model = {}
class_item = classes[0]
# Select a 7-day period in the middle of each season
winter_start = pd.to_datetime('2016-02-10')
spring_start = pd.to_datetime('2016-05-15')
summer_start = pd.to_datetime('2016-08-15')
fall_start = pd.to_datetime('2016-10-15')

# Define the date ranges for each season
date_ranges = [
    (winter_start, winter_start + pd.DateOffset(days=6)),
    (spring_start, spring_start + pd.DateOffset(days=6)),
    (summer_start, summer_start + pd.DateOffset(days=6)),
    (fall_start, fall_start + pd.DateOffset(days=6))
]

for class_item in classes:
    column_name = class_item.column_name
    #ADD LAG COLUMNS AND ITS DONE....
    # Check if the column name exists in the real_data list
    if column_name in real_data:
        # Get the column data from the DataFrame
        column_data = df[column_name]  # Replace 'column_name' with the appropriate column name in your DataFrame
        
        # Save the models
        ANN_model[column_name] = class_item.model_dict[f'{column_name}_ANN']
        LR_model[column_name] = class_item.model_dict[f'{column_name}_LR']
        LGBM_model[column_name] = class_item.model_dict[f'{column_name}_LGBM']

        # Save the scalers
        ANN_model[column_name + '_scaler'] = class_item.scaler_dict[f'{column_name}_ANN']
        LR_model[column_name + '_scaler'] = class_item.scaler_dict[f'{column_name}_LR']
        LGBM_model[column_name + '_scaler'] = class_item.scaler_dict[f'{column_name}_LGBM']

        # Get the factor columns for each model
        ANN_factors = class_item.scaler_dict[f'{column_name}_ANN'].get_feature_names_out().tolist()
        LR_factors = class_item.scaler_dict[f'{column_name}_LR'].get_feature_names_out().tolist()
        LGBM_factors = class_item.scaler_dict[f'{column_name}_LGBM'].get_feature_names_out().tolist()
    else:
        print(f"Column '{column_name}' not found in the DataFrame.")



    # Iterate over the date ranges and make predictions for each season
    for start, end in date_ranges:
        # Slice the data for the current season
        season_data = real_data_X[(real_data_X.index >= start) & (real_data_X.index <= end)]
        
        # Get the feature columns for the current season
        cols_atm_aviable = ...  # Replace with the appropriate code to obtain the feature columns
        
        # Make predictions for the current season using LR model
        y_pred_LR = predict_general2(LR_model[column_name], LR_model[column_name + '_scaler'], season_data[cols_atm_aviable])
        
        # Make predictions for the current season using ANN model
        y_pred_ANN = predict_general(ANN_model[column_name], ANN_model[column_name + '_scaler'], season_data[cols_atm_aviable])
        
        # Make predictions for the current season using LGBM model
        y_pred_LGBM = predict_general(LGBM_model[column_name], LGBM_model[column_name + '_scaler'], season_data[cols_atm_aviable])
        
        column_data = df[column_name]
        # Print or use the predictions for further processing
        print(f"Predictions for {start.strftime('%Y-%m-%d')} to {end.strftime('%Y-%m-%d')}:")
        print("LR Model:", y_pred_LR)
        print("ANN Model:", y_pred_ANN)
        print("LGBM Model:", y_pred_LGBM)

import statsmodels.api as sm



#%%



pred_df

def calculate_proportion_below_real(df, real_column, prediction_column):
    total_count = len(df)
    count_below = sum(df[prediction_column] < df[real_column])

    proportion = count_below / total_count

    return proportion

def plot_residuals(df, main_column_A, main_column_B, prediction_A, prediction_B, x_label = "test", y_label="test", save_name="test.svg"):
    # Calculate the residuals for prediction A and prediction B
    dataframe = df.copy()
    dataframe['residuals_A'] = dataframe[prediction_A] - dataframe[main_column_A]
    dataframe['residuals_B'] = dataframe[prediction_B] - dataframe[main_column_B]
    plt.figure(figsize=(4, 4))
    sns.set(style="whitegrid", font="Times New Roman", rc={"grid.linewidth": 0, "axes.grid": False})
    # Create the KDE plot of residuals
    g = sns.kdeplot(data=dataframe, x='residuals_A', y='residuals_B', kind='kde', fill=False, cmap='coolwarm')
    # Set x-axis and y-axis labels
    plt.axvline(x=0, linestyle='--',alpha = 0.5)
    plt.axhline(y=0, linestyle='--',alpha = 0.5)
    # Remove grid lines
    # Set font to Times New Roman
    g.set_xlabel(x_label)
    g.set_ylabel(y_label)


    #g.ax_joint.set_xlim(-2, 2)
    #g.ax_joint.set_ylim(-2, 2)

    # Save the plot as SVG file
    plt.savefig(save_name, format='svg', bbox_inches='tight')
    plt.show()

    # Set 'Times New Roman' as the default font
    mpl.rcParams['font.family'] = 'Times New Roman'
    mpl.rcParams['font.size'] = 16


 
plot_residuals(pred_df, "external_cons","total_steam","external_cons_pred","total_steam_pred", x_label = "External consumpion residua [t/h]", y_label="total steam residua [t/h]", save_name="ext.svg")

plot_residuals(pred_df, "W_trans_0_ANN","total_steam","W_trans_pred","total_steam_pred", x_label = "Turbine inflow residua [t/h]", y_label="total steam residua [t/h]", save_name="turb.svg")

plot_residuals(pred_df, "total_steam","W_trans_0_ANN","total_steam_pred_min","W_trans_pred")

calculate_proportion_below_real(pred_df, "total_steam", "total_steam_pred_min")
calculate_proportion_below_real(pred_df, "W_gen_0_ANN", "W_plan")

plot_residuals(pred_df, "W_gen_0_ANN","W_trans_0_ANN","W_plan","W_trans_pred")
plot_residuals(pred_df, "W_gen_0_ANN","W_trans_0_ANN","W_plan","W_trans_pred")



plot_residuals(pred_df, "mst_TG","W_gen_0_ANN","mst_TG_pred","W_gen_min_ANN")
#%%
def predict_general2(model, scaler, real_data_df):
    # Get the columns from the dataframe that were used in training
    cols_used_in_training = scaler.get_feature_names_out().tolist()
    X = real_data_df[cols_used_in_training].reset_index(drop=True)
    X_scaled = scaler.transform(X)
    # Convert the scaled data back into a DataFrame
    X_scaled = pd.DataFrame(X_scaled, columns=cols_used_in_training)
    X_scaled["const"] = 1
    
    # Make predictions
    y_pred = model.predict(X_scaled)


    return y_pred
real_data_X
def add_lagged_columns(df, column_name, start=24, end=48):
    for i in range(start, end+1):
        col_name = f"{column_name}_lag{i}"
        df[col_name] = df[column_name].shift(periods=i)
    return df
val_start = '2015-01-10'
val_end = '2015-12-30'
model_names = ["_LGBM","_ANN","_LR"]

df_valid_plot = df_ready["2016"].copy()
df_results = pd.DataFrame()

for cls in classes[:-2]:
    df_valid_plot = add_lagged_columns(df_ready, cls.column_name, start=24, end=48)
df_valid_plot = df_valid_plot.dropna()
for cls in classes[:]:
    df_results[cls.column_name] = df_valid_plot[cls.column_name]
    for model_name in ["_LGBM", "_ANN", "_LR"]:
        #for model_name in ["_LGBM", "_ANN", "_LR"]:
        col_name = cls.column_name + model_name
        print(col_name)
        scaler = cls.scaler_dict[col_name]
        model = cls.model_dict[col_name]
        if '_ANN' in model_name:
            y_pred = predict_general(model,scaler, df_valid_plot)
            df_results[col_name] = y_pred
        if '_LR' in model_name:
            if col_name == "W_gen_LR":
                y_pred = predict_general(model,scaler, df_valid_plot)
                df_results[col_name] = y_pred.tolist()
            else:
                y_pred = predict_general2(model,scaler, df_valid_plot)
                df_results[col_name] = y_pred.tolist()
        if '_LGBM' in model_name:
            y_pred = predict_general(model,scaler, df_valid_plot)
            df_results[col_name] = y_pred.tolist()


for cls in classes[:-2]:
    df_results[cls.column_name] = df_ready[val_start:val_end][cls.column_name]

for cls in classes[-2:]:
    print(cls.column_name)








model_A = "total_steam_LGBM_min"
model_B = "total_steam_LR_max"
model_C = "total_steam_ANN_max"

df_results["total_steam_LGBM_min"] = predict_general(classes[0].model_dict[model_A],classes[0].scaler_dict[model_A], df_valid_plot)
df_results["total_steam_ANN_min"] = predict_general(classes[0].model_dict[model_C],classes[0].scaler_dict[model_C], df_valid_plot)
df_results["total_steam_LR_min"] = predict_general2(classes[0].model_dict[model_B],classes[0].scaler_dict[model_B], df_valid_plot).tolist()

columns_dict = {
    "total_steam": {
        "predictions": ["total_steam_LGBM", "total_steam_ANN", "total_steam_LR"],
        "units": r'$\mathregular{total\;steam\;generated\;[t/h]}$',
        "legend_lables": [r'$\mathregular{\dot{m}_{st,total}}$', r'$\mathregular{\dot{m}_{st,total,LGBM}}$', r'$\mathregular{\dot{m}_{st,total,ANN}}$', r'$\mathregular{\dot{m}_{st,total,LR}}$'],
        "y_limit" : (95,150)
    },
    "external_cons": {
        "predictions": ["external_cons_LGBM", "external_cons_ANN", "external_cons_LR"],
        "units": r'$\mathregular{external\;consumption\;[t/h]}$',
        "legend_lables": [r'$\mathregular{\dot{m}_{ext}}$', r'$\mathregular{\dot{m}_{ext,LGBM}}$', r'$\mathregular{\dot{m}_{ext,ANN}}$', r'$\mathregular{\dot{m}_{ext,LR}}$'],
        "y_limit" : (4,10)
    },
    "m_st_bo": {
        "predictions": ["m_st_bo_LGBM", "m_st_bo_ANN", "m_st_bo_LR"],
        "units": r'$\mathregular{blow-off\;steam\;[t/h]}$',
        "legend_lables": [r'$\mathregular{\dot{m}_{st,bo}}$', r'$\mathregular{\dot{m}_{st,bo,LGBM}}$', r'$\mathregular{\dot{m}_{st,bo,ANN}}$', r'$\mathregular{\dot{m}_{st,bo,LR}}$'],
        "y_limit" : (3,5.5)
    },
    "de_air": {
        "predictions": ["de_air_LGBM", "de_air_ANN", "de_air_LR"],
        "units": r'$\mathregular{deaerator\;[t/h]}$',
        "legend_lables": [r'$\mathregular{\dot{m}_{de,air}}$', r'$\mathregular{\dot{m}_{de,air,LGBM}}$', r'$\mathregular{\dot{m}_{de,air,ANN}}$', r'$\mathregular{\dot{m}_{de,air,LR}}$'],
        "y_limit" : (3,7)
    },
    "mst_sf_demi": {
        "predictions": ["mst_sf_demi_LGBM", "mst_sf_demi_ANN", "mst_sf_demi_LR"],
        "units": r'$\mathregular{pre-heating\;[t/h]}$',
        "legend_lables": [r'$\mathregular{\dot{m}_{st,demi}}$', r'$\mathregular{\dot{m}_{st,demi,LGBM}}$', r'$\mathregular{\dot{m}_{st,demi,ANN}}$', r'$\mathregular{\dot{m}_{st,demi,LR}}$'],
        "y_limit" : (2,5)
    }
}



# Set the font to be Times New Roman for all text
mpl.rcParams['font.family'] = 'Times New Roman'


def plot_predictions(df, main_column, prediction_columns, labels, start_date, end_date, value_axis, legend_labels, line_plot_y_range):
    # Filter the dataframe for the specific interval for plotting
    df_plot = df.loc[start_date:end_date].copy()
    df_plot = df_plot.dropna(subset=[main_column])

    # Replace Inf values with a large number
    df_plot = df_plot.replace([np.inf, -np.inf], np.finfo(np.float64).max)
    # Set up the grid
    gs = gridspec.GridSpec(2, 3, height_ratios=[2, 1])
    fig = plt.figure(figsize=(8, 4.5))

    # Plot the main column and prediction columns in the first subplot
    ax0 = plt.subplot(gs[0, :])  # spans all columns of the first row
    colors = ['#000000', '#568CD6', '#f1b460', '#ff4040']
    ax0.plot(df_plot.index, df_plot[main_column], label=legend_labels[0], color=colors[0], linewidth=1.5, linestyle='--')
    for i, pred_col in enumerate(prediction_columns):
        ax0.plot(df_plot.index, df_plot[pred_col], label=legend_labels[i+1], color=colors[i+1],alpha = 0.70)
    ax0.set_ylabel(value_axis)
    # Set y-axis limits for "coldA" column using quantiles
    quantiles = df_plot[main_column].quantile([0.01, 0.99])  # Example quantiles: 5th and 95th percentiles
    ax0.set_ylim(line_plot_y_range[0], line_plot_y_range[1])  # Set y-axis limits
    ax0.legend(loc="lower left")
    ax0.grid(False)  # Removes grid lines

    # Calculate prediction metrics for the entire dataframe
    metrics = [r'$\mathregular{R^2\;Score}$', 'MAE', 'MRE']
    metric_values = []
    true_values = df[main_column]
    for pred_col in prediction_columns:
        pred_values = df[pred_col]
        r2 = np.round(r2_score(true_values, pred_values), 3)
        mae = np.round(mean_absolute_error(true_values, pred_values), 3)
        relative_error = np.round(np.mean(np.abs((true_values - pred_values) / true_values)), 3)
        metric_values.append([r2, mae, relative_error])

    # Create a subplot for each metric
    for i, metric in enumerate(metrics):
        ax = plt.subplot(gs[1, i])  # one subplot per column in the second row
        bars = ax.bar(labels, [value[i] for value in metric_values], color=colors[1:])
        ax.set_ylabel(metric)
        ax.grid(False)
        ax.spines['top'].set_visible(False)  # Hide top spine
        ax.spines['right'].set_visible(False)  # Hide right spine
        ax.spines['bottom'].set_visible(True)  # Hide bottom spine
        ax.spines['left'].set_visible(True)  # Hide left spine
        ax.tick_params(axis='both', which='both', length=0)  # Hide tick lines
        plt.xticks(rotation=0)  # Rotate x-axis labels
        if i == 0:  # Only for the first bar
            ax.set_ylim(0, 1)  # Set y-axis limits for the first bar
        for bar in bars:
            height = bar.get_height()
            ax.text(bar.get_x() + bar.get_width() / 2, height, f'{height:.2f}',
                    ha='center', va='bottom', color='black', fontsize=10)

    # Adjust layout
    plt.tight_layout()
    save_name = f"_{main_column}_score.svg"
    plt.savefig(save_name, format='svg', bbox_inches='tight')
    # Show the plot
    plt.show()



# Example usage
start_date = '2016-04-14'
end_date = '2016-04-20'
main_column = 'total_steam'
prediction_columns = ['total_steam_pred', 'total_steam_pred_max', 'total_steam_pred_min']
labels = ['LGBM', 'ANN', 'LR']  # Custom labels for bar plot

#plot_predictions(pred_df, main_column, prediction_columns, labels, start_date, end_date)


#for main_column, prediction_columns in columns_dict.items():
#    print("Main Column:", main_column)
#    print("Prediction Columns:", prediction_columns)
#    plot_predictions(df_results, main_column, prediction_columns, labels, start_date, end_date)



for main_column, info in columns_dict.items():
    prediction_columns = info["predictions"]
    units = info["units"]
    legend_lables = info["legend_lables"]
    y_limit = info["y_limit"]
    print("Main column:", main_column)
    print("Prediction columns:", prediction_columns)
    print("Units:", units)
    plot_predictions(df_results, main_column, prediction_columns, labels, start_date, end_date,units,legend_lables,y_limit)

columns_dict = {
"total_steam": {
    "predictions": ["total_steam_LGBM_min", "total_steam_ANN_min", "total_steam_LR_min"],
    "units": r'$\mathregular{total\;steam\;generated\;[t/h]}$',
    "legend_lables": [r'$\mathregular{\dot{m}_{st,total}}$',r'$\mathregular{\dot{m}_{st,total,LGBM}}$', r'$\mathregular{\dot{m}_{st,total,ANN}}$', r'$\mathregular{\dot{m}_{st,total,LR}}$'],
    "y_limit" : (75,110)
}
}   

for main_column, info in columns_dict.items():
    prediction_columns = info["predictions"]
    units = info["units"]
    legend_lables = info["legend_lables"]
    y_limit = info["y_limit"]
    print("Main column:", main_column)
    print("Prediction columns:", prediction_columns)
    print("Units:", units)
    plot_predictions(df_results, main_column, prediction_columns, labels, start_date, end_date,units,legend_lables,y_limit)


for main_column, prediction_columns in columns_dict.items():
    print("Main Column:", main_column)
    print("Prediction Columns:", prediction_columns)
    plot_predictions(df_results, main_column, prediction_columns, labels, start_date, end_date)




columns_dict = {
    "W_gen": {
        "predictions": ["W_gen_LGBM", "W_gen_ANN", "W_gen_LR"],
        "units": r'$\mathregular{power\;generated\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{gen}}$', r'$\mathregular{W_{gen,LGBM}}$', r'$\mathregular{W_{gen,ANN}}$', r'$\mathregular{W_{gen,LR}}$'],
        "y_limit" : (5,15.5)
    },
    "W_trans": {
        "predictions": ["W_trans_LGBM", "W_trans_ANN", "W_trans_LR"],
        "units": r'$\mathregular{power\;tranfered\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{trans}}$', r'$\mathregular{W_{trans,LGBM}}$', r'$\mathregular{W_{trans,ANN}}$', r'$\mathregular{W_{trans,LR}}$'],
        "y_limit" : (4.5,12.5)
    }
}   

for main_column, info in columns_dict.items():
    prediction_columns = info["predictions"]
    units = info["units"]
    legend_lables = info["legend_lables"]
    y_limit = info["y_limit"]
    print("Main column:", main_column)
    print("Prediction columns:", prediction_columns)
    print("Units:", units)
    plot_predictions(df_results, main_column, prediction_columns, labels, start_date, end_date,units,legend_lables,y_limit)


columns_dict = {
    "W_gen": {
        "predictions": ["W_gen_LGBM", "W_gen_ANN", "W_gen_LR"],
        "units": r'$\mathregular{power\;generated\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{gen}}$', r'$\mathregular{W_{gen,LGBM}}$', r'$\mathregular{W_{gen,ANN}}$', r'$\mathregular{W_{gen,LR}}$']
    },
    "W_trans": {
        "predictions": ["W_trans_LGBM", "W_trans_ANN", "W_trans_LR"],
        "units": r'$\mathregular{power\;tranfered\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{trans}}$', r'$\mathregular{W_{trans,LGBM}}$', r'$\mathregular{W_{trans,ANN}}$', r'$\mathregular{W_{trans,LR}}$']
    },
    "W_gen_0_ANN": {
        "predictions": ["W_plan", "W_gen_0_ANN_pred", "W_gen_min_ANN"],
        "units": r'$\mathregular{power\;tranfered\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{trans}}$', r'$\mathregular{W_{trans,LGBM}}$', r'$\mathregular{W_{trans,ANN}}$', r'$\mathregular{W_{trans,LR}}$']
    }
}   

columns_dict = {    "W_gen_0_ANN": {
        "predictions": ["W_plan", "W_gen_0_ANN_pred", "W_gen_min_ANN"],
        "units": r'$\mathregular{power\;tranfered\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{trans}}$', r'$\mathregular{W_{trans,LGBM}}$', r'$\mathregular{W_{trans,ANN}}$', r'$\mathregular{W_{trans,LR}}$']
    }
}   
for main_column, info in columns_dict.items():
    prediction_columns = info["predictions"]
    units = info["units"]
    legend_lables = info["legend_lables"]
    print("Main column:", main_column)
    print("Prediction columns:", prediction_columns)
    print("Units:", units)
    plot_predictions(pred_df, main_column, prediction_columns, labels, start_date, end_date,units,legend_lables)



calculate_proportion_below_real(pred_df, "W_gen_0_ANN", "W_plan")


#%%
import pandas as pd
import locale
def calculate_difference(df, colA, colB, coef):
    # Calculate the sum of colA and colB
    sum_values = df[colA] + df[colB]
    
    # Calculate the difference between colA and colB
    difference = df[colA] - df[colB]
    
    # Multiply the difference by the coefficient
    result = difference * coef
    
    # Sum the result
    total = result.sum()
    
    # Set locale to Czech for financial formatting
    locale.setlocale(locale.LC_ALL, 'cs_CZ.utf8')
    
    # Format the total as financial value in CZK with separation of thousands and millions
    formatted_total = locale.format_string('%.0f Kč', total, grouping=True)
    
    # Return the formatted total
    return formatted_total

result = calculate_difference(pred_df, 'W_gen_0_ANN', 'W_plan', df_cash["average_price (Kč/MWh)"].mean())
print(result)
result = calculate_difference(pred_df, 'W_plan', 'W_gen_min_ANN', df_cash["average_price (Kč/MWh)"].mean())
print(result)


import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
sns.set(font='Times New Roman')

# Calculate the quantiles
quantiles = [0.33, 0.66, 0.99]
quantile_names = ['Q33', 'Q66', 'Q99']
colors = ['#568CD6', '#f1b460']

results = []

for q, q_name in zip(quantiles, quantile_names):
    quantile_value = preds_df['W_trans_0_ANN'].quantile(q)
    mask = preds_df['W_trans_0_ANN'] >= quantile_value

    mae_ANN = mean_absolute_error(preds_df.loc[mask, 'W_trans_0_ANN'], preds_df.loc[mask, 'W_trans_0_ANN'])
    mae_LGBM = mean_absolute_error(preds_df.loc[mask, 'W_trans_0_ANN'], preds_df.loc[mask, 'W_trans_0_LGBM'])

    results.append({'quantile': q_name, 'model': 'W_trans_0_ANN', 'mae': mae_ANN})
    results.append({'quantile': q_name, 'model': 'W_trans_0_LGBM', 'mae': mae_LGBM})

mae_quantiles_df = pd.DataFrame(results)

# Set style and font
sns.set(style="whitegrid", font='Times New Roman')

plt.figure(figsize=(8, 4))

# Create a bar plot
mae_quantiles_plot = sns.barplot(x='quantile', y='mae', hue='model', data=mae_quantiles_df,palette=colors)

# Set plot labels and title
plt.xlabel('Quantile', fontsize=14)
plt.ylabel('Mean Absolute Error (MAE)', fontsize=14)
plt.title('Comparison of ANN and LGBM Model Performance at Different Quantiles of LR Model', fontsize=16, fontweight='bold')
plt.grid(False)
# Save the plot as SVG
plt.savefig('FigTwo.svg', format='svg', bbox_inches='tight')

# Show the plot
plt.show()


#%%


import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
sns.set(font='Times New Roman')

# Split LR values into three quantiles
q33 = pred_df['W_trans_0_ANN'].quantile(0.33)
q66 = pred_df['W_trans_0_ANN'].quantile(0.66)
q99 = pred_df['W_trans_0_ANN'].quantile(0.99)
colors = ['#568CD6', '#f1b460']

# Filter data for each quantile
data_q33 = pred_df[pred_df['W_trans_0_ANN'] <= q33]
data_q66 = pred_df[(pred_df['W_trans_0_ANN'] > q33) & (pred_df['W_trans_0_ANN'] <= q66)]
data_q99 = pred_df[pred_df['W_trans_0_ANN'] > q66]

# Calculate MSE for each quantile and model
mse_results = []
for quantile, data in zip(['Q33', 'Q66', 'Q99'], [data_q33, data_q66, data_q99]):
    for model in ['W_trans_0_LGBM', 'W_trans_0_ANN']:
        errors = data[model] - data['W_trans_0_ANN']
        mse = np.mean(errors)
        mse_results.append({'quantile': quantile, 'model': model, 'mse': mse})

# Convert the results to a DataFrame
mse_df = pd.DataFrame(mse_results)

# Plotting
sns.set(style="whitegrid")

plt.figure(figsize=(6, 4))
ax = sns.barplot(x='quantile', y='mse', hue='model', data=mse_df, palette=colors)
ax.axhline(0, color='gray', linestyle='--')  # Add a grid line at 0
ax.set_title('Signed Mean  Error for ANN and LGBM in Different Quantiles\nPower Transferred', fontweight='bold', fontname='Times New Roman')

ax.set_ylabel('SME', fontname='Times New Roman')
ax.set_xlabel('Quantiles', fontname='Times New Roman')
plt.legend(loc='upper right', prop={'family': 'Times New Roman'})
# Remove the original legend
ax.get_legend().remove()

# Add a custom legend with specified labels
handles, _ = ax.get_legend_handles_labels()
ax.legend(handles, ["LGBM Model", "ANN model"], loc='upper right', prop={'family': 'Times New Roman'})

plt.grid(False)
plt.savefig('FigTwo.svg', format='svg',bbox_inches='tight')
plt.show()













columns_dict = {    "W_gen_0_LR": {
        "predictions": ["W_gen_0_ANN_pred", "W_gen_min_ANN"],
        "units": r'$\mathregular{power\;tranfered\;[MWh]}$',
        "legend_lables": [r'$\mathregular{W_{trans}}$', r'$\mathregular{W_{trans,LGBM}}$', r'$\mathregular{W_{trans,ANN}}$', r'$\mathregular{W_{trans,LR}}$']
    }
}   
for main_column, info in columns_dict.items():
    prediction_columns = info["predictions"]
    units = info["units"]
    legend_lables = info["legend_lables"]
    print("Main column:", main_column)
    print("Prediction columns:", prediction_columns)
    print("Units:", units)
    plot_predictions(pred_df, main_column, prediction_columns, labels, start_date, end_date,units,legend_lables)



def plot_predictions(df, ann_col, lr_col, lgbm_col, title):
    fig, ax = plt.subplots(figsize=(8, 4))
    ann_color = '#568CD6'
    lr_color = '#D656CC'
    lgbm_color = '#f1b460'
    # Plot the predictions for each model
    sns.lineplot(data=df, x=df.index, y=df[ann_col], label='ANN model', color=ann_color, alpha=0.7)
    sns.lineplot(data=df, x=df.index, y=df[lr_col], label='LR model', color=lr_color, linestyle='--', linewidth=2)
    sns.lineplot(data=df, x=df.index, y=df[lgbm_col], label='LGBM model', color=lgbm_color, alpha=0.7)

    # Set plot labels and title
    ax.set_ylabel('Power transfered [MWh]', fontsize=12, fontname='Times New Roman')
    ax.set_title(title, fontsize=14, fontname='Times New Roman')
    ax.grid(False)
    
    # Add legend and show the plot
    ax.legend(fontsize=12)
    plt.savefig('Figtrwew.svg', format='svg', bbox_inches='tight')
    plt.show()

# Example usage:
# Replace 'pred_df' with your DataFrame containing the predictions, and
# adjust the column names accordingly
plot_predictions(pred_df[1000:1300], 'W_trans_0_ANN', 'W_trans_0_ANN', 'W_trans_0_LGBM',  r'$\mathregular{Power transfered\; (W_{transfered})\;model\;comparison}$')




import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from iapws import IAPWS95
import statsmodels.api as sm
from matplotlib.dates import DateFormatter, AutoDateLocator
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose
#ff4040
colors = ['#56D660', '#568CD6', '#ff4040', '#f1b460']
def Q_DHS_plot(df, Q_req, Q_DHS, m_st_bp, x_label, y_label_1, y_label_2, title, fig_height=3.3):
    # Calculate figure width for A4 paper with standard margins
    # A4 width = 8.27 inches, assuming 1 inch for margins
    fig_width = 8.27 - 1 
    T_st_T = 70 + 273.15  # convert temperature to Kelvin (K)
    water_T = IAPWS95(T=T_st_T, x=0)
    steam_T = IAPWS95(T=235+273.15, P=1.22)
    h_steam = steam_T.h - water_T.h
    # Set seaborn style
    sns.set(style="whitegrid", font="Times New Roman", rc={"grid.linewidth": 0, "axes.grid": False})
    plt.rcParams['xtick.bottom'] = True
    sns.despine()
    # Compute moving averages
    Q_req_ma = df[Q_req].rolling(window=24).mean()
    Q_DHS_ma = df[Q_DHS].rolling(window=24).mean()

    m_st_bp_MW = df[m_st_bp] * h_steam * 0.27778 / 1000
    m_st_bp_ma = m_st_bp_MW.rolling(window=24).mean()
    
    Q_DHS_adjusted = Q_DHS_ma - m_st_bp_ma

    # Create figure and axis objects with subplots()
    fig,ax = plt.subplots(figsize=(fig_width, fig_height))
    ax.plot(df.index, Q_req_ma, label=r'$\mathregular{\dot{Q}_{contract,extrapolated}}$', color='#56d660')
    ax.plot(df.index, Q_DHS_adjusted, label=r'$\mathregular{\dot{Q}_{DHS,adjusted}}$', color='#568CD6')

    # Make a plot with a separate y-axis for m_st_bp using twinx()
    ax2=ax.twinx()
    ax2.plot(df.index, m_st_bp_ma, label=r'$\mathregular{\dot{Q}_{bypass}}$', color='#f1b460', linestyle='dashed')

    x_label = r''
    y_label_1 = r'Heat [MWh]'
    y_label_2 = r'Bypass heat contribution [MWh]'
    # Set labels and title
    ax.set_xlabel(x_label, fontname='Times New Roman', fontsize=12)
    ax.set_ylabel(y_label_1, fontname='Times New Roman', fontsize=12)
    ax2.set_ylabel(y_label_2, fontname='Times New Roman', fontsize=12)
    plt.title("Extrapolated Heat contract and Adjusted heat delivery", fontname='Times New Roman', fontsize=16)

    # Legend 
    fig.legend(loc="upper left", bbox_to_anchor=(0,1), bbox_transform=ax.transAxes)
    ax.tick_params(axis='x', rotation=15)
    # Rotate x-axis labels
    plt.xticks(rotation=45)

    # Save and show the plot
    plt.savefig("EST.svg", format='svg', bbox_inches='tight')
    plt.show()

#Q_DHS_plot(df_ready["2016"][-6250:2000], "Q_req","Q_DHS", "lin_bp", "Hello kitti", "Test me ", "Im so pretty", "title")


def plot_boiler_activity(df, threshold=20):
    # For the boxplot, filter rows where boiler values are above the threshold
    active_boilers = df[['Boiler_1', 'Boiler_2', 'Boiler_3', 'Boiler_4']].apply(lambda x: x[x > threshold])
    fig_width = 8.27 - 1 
    fig_height = 2.73
    plt.figure(figsize=(fig_width, fig_height))
    colors = ['#568CD6', '#56D660', '#ff4040', '#f1b460']
    boxprops = {'edgecolor': 'white'}  # Set edgecolor to white
    meanprops = {'marker': 'o', 'markerfacecolor': 'white', 'markeredgecolor': 'white'}  # Set mean properties to white
    sns.boxplot(data=active_boilers, orient='v', palette=colors, boxprops=boxprops, meanprops=meanprops)  # Pass boxprops and meanprops arguments

    plt.ylabel(r'$\mathregular{\dot{m}_{st,production}}$[t/h]', fontsize=12)
    plt.title('Boiler Values Above Threshold', fontsize=12)
    
    # Set the X axis labels
    plt.xticks(range(len(active_boilers.columns)), ['Boiler 1', 'Boiler 2', 'Boiler 3', 'Boiler 4'], fontsize=12)
    
    plt.savefig("tvoje_ma22ma.svg", format='svg', bbox_inches='tight')
    plt.show()

    # Ignore the rows where boilers_active is 0
    df = df[df['boilers_active'] > 0]

    # Create pie charts for each season
    fig, axs = plt.subplots(2, 2, figsize=(fig_width, fig_width))  # Create 2x2 subplots
    seasons = ['Winter', 'Spring', 'Summer', 'Fall']
    season_colors = {'Winter': '#8bb3e6', 'Spring': '#90d17f', 'Summer': '#ff9a87', 'Fall': '#ffa04c'}
    
    for ax, season in zip(axs.flat, seasons):
        season_df = df[df[season] == 1]  # Filter data for the current season
        counts = season_df['boilers_active'].value_counts(normalize=True)  # Get percentage of each count
        ax.pie(counts, colors=colors, startangle=140)
        ax.set(ylabel='', title=f'{season} - % of Activation Cases', aspect='equal')  # Equal aspect ratio ensures that pie is drawn as a circle.
        ax.title.set_size(14)  # Set the title's font size
    ax.legend(counts.index, title="Boilers Active", loc="best")  # Add legend
    
    plt.tight_layout()
    plt.savefig("tvoje_ma233ma.svg", format='svg', bbox_inches='tight')
    plt.show()
plot_boiler_activity(df_ready)

def plot_seasonal_decomposition(df, columns):
    fig_width = 8.27 - 1
    for i, col in enumerate(columns):
        # Convert column to time series
        ts = pd.Series(df[col].values, index=df.index)

        # Perform seasonal decomposition
        result = sm.tsa.seasonal_decompose(ts, model='additive', period=24)

        # Plot decomposition
        plt.figure(figsize=(fig_width, 16))
        plt.subplot(611)
        plt.plot(ts, label='Original', color='#568CD6')
        plt.legend(loc='upper left', fontsize=14)
        plt.tick_params(axis='both', which='both', labelsize=14)
        plt.gca().axes.get_xaxis().set_visible(False)  # Hide X-axis
        plt.subplot(612)
        plt.plot(result.trend, label='Trend', color='#568CD6')
        plt.legend(loc='upper left', fontsize=14)
        plt.tick_params(axis='both', which='both', labelsize=14)
        plt.gca().axes.get_xaxis().set_visible(False)  # Hide X-axis
        plt.subplot(613)
        plt.plot(result.seasonal, label='Seasonality', color='#568CD6')
        plt.legend(loc='upper left', fontsize=14)
        plt.tick_params(axis='both', which='both', labelsize=12)
        plt.gca().axes.get_xaxis().set_visible(False)  # Hide X-axis
        plt.subplot(614)
        plt.plot(result.resid, label='Residuals', color='#568CD6')
        plt.legend(loc='upper left', fontsize=14)
        plt.tick_params(axis='both', which='both', labelsize=12)
        plt.gca().xaxis.set_major_formatter(DateFormatter('%d-%m'))
        plt.savefig(f'{chr(67 + i)}.svg', format='svg',bbox_inches='tight')
        plot_acf(ts, lags=48, alpha=0.05)  # Plot autocorrelation with lag 24
        plt.tick_params(axis='both', which='both', labelsize=14)
        plt.savefig(f'{chr(66 + i)}.svg', format='svg',bbox_inches='tight')
        plot_pacf(ts, lags=48, alpha=0.05)  # Plot partial autocorrelation with lag 24
        plt.tick_params(axis='both', which='both', labelsize=14, rotation=0)
        plt.suptitle(col, fontsize=14)
        
        # Format X-axis as dd-mm
        plt.gca().xaxis.set_major_formatter(DateFormatter('%d-%m'))
        plt.gca().xaxis.set_major_locator(AutoDateLocator(maxticks=5))

        # Save the figures as A.svg, B.svg, C.svg, etc.
        plt.savefig(f'{chr(65 + i)}.svg', format='svg',bbox_inches='tight')

# Usage
plot_seasonal_decomposition(df_ready["5/1/2016":"5/7/2016"], main_columns[:1])
plot_seasonal_decomposition(df_ready["5/1/2016":"5/7/2016"], main_columns[2:3])

import matplotlib.pyplot as plt
sns.set(style="whitegrid", font="Times New Roman", rc={"grid.linewidth": 0, "axes.grid": False})
def plot_train_validation_test_splits(df, column_name, train_end, val_start, val_end, test_start):
    # Split the dataframe into train, validation, and test sets
    train = df[:train_end]
    validation = df[val_start:val_end]
    test = df[test_start:]

    # Set up the figure and axes
    plt.figure(figsize=(8.27, 4))
    ax = plt.gca()

    # Plot the train, validation, and test sets with colors
    train[column_name].plot(ax=ax, label='Train', color='#568CD6')
    validation[column_name].plot(ax=ax, label='Validation', color='#56D660')
    test[column_name].plot(ax=ax, label='Test', color='#f1b460')

    # Remove colors for data in the year 2017

    # Set the color for data in the year 2017
    test['2017'][column_name].plot(ax=ax, label='Removed', color='#ff4040')

    # Set labels and title
    plt.xlabel('', fontsize=12)
    plt.ylabel("Power transferred [MWh]", fontsize=12)
    plt.title('Train, Validation, and Test Splits')

    # Set legend font size
    plt.legend(fontsize=12,ncol = 2)
    plt.tick_params(axis='both', which='both', labelsize=12, rotation=0)

    # Show the plot
    plt.savefig("traintestsplit.svg", format='svg', bbox_inches='tight')
    plt.show()

# Usage
# Define train end, validation start and end, and test start dates
train_end = '2015-01-06'
val_start = '2015-01-07'
val_end = '2015-12-31'
test_start = '2016-01-01'
# Call the function with these dates
plot_train_validation_test_splits(df_ready, "W_trans", train_end, val_start, val_end, test_start)














