sf_runoff_lt.py

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

from sklearn.svm import SVR
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.compose import TransformedTargetRegressor
from sklearn.model_selection import GridSearchCV,TimeSeriesSplit
from sklearn.metrics import mean_squared_error


import matplotlib.pyplot as plt

import os

import pdb

def shift_series_(s, shift_range,t_unit):

    s_shifts = [s.shift(-t_unit * shift, freq='D').rename(f'{s.name}_{shift}') for shift in range(*shift_range)]
    return pd.concat(s_shifts, axis=1)


def create_it_matrix_lt(daily_input, t_length,t_unit, lead_t):

    # This function takes as input the daily temperature, precipitation and runoff and generates the input-target matrix

    # Read the daily input and extract runoff, evaporation, temperature and precipitation dataframe
    if isinstance(daily_input, str):
        daily_input = pd.read_csv(daily_input, index_col=0, parse_dates=True)
    runoff = daily_input[['Q']]
    temp = daily_input[[c for c in daily_input.columns if c[0] == 'T']]
    prec = daily_input[[c for c in daily_input.columns if c[0] == 'P']]
    evap = daily_input[[c for c in daily_input.columns if c[0] == 'E']]


    # Compute the t_unit days average runoff
    runoff_t_unit = runoff.rolling(t_unit, min_periods=t_unit).mean()
    #runoff_t_unit_shifted = runoff_t_unit.shift(t_unit * lead_t,freq='D')
    #pdb.set_trace()
    
    # Compute the t_unit days average temperature
    if not temp.empty:
        temp_t_unit = temp.rolling(t_unit, min_periods=t_unit).mean()
        temp_t_unit = pd.concat([shift_series_(temp_t_unit.loc[:, col], (-t_length + 1, 1-lead_t),t_unit) for col in temp_t_unit], axis=1)

    # Compute the t_unit days sum precipitation
    if not prec.empty:
        prec_t_unit = prec.rolling(t_unit, min_periods=t_unit).sum()
        prec_t_unit = pd.concat([shift_series_(prec_t_unit.loc[:, col], (-t_length + 1, 1-lead_t),t_unit) for col in prec_t_unit], axis=1)
    
    # Compute the t_unit days sum evapotranspiration
    if not evap.empty:
        evap_t_unit = evap.rolling(t_unit, min_periods=t_unit).sum()
        evap_t_unit = pd.concat([shift_series_(evap_t_unit.loc[:, col], (-t_length + 1, 1-lead_t),t_unit) for col in evap_t_unit], axis=1)

    result = pd.concat([runoff_t_unit, temp_t_unit, prec_t_unit, evap_t_unit], axis=1).dropna()
    #result = result.shift(-t_unit * lead_t,freq='D')
    # Create the input-target matrix
    return result




"""
def sf_input_matrix(input_matrix, seasonal_forecast, lead_time, member):

    X = input_matrix.copy()
    X = X.loc[X.index.is_month_end, :]

    for offset in range(0, -lead_time, -1):
        sf_dates = X.index + pd.offsets.MonthBegin(offset - 1)
        X_columns = [c for c in X.columns if int(c.split('_')[1]) == offset]
        sf_columns = [f'SF_{c.split("_")[0]}_lt{lead_time + offset}_m{member}' for c in X_columns]
        X.loc[:, X_columns] = seasonal_forecast.loc[sf_dates, sf_columns].values

    return X
"""


