<jupyter_start><jupyter_text>def model_calibration(period_name, period_instance):         use meteo data    _meteo = meteo[                   : period_instance[1]]         runoff data as a target    Qobs   = runoff[period_instance[0]: period_instance[1]]         Define an objective function we want to minimize        def objective_function(params):                 calculate runoff with updated parameters        Qsim = pd.DataFrame({"Qsim": gr4j.simulation(_meteo, params)}, index = _meteo.index)                 Wrap-up simulated runoff with the observed one         Drop all missing values for consistency        Q_df = pd.concat([Qsim, Qobs], axis=1).dropna()                 Calculate objective function value         1 - NS        return 1 - NS(Q_df["Qobs"].values, Q_df["Qsim"].values)         The core of calibration algorithm     Specify objective function     Specify parameters' boundaries within lookup is running     Specify parameters of the algorithm     maxiter - iterations of calibration algorithm run     polish  - local search of the best solution     disp    - displaying progress     seed    - random number seed for reproducibility        result = de(objective_function, gr4j.bounds(), maxiter=50, polish=True, disp=True, seed=142)         return optimal parameters we search        return result.x  Create a holder for optimal model parametersopt_pars = {} Run calibration routine for every period we have definedfor period_name, period_instance in periods.items():        print(period_name)        opt_pars[period_name] = model_calibration(period_name, period_instance)    print("Done!") opt_pars<jupyter_code>opt_pars = {'first' : np.array([2.83722215e+02, 7.05682184e-01, 1.78846738e+01, 2.13065737e+00, 1.69027623e-01, 3.67987870e+00]), 
            'second': np.array([6.66034958e+02, 1.81201518e-01, 8.27679790e+00, 2.04745081e+00, 7.43708427e-01, 5.96813778e+00]), 
            'full'  : np.array([390.87276986,   0.46892242,  13.37185012,   2.06929369, 0.51522738,   4.70385877])}
# Calculate a simulated runoff and place it in the one pandas DataFrame
qsim = pd.concat([pd.DataFrame({"Qsim_{}".format(period) : gr4j.simulation(meteo, opt_pars[period])}, index = meteo.index) for period in opt_pars.keys()], axis=1)
# Create a helper pandas DataFrame for holding all the data together
qdf = pd.concat([qsim, runoff], axis=1)
# Calculate a calibration/validation matrix

## Create a holder
validation_array = np.zeros((3, 3))

## Calibration periods in rows
for i, period_name in enumerate(periods.keys()):
    
    ## Validation periods in columns
    for j, period_instance in enumerate(periods.values()):
        
        Qo = runoff[period_instance[0]: period_instance[1]]
        
        Qs = qdf["Qsim_"+period_name][period_instance[0]: period_instance[1]]
        
        Q_df = pd.concat([Qs, Qo], axis=1).dropna()
        
        validation_array[i, j] = NS(Q_df["Qobs"].values, Q_df["Qsim_"+period_name].values)
        

# Convert to a pandas DataFrame for a nice look
validation_table = pd.DataFrame(validation_array, columns=periods.keys(), index=periods.keys())

# Plot the results
validation_table
# Read in HYDROGRAPH model results
qsim_hydrograph = pd.read_csv("../data/runoff_nelka_hydrograph.csv", header=0, index_col=0, parse_dates=True, names=["Hydrograph"])

qsim_hydrograph = qsim_hydrograph.shift(-1)

# Convert the baseline solution back to m3/s and concatenate with HYdrograph model results
qdf = pd.concat([qdf / coef, qsim_hydrograph], axis=1)

qdf = qdf.dropna()
# Year-wise NS BASELINE
baseline_NS = [NS(qdf["Qobs"][str(year)], qdf["Qsim_full"][str(year)]) for year in range(1980, 1986)]

hydrograph_NS = [NS(qdf["Qobs"][str(year)], qdf["Hydrograph"][str(year)]) for year in range(1980, 1986)]
print("baseline, full period NS:", NS(qdf["Qobs"]["1980":"1985"], qdf["Qsim_full"]["1980":"1985"]))
print("hydrograph, full period NS:", NS(qdf["Qobs"]["1980":"1985"], qdf["Hydrograph"]["1980":"1985"]))
plt.figure(figsize=(12, 8))
plt.bar([i for i in range(1980, 1986)], hydrograph_NS, alpha=0.85, label="HYDROGRAPH")
plt.bar([i for i in range(1980, 1986)], baseline_NS, label="Baseline")
plt.legend()
plt.ylabel("NS")
qdf[["Qobs", "Hydrograph", "Qsim_full"]]["1981"].plot(style=["k-", "r-", "b--"], figsize=(16,10))
qdf[["Qobs", "Hydrograph", "Qsim_full"]]["1985"].plot(style=["k-", "r-", "b--"], figsize=(16,10))<jupyter_output><empty_output><jupyter_text>$$ NSE = 1-\frac{\sum_{i=1}^{N}(obs_{i}-sim_{i})^2}{\sum_{i=1}^{N}(obs_{i}-\bar{obs})^2} $$_______________$$ NSE_{gain} = 1-\frac{\sum_{i=1}^{N}(obs_{i}-sim_{i})^2}{\sum_{i=1}^{N}(obs_{i}-baseline_{i})^2} $$<jupyter_code>def NS_gain(obs, sim, baseline):
    
    numerator = np.nansum( (obs - sim)**2 )
    
    denominator = np.nansum( (obs - baseline)**2 )
    
    return 1 - numerator / denominator
gain_NS = [NS_gain(qdf["Qobs"][str(year)], qdf["Hydrograph"][str(year)], qdf["Qsim_full"][str(year)]) for year in range(1980, 1986)]
plt.figure(figsize=(12, 8))
plt.bar([i for i in range(1980, 1986)], gain_NS, label="Efficiency gain under the baseline")
plt.ylabel("NS gain")
plt.legend()
print("Full period efficiency gain: ", NS_gain(qdf["Qobs"]["1980":"1985"], qdf["Hydrograph"]["1980":"1985"], qdf["Qsim_full"]["1980":"1985"]))<jupyter_output>Full period efficiency gain:  0.4767395780773207
{'fraction_non_alphanumeric': 0.1574492099322799, 'fraction_numerical': 0.07072987208427389, 'mean_word_length': 4.4421699078812695, 'pattern_counts': {'":': 8, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '280130', 'license': 'MIT', 'n_tokens_mistral': 2073, 'n_tokens_neox': 1794, 'n_words': 419, 'path': 'notebooks/nelka.ipynb', 'repo_name': 'hydrogo/2018_VinoRead_wrkshp', 'size': '5316'}