def svr_gridSearch(daily_input, t_length,t_unit, C_range=np.logspace(-3, 1, 5), epsilon_range=np.logspace(-5, 0, 5),
                   plot=False, n_splits=8):

    # svr_gridSearch run the grid search on C and epsilon svr parameters.
    
    # Read the input-target matrix 
    it_matrix = create_it_matrix(daily_input, t_length)
    X = it_matrix.drop(columns='Q')
    y = it_matrix['Q']
    
    # Set up the splits respecting of the time-series nature of the dataset

    tscv = TimeSeriesSplit(gap=t_unit,n_splits=n_splits, test_size=365)
    tscv.split(X)
    
    #pdb.set_trace()
     
    # Set up the grid search parameters
    svr_estimator = SVR(kernel='rbf', gamma='scale', cache_size=1000)
    svr_estimator = make_pipeline(StandardScaler(),
                                  TransformedTargetRegressor(regressor=svr_estimator, transformer=StandardScaler()))
    parameters = {'transformedtargetregressor__regressor__C': C_range,
                  'transformedtargetregressor__regressor__epsilon': epsilon_range}
    svr_model = GridSearchCV(svr_estimator, cv=tscv, param_grid=parameters, n_jobs=-1, verbose=1, refit=True)


    # execute the grid search
    svr_model.fit(X, y)
    
    
    # Select the best C and epsilon
    best_C = svr_model.best_params_['transformedtargetregressor__regressor__C']
    best_epsilon = svr_model.best_params_['transformedtargetregressor__regressor__epsilon']
    print(f'For {t_length} months of data input: Best estimator: C={best_C}, epsilon={best_epsilon}')

    # Check if the best C (or epsion) is in the border of the grid
    if best_C == max(C_range) or best_C == min(C_range):
        print(f'Warning: best C found on the grid limit: C = {best_C}')
    if best_epsilon == max(epsilon_range) or best_epsilon == min(epsilon_range):
        print(f'Warning: best epsilon found on the grid limit: epsilon = {best_epsilon}')
    print()
    
    if plot:

        # Gridsearch plot
        plt.figure()
        scores = svr_model.cv_results_['mean_test_score'].reshape(
            len(svr_model.param_grid['transformedtargetregressor__regressor__C']),
            len(svr_model.param_grid['transformedtargetregressor__regressor__epsilon']))
        plt.imshow(scores)
        plt.colorbar()
        plt.xticks(np.arange(len(svr_model.param_grid['transformedtargetregressor__regressor__epsilon'])),
                   ['%.3f' % a for a in svr_model.param_grid['transformedtargetregressor__regressor__epsilon']], rotation=45)
        plt.yticks(np.arange(len(svr_model.param_grid['transformedtargetregressor__regressor__C'])),
                   ['%.3f' % a for a in svr_model.param_grid['transformedtargetregressor__regressor__C']])
        plt.xlabel('epsilon')
        plt.ylabel('C')
        plt.title(f'Heatmap for {t_length} months of input data')
        plt.tight_layout()

        # Scatterplot true vs. estimated training
        plt.figure()
        true_values = y
        est_values = svr_model.predict(X)
        xy = np.vstack([true_values, est_values])
        z = gaussian_kde(xy)(xy)
        idx = z.argsort()
        x, y, z = true_values[idx], est_values[idx], z[idx]
        plt.scatter(x, y, c=z)
        plt.plot([min(true_values), max(true_values)], [min(true_values), max(true_values)], color='black')
        plt.xlabel('True values')
        plt.ylabel('Estimated values')
        plt.title(f'Heatmap for {t_length} months of input data')
        plt.gca().set_aspect('equal', adjustable='box')
        plt.draw()
    return {'C': best_C, 'epsilon': best_epsilon, 'score': svr_model.best_score_,
            'best_estimator': svr_model.best_estimator_}


def feature_sel(daily_input):

    # Find the best t_length and temperature/precipitation stations. Return score for each combination

    gridSearch_results = []
    for t_len in range(1, 12):
        gridSearch_results.append(svr_gridSearch(daily_input, 
                                                 t_len,
                                                 C_range=np.logspace(-2, 2, 6), 
                                                 epsilon_range=np.logspace(-3, 0, 6),
                                                 plot=True))

    return pd.DataFrame(data=gridSearch_results, index=range(1, 12))

    # if isinstance(daily_input, str):
    #     daily_input = pd.read_csv('/home/mcallegari@eurac.edu/SECLI-FIRM/Mattia/SF_runoff/Zoccolo/inputs/daily_input.csv', index_col=0, parse_dates=True)
    #
    # temp_col = [c for c in daily_input.columns if c[0] == 'T']
    # prec_col = [c for c in daily_input.columns if c[0] == 'P']
    #
    # for feat in itertools.product(temp_col, prec_col):
    #     for t_len in range(1, 12):
    #         svr_gridSearch(daily_input.loc[:, ('Q',) + feat], t_len, plot=False)


def training(daily_input, t_length=None, svr_C=None, svr_epsilon=None):

    # This function takes as input the daily temperature, precipitation and runoff and train a model to predict the
    # monthly mean runoff using a time series of temperature and precipitation of length t_length as input features

    # Execute the feature selection if t_length is None
    if t_length is None:

        print('Finding the best input time series length...')
        fs = feature_sel(daily_input)
        fs_best = fs.loc[fs['score'].idxmax()]
        svr_C, svr_epsilon = fs_best['C'], fs_best['epsilon']
        t_length = fs_best.name
        print(f'Best time series length: {t_length}. '
              f'Best C={svr_C}, epsilon={svr_epsilon} (score={fs_best["score"]})')
        svr_model = fs_best['best_estimator']

    else:

        it_matrix = create_it_matrix(daily_input, t_length)
        X = it_matrix.drop(columns='Q')
        y = it_matrix['Q']

        # Fit the SVR with standardized input and target
        svr_estimator = SVR(kernel='rbf', C=svr_C, epsilon=svr_epsilon, gamma='scale', cache_size=1000)
        svr_model = make_pipeline(StandardScaler(),
                                  TransformedTargetRegressor(regressor=svr_estimator, transformer=StandardScaler()))
        svr_model.fit(X, y)

    return svr_model


def monthly_climatology(daily_input,t_unit):

    if isinstance(daily_input, str):
        daily_input = pd.read_csv(daily_input, index_col=0, parse_dates=True)

    monthly_mean_columns = [c for c in daily_input.columns if c[0] in ['Q', 'T']]
    monthly_mean = daily_input.loc[:, monthly_mean_columns].groupby(by=daily_input.index.month).mean()
    #remember to add ['E', 'P']
    monthly_sum_columns = [c for c in daily_input.columns if c[0] in ['P','E']]
    monthly_sum = daily_input.loc[:, monthly_sum_columns].groupby(by=daily_input.index.month).mean() * t_unit
    #pdb.set_trace()
    return pd.concat([monthly_mean, monthly_sum], axis=1)#[monthly_mean_columns,monthly_sum_columns]

def daily_climatology(daily_input,t_unit):
    
    runoff = daily_input[['Q']]
    temp = daily_input[[c for c in daily_input.columns if c[0] == 'T']]
    prec = daily_input[[c for c in daily_input.columns if c[0] == 'P']]
    evap = daily_input[[c for c in daily_input.columns if c[0] == 'E']]


    # Compute the t_unit days average runoff
    runoff_t_unit = runoff.rolling(t_unit, min_periods=t_unit).mean()

    
    # Compute the t_unit days average temperature
    if not temp.empty:
        temp_t_unit = temp.rolling(t_unit, min_periods=t_unit).mean()
        #temp_t_unit = pd.concat([shift_series_t_unitdays(temp_t_unit.loc[:, col], (-t_length + 1, 1)) for col in temp_t_unit], axis=1)

    # Compute the t_unit days sum precipitation
    if not prec.empty:
        prec_t_unit = prec.rolling(t_unit, min_periods=t_unit).sum()
        #prec_t_unit = pd.concat([shift_series_t_unitdays(prec_t_unit.loc[:, col], (-t_length + 1, 1)) for col in prec_t_unit], axis=1)
    
    # Compute the t_unit days sum evapotranspiration
    if not evap.empty:
        evap_t_unit = evap.rolling(t_unit, min_periods=t_unit).sum()
        #evap_t_unit = pd.concat([shift_series_t_unitdays(evap_t_unit.loc[:, col], (-t_length + 1, 1)) for col in evap_t_unit], axis=1)

    daily_t_unit = pd.concat([runoff_t_unit, temp_t_unit, prec_t_unit, evap_t_unit], axis=1)
    daily_mean = daily_t_unit.groupby(by=daily_t_unit.index.day_of_year).mean()

    #pdb.set_trace()
    return daily_mean




def loyo_cv_lc_nofor(daily_input,output_folder, t_length=None, svr_C=None, svr_epsilon=None,
               lead_time=range(1, 8)):

    # Compute the climatology for all the inputs
    if isinstance(daily_input, str):
        daily_input = pd.read_csv(daily_input, index_col=0, parse_dates=True)
    monthly_clim = monthly_climatology(daily_input)


    # Create the total input-target matrix and select the start and end year for the loo-cv
    it_matrix = create_it_matrix(daily_input, t_length)
    year_start = min(it_matrix.index[(it_matrix.index.month == 1) & (it_matrix.index.day == 31)].year)
    year_end = max(it_matrix.index[(it_matrix.index.month == 12) & (it_matrix.index.day == 31)].year)

    # For each training length and each year train a different model and test with different configurations
    for Nyears_training in range(1, year_end-year_start+1):
        prediction = []
        for year in range(year_start, year_end+1):
            print(f'Testing on year {year}; Training years: {Nyears_training} ...')

            # Drop the years that should not be considered in the training and train the model
            daily_input_loo = daily_input.copy()
            training_years = (list(range(year+1, year_end+1)) + list(range(year_start, year)))[-Nyears_training:]
            daily_input_loo.loc[[i for i in daily_input.index if i.year not in training_years], 'Q'] = np.nan
            # daily_input_loo.loc[daily_input.index.year == year, 'Q'] = np.nan
            svr_model = training(daily_input_loo, t_length, svr_C, svr_epsilon)

            # Select the dates on which to execute the prediction
            test_dates = daily_input.index[daily_input.index.is_month_end & (daily_input.index.year == year)]

            # Save the true runoff values
            y = {}
            y['true_runoff'] = daily_input.loc[daily_input.index.year == year, 'Q'].resample("M").mean().values

            # Compute runoff monthly climatology considering the subset used for training the svr
            y['runoff_clim'] = [daily_input_loo.loc[daily_input_loo.index.month == m, 'Q'].mean() for m in range(1, 13)]

            # Predict using true temperature and precipitation (no forecast)
            X_trueTP = it_matrix.loc[test_dates, :].drop(columns='Q')
            y['trueTP'] = svr_model.predict(X_trueTP)

             # Predict using temperature and precipitation climatology (no forecast) for all lead times
            X_climTP = X_trueTP.copy()
            for lt in lead_time:
                change_dest = [c for c in X_climTP.columns if c.split('_')[1] == str(-lt + 1)]
                change_source = [c.split('_')[0] for c in change_dest]
                X_climTP.loc[:, change_dest] = monthly_clim.loc[(test_dates - pd.offsets.MonthBegin(lt)).month, change_source].values
                y[f'climTP_lt{lt}'] = svr_model.predict(X_climTP)
                
            '''
            # Predict using seasonal forecast with different lead times
            for lt in lead_time:
                for m in range(1, 26):
                    X_sf = sf_input_matrix(X_trueTP, seasonal_forecast, lt, m)
                    y[f'sfTP_m{m}_lt{lt}'] = svr_model.predict(X_sf)
            '''
            # Store the results in a dataframe
            prediction.append(pd.DataFrame(data=y, index=test_dates-pd.offsets.MonthBegin()))
            
        pd.concat(prediction, axis=0).to_csv(
            os.path.join(output_folder, f'Runoff_forecast_{Nyears_training}_trainingyears.csv'))

def loyo_cv_lc(daily_input, seasonal_forecast, output_folder, t_length=None, svr_C=None, svr_epsilon=None,
               lead_time=range(1, 8)):

    # Compute the climatology for all the inputs
    if isinstance(daily_input, str):
        daily_input = pd.read_csv(daily_input, index_col=0, parse_dates=True)
    monthly_clim = monthly_climatology(daily_input)

    # Read the seasonal forecast of temperature and precipitation
    if isinstance(seasonal_forecast, str):
        seasonal_forecast = pd.read_csv(seasonal_forecast, index_col=0, parse_dates=True)

    # Create the total input-target matrix and select the start and end year for the loo-cv
    it_matrix = create_it_matrix(daily_input, t_length)
    year_start = min(it_matrix.index[(it_matrix.index.month == 1) & (it_matrix.index.day == 31)].year)
    year_end = max(it_matrix.index[(it_matrix.index.month == 12) & (it_matrix.index.day == 31)].year)

    # For each training length and each year train a different model and test with different configurations
    for Nyears_training in range(1, year_end-year_start+1):
        prediction = []
        for year in range(year_start, year_end+1):
            print(f'Testing on year {year}; Training years: {Nyears_training} ...')

            # Drop the years that should not be considered in the training and train the model
            daily_input_loo = daily_input.copy()
            training_years = (list(range(year+1, year_end+1)) + list(range(year_start, year)))[-Nyears_training:]
            daily_input_loo.loc[[i for i in daily_input.index if i.year not in training_years], 'Q'] = np.nan
            # daily_input_loo.loc[daily_input.index.year == year, 'Q'] = np.nan
            svr_model = training(daily_input_loo, t_length, svr_C, svr_epsilon)

            # Select the dates on which to execute the prediction
            test_dates = daily_input.index[daily_input.index.is_month_end & (daily_input.index.year == year)]

            # Save the true runoff values
            y = {}
            y['true_runoff'] = daily_input.loc[daily_input.index.year == year, 'Q'].resample("M").mean().values

            # Compute runoff monthly climatology considering the subset used for training the svr
            y['runoff_clim'] = [daily_input_loo.loc[daily_input_loo.index.month == m, 'Q'].mean() for m in range(1, 13)]

            # Predict using true temperature and precipitation (no forecast)
            X_trueTP = it_matrix.loc[test_dates, :].drop(columns='Q')
            y['trueTP'] = svr_model.predict(X_trueTP)

            # Predict using temperature and precipitation climatology (no forecast) for all lead times
            X_climTP = X_trueTP.copy()
            for lt in lead_time:
                change_dest = [c for c in X_climTP.columns if c.split('_')[1] == str(-lt + 1)]
                change_source = [c.split('_')[0] for c in change_dest]
                X_climTP.loc[:, change_dest] = monthly_clim.loc[(test_dates - pd.offsets.MonthBegin(lt)).month, change_source].values
                y[f'climTP_lt{lt}'] = svr_model.predict(X_climTP)
                

            # Predict using seasonal forecast with different lead times
            for lt in lead_time:
                for m in range(1, 26):
                    X_sf = sf_input_matrix(X_trueTP, seasonal_forecast, lt, m)
                    y[f'sfTP_m{m}_lt{lt}'] = svr_model.predict(X_sf)

            # Store the results in a dataframe
            prediction.append(pd.DataFrame(data=y, index=test_dates-pd.offsets.MonthBegin()))

        pd.concat(prediction, axis=0).to_csv(
            os.path.join(output_folder, f'Runoff_forecast_{Nyears_training}_trainingyears.csv'))


# -----------------------------------------------------------------------

# Read and plot results


def root_mean_squared_error(y_true, y_pred):
    return np.sqrt(mean_squared_error(y_true, y_pred))

def smape(A, F):
    return 100/len(A) * np.sum(2 * np.abs(F - A) / (np.abs(A) + np.abs(F)))


def monthly_rmse(fileName, plot=True):

    # Open the result file
    runoff = pd.read_csv(fileName, index_col=0, parse_dates=True)

    # Create the ensamble mean
    for lt in range(1, 8):
        columns = [c for c in runoff.columns if 'sfTP_m' in c and f'lt{lt}' in c]
        runoff[f'sfTP_em_lt{lt}'] = runoff.loc[:, columns].mean(axis=1)

    # Compute the RMSE for each month
    runoff_error = []
    for m in range(1, 12):
        runoff_m = runoff.loc[runoff.index.month == m, :]
        runoff_error.append(
            runoff_m.apply(lambda y_pred: root_mean_squared_error(runoff_m['true_runoff'], y_pred), axis=0)
        )

    runoff_error = pd.DataFrame(data=runoff_error, index=range(1, 12))

    if plot:
        plt.figure()
        for lt in range(1, 8):
            runoff_error.loc[:, f'sfTP_em_lt{lt}'].plot(marker='o', label=f'lead_time={lt}')
        runoff_error.loc[:, 'runoff_clim'].plot(marker='o', color='black', label='climatology', linewidth=3)
        runoff_error.loc[:, 'trueTP'].plot(marker='o', color='red', label='era5', linewidth=3)
        plt.legend()
        plt.ylabel('RMSE ($m^3/s$)')
        plt.xlabel('Month')


def plot_it_matrix(daily_input, var, common_ylim=True):

    # ## Plot each variable of the input-target matrix

    # Create the input-target matrix from the daily_input
    it_matrix = create_it_matrix(daily_input, 1).rename(columns=lambda c: c[:2])

    # Create the ylabel dictionary
    plt_ylabel = {'P': 'Total precipitation (m)', 'T': 'Mean temperature (K)', 'Q': 'Runoff ($m^3/s$)'}

    # Set the ylim
    if common_ylim:
        selected_vars = it_matrix.loc[:, [c for c in it_matrix.columns if c[0] == var[0]]].values
        b = (selected_vars.max() - selected_vars.min()) * 0.03
        plt_ylim = (selected_vars.min()-b, selected_vars.max()+b)

    # Plot each year with a different color using the day of the year as x axis
    plt.figure()
    for y in range(it_matrix.index.year.min(), it_matrix.index.year.max() + 1):
        curr = it_matrix[it_matrix.index.year == y]
        curr.set_index(curr.index.dayofyear, inplace=True)
        curr[var].plot(label='_nolegend_')

    # Plot the daily climatology
    clim = it_matrix[var].groupby(it_matrix.index.dayofyear).mean().loc[1:365]
    clim.plot(label='Mean', color='black', linewidth=5)

    # Set the figure properties
    if common_ylim:
        plt.ylim(plt_ylim)
    plt.xlabel('Day')
    plt.ylabel(plt_ylabel[var[0]])
    plt.legend()


def lead_time_rmse(fileName):

    # fileName = '/home/mcallegari@eurac.edu/SECLI-FIRM/Mattia/SF_runoff/Zoccolo/Results/Learning_curve/Runoff_forecast_26_trainingyears.csv'

    # Open the result file
    runoff = pd.read_csv(fileName, index_col=0, parse_dates=True)

    # Create the ensamble mean
    for lt in range(1, 8):
        columns = [c for c in runoff.columns if 'sfTP_m' in c and f'lt{lt}' in c]
        runoff[f'sfTP_em_lt{lt}'] = runoff.loc[:, columns].mean(axis=1)

    # Compute the RMSE
    runoff_error = runoff.apply(lambda y_pred: root_mean_squared_error(runoff['true_runoff'], y_pred), axis=0)

    plt.figure()
    plt.plot([1, 7], [runoff_error['runoff_clim']]*2, color='black', label='runoff climatology')
    plt.plot([1, 7], [runoff_error['trueTP']] * 2, color='red', label='era5')
    plt.plot(range(1, 8), runoff_error[[f'climTP_lt{lt}' for lt in range(1, 8)]], marker='o', label='era5 climatology')
    plt.plot(range(1, 8), runoff_error[[f'sfTP_em_lt{lt}' for lt in range(1, 8)]], marker='o', label='SEAS5 ensamble mean')
    plt.plot(range(1, 8), runoff_error[[f'sfTP_m1_lt{lt}' for lt in range(1, 8)]], color='C1', alpha=0.5, label='SEAS5 members')
    for m in range(2, 26):
        plt.plot(range(1, 8), runoff_error[[f'sfTP_m{m}_lt{lt}' for lt in range(1, 8)]], color='C1', alpha=0.5)
    plt.legend()
    plt.xlabel('Lead time (months)')
    plt.ylabel('RMSE ($m^3/s$)')

    
def spatial_avg_daily_input(daily_input):
    t_columns = [c for c in daily_input.columns if c[0] =='T']
    daily_input['T'] = daily_input[t_columns].mean(axis=1)
    daily_input=daily_input.drop(columns = t_columns)

    e_columns = [c for c in daily_input.columns if c[0] =='E']
    daily_input['E'] = daily_input[e_columns].mean(axis=1)
    daily_input=daily_input.drop(columns = e_columns)

    p_columns = [c for c in daily_input.columns if c[0] =='P']
    daily_input['P'] = daily_input[p_columns].mean(axis=1)
    daily_input=daily_input.drop(columns = p_columns)
    
    return daily_input;


def spatial_stats_daily_input(daily_input):
    
    new_daily_input=pd.DataFrame(daily_input.Q)
    
    t_columns = [c for c in daily_input.columns if c[0] =='T']
    t_vars=daily_input[t_columns]
    new_daily_input.loc[:,'T'] = t_vars.mean(axis=1)
    new_daily_input.loc[:,'T5']=t_vars.quantile(q=0.05,axis=1)
    new_daily_input.loc[:,'T25']=t_vars.quantile(q=0.25,axis=1)
    new_daily_input.loc[:,'T75']=t_vars.quantile(q=0.75,axis=1)
    new_daily_input.insert(loc=5,column='T95',value=t_vars.quantile(q=0.95,axis=1))
    
    e_columns = [c for c in daily_input.columns if c[0] =='E']
    e_vars=daily_input[e_columns]
    new_daily_input.loc[:,'E'] = e_vars.mean(axis=1)
    new_daily_input.loc[:,'E5']=e_vars.quantile(q=0.05,axis=1)
    new_daily_input.loc[:,'E25']=e_vars.quantile(q=0.25,axis=1)
    new_daily_input.loc[:,'E75']=e_vars.quantile(q=0.75,axis=1)
    new_daily_input.loc[:,'E95']=e_vars.quantile(q=0.95,axis=1)
    
    p_columns = [c for c in daily_input.columns if c[0] =='P']
    p_vars=daily_input[p_columns]
    new_daily_input.loc[:,'P'] = p_vars.mean(axis=1)
    new_daily_input.loc[:,'P5']=p_vars.quantile(q=0.05,axis=1)
    new_daily_input.loc[:,'P25']=p_vars.quantile(q=0.25,axis=1)
    new_daily_input.loc[:,'P75']=p_vars.quantile(q=0.75,axis=1)
    new_daily_input.loc[:,'P95']=p_vars.quantile(q=0.95,axis=1)
    
    return new_daily_input



def create_gap(train_index,test_index,gap):
    right=((train_index+1 == test_index[0]).sum()==1) and ((train_index-1 == test_index[-1]).sum()==0)
    centre=((train_index+1 == test_index[0]).sum()==1) and ((train_index-1 == test_index[-1]).sum()==1)
    left = ((train_index+1 == test_index[0]).sum()==0) and ((train_index-1 == test_index[-1]).sum()==1)
    if right:
        train_index=train_index[0:-gap]

    if left:
        train_index=train_index[gap:]

    if centre:
        pos = np.where(train_index+1 == test_index[0])[0][0]
        train_index=np.concatenate((train_index[:pos-gap],train_index[pos+gap:]),axis=0)
        
    return train_index;


def compute_anomalies(climatologies,pred):
    
    #compute real climatology
    anomalies=pd.DataFrame(pred.true_runoff-pred.runoff_clim,columns=['true_runoff'])
    
    #get the climatology of prediction on the wanted days
    clim_on_test_dates = pd.DataFrame(climatologies.loc[anomalies.index.day_of_year])
    
    #create an array with the proper shape
    repeated_clim=np.repeat((np.array(clim_on_test_dates.prediction)[...,np.newaxis]),
                    pred.shape[1]-2,
                    axis=1)

    #subtract clim to predictions
    anomalies_pred= pred.iloc[:,2:]- repeated_clim

    #put together real and predicted anomalies
    anomalies=pd.concat([anomalies,anomalies_pred],axis=1)

    return anomalies