From 5337deec3222f4f009b253ced96ed593ddab66b7 Mon Sep 17 00:00:00 2001 From: ssheeder_rtiintl Date: Thu, 23 May 2024 13:37:58 -0400 Subject: [PATCH] Updated to run Example 4, which includes the CMA-ES optimization routine in the calibrate method of the marrmot_model class --- .vscode/launch.json | 3 +- examples/example_3.py | 89 ++++++++ examples/example_4.py | 198 ++++++++++++++++++ outcmaes/axlen.dat | 2 + outcmaes/axlencorr.dat | 2 + outcmaes/axlenprec.dat | 2 + outcmaes/fit.dat | 2 + outcmaes/sigvec.dat | 2 + outcmaes/stddev.dat | 2 + outcmaes/xmean.dat | 2 + outcmaes/xrecentbest.dat | 2 + .../Functions/Objective_functions/__init__.py | 13 ++ .../Objective_functions/of_inverse_nse.py | 2 +- .../Functions/Objective_functions/of_kge.py | 4 +- .../of_mean_hilo_root5_kge.py | 4 +- .../solver_functions/NewtonRaphson.py | 1 - .../Models/Flux/interflow/interflow_10.py | 30 ++- .../Models/Flux/recharge/recharge_2.py | 4 +- .../Models/Flux/recharge/recharge_3.py | 26 +-- .../Models/Flux/recharge/recharge_5.py | 4 +- .../Models/Flux/recharge/recharge_6.py | 1 + .../Models/Flux/recharge/recharge_7.py | 3 +- .../Models/Flux/recharge/recharge_8.py | 2 + .../Models/Models/m_01_collie1_1p_1s.py | 12 +- .../Models/Models/m_02_wetland_4p_1s.py | 4 +- .../Models/Models/m_03_collie2_4p_1s.py | 3 +- .../Models/Models/m_04_newzealand1_6p_1s.py | 3 +- .../Models/Models/m_05_ihacres_7p_1s.py | 3 +- .../Models/Models/m_06_alpine1_4p_2s.py | 1 + .../Models/Models/m_07_gr4j_4p_2s.py | 3 +- src/pymarrmot/Models/Models/m_08_us1_5p_2s.py | 2 +- .../Models/Models/m_09_susannah1_6p_2s.py | 1 + .../Models/Models/m_10_susannah2_6p_2s.py | 1 + .../Models/Models/m_11_collie3_6p_2s.py | 1 + .../Models/Models/m_12_alpine2_6p_2s.py | 1 + .../Models/Models/m_13_hillslope_7p_2s.py | 1 + .../Models/Models/m_14_topmodel_7p_2s.py | 1 + .../Models/Models/m_15_plateau_8p_2s.py | 1 + .../Models/Models/m_16_newzealand2_8p_2s.py | 1 + .../Models/Models/m_17_penman_4p_3s.py | 1 + .../Models/Models/m_18_simhyd_7p_3s.py | 1 + .../Models/Models/m_19_australia_8p_3s.py | 1 + .../Models/Models/m_20_gsfb_8p_3s.py | 1 + .../Models/Models/m_21_flexb_9p_3s.py | 1 + .../Models/Models/m_28_xinanjiang_12p_4s.py | 1 + .../Models/Models/m_30_mopex2_7p_5s.py | 2 +- .../Models/Models/m_32_mopex4_10p_5s.py | 4 +- .../Models/Models/m_33_sacramento_11p_5s.py | 3 +- .../Models/Models/m_36_modhydrolog_15p_5s.py | 4 +- .../Models/Models/m_37_hbv_15p_5s.py | 4 +- .../Models/Models/m_38_tank2_16p_5s.py | 3 +- .../Models/Models/m_39_mcrm_16p_5s.py | 2 +- .../Models/Models/m_40_smar_8p_6s.py | 2 +- .../Models/Models/m_41_nam_10p_6s.py | 1 + .../Models/Models/m_42_hycymodel_12p_6s.py | 3 +- .../Models/Models/m_44_echo_16p_6s.py | 3 +- .../Models/Models/m_45_prms_18p_7s.py | 10 +- .../Models/Models/m_46_classic_12p_8s.py | 2 +- src/pymarrmot/Models/Models/marrmot_model.py | 55 +++-- src/pymarrmot/Models/Models/temp.py | 18 ++ 60 files changed, 476 insertions(+), 85 deletions(-) create mode 100644 examples/example_3.py create mode 100644 examples/example_4.py create mode 100644 outcmaes/axlen.dat create mode 100644 outcmaes/axlencorr.dat create mode 100644 outcmaes/axlenprec.dat create mode 100644 outcmaes/fit.dat create mode 100644 outcmaes/sigvec.dat create mode 100644 outcmaes/stddev.dat create mode 100644 outcmaes/xmean.dat create mode 100644 outcmaes/xrecentbest.dat create mode 100644 src/pymarrmot/Functions/Objective_functions/__init__.py create mode 100644 src/pymarrmot/Models/Models/temp.py diff --git a/.vscode/launch.json b/.vscode/launch.json index 7774467..889dcfc 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -10,7 +10,8 @@ "type": "debugpy", "request": "launch", "program": "${file}", - "console": "integratedTerminal" + "console": "integratedTerminal", + "justMyCode": true } ] } \ No newline at end of file diff --git a/examples/example_3.py b/examples/example_3.py new file mode 100644 index 0000000..86865c8 --- /dev/null +++ b/examples/example_3.py @@ -0,0 +1,89 @@ +""" +Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter +This file is part of the Modular Assessment of Rainfall-Runoff Models +Toolbox (MARRMoT). +MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY +WARRANTY. See for details. + +Contact: l.trotter@unimelb.edu.au + +This example workflow contains an example application of multiple models +to a single catchment. +It includes 6 steps: + +1. Data preparation +2. Model choice and setup +3. Model solver settings +For each model in list + 4. Model generation and set-up + 5. Model runs +6. Output visualization +""" + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +from pymarrmot.models.models.m_29_hymod_5p_5s import m_29_hymod_5p_5s +from pymarrmot.models.models.m_01_collie1_1p_1s import m_01_collie1_1p_1s +from pymarrmot.models.models.m_27_tank_12p_4s import m_27_tank_12p_4s + +# 1. Prepare data +df = pd.read_csv('c:/users/ssheeder/repos/pymarrmot/examples/Example_DataSet.csv') + +# Create a climatology data input structure +input_climatology = { + 'dates': df['Date'].to_numpy(), # Daily data: date in ??? format + 'precip': df['Precip'].to_numpy(), # Daily data: P rate [mm/d] + 'temp': df['Temp'].to_numpy(), # Daily data: mean T [degree C] + 'pet': df['PET'].to_numpy(), # Daily data: Ep rate [mm/d] +} + +# 2. Define the model settings +model_list = [m_29_hymod_5p_5s, + m_01_collie1_1p_1s, + m_27_tank_12p_4s] + +# 3. Define the solver settings +input_solver_opts = { + 'resnorm_tolerance': 0.1, + 'resnorm_maxiter': 6 +} + +# Create and run all model objects +results_sampling = [] + +results_sampling.append(['model name', 'parameter_values', 'output_ex', 'output_in', 'output_ss', 'output_wb']) + +for model in model_list: + print(f"Now starting model {model}.") + m = model() + m.delta_t = 1 + model_range = m.parRanges + numPar = m.numParams + numStore = m.numStores + + input_theta = model_range[:, 0] + np.random.rand(numPar) * (model_range[:, 1] - model_range[:, 0]) + input_s0 = np.zeros(numStore) + + output_ex, output_in, output_ss, output_waterbalance = m.get_output(4,input_climatology, input_s0, input_theta, input_solver_opts) + + results_sampling.append([model, input_theta, output_ex, output_in, output_ss, output_waterbalance]) + +# 6. Analyze the outputs +t = input_climatology['dates'] +streamflow_observed = df['Q'].to_numpy() + +plt.figure(figsize=(10, 6)) +plt.plot(t, streamflow_observed, 'k', linewidth=2) +for i in range(len(model_list)): + plt.plot(t, results_sampling[i + 1][2]['Q']) + +plt.legend(['Observed'] + model_list, loc='upper right') +plt.title('Model sampling results') +plt.ylabel('Streamflow [mm/d]') +plt.xlabel('Time [d]') +plt.grid(True) +plt.xticks(rotation=45) +plt.ylim([0, 70]) +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/examples/example_4.py b/examples/example_4.py new file mode 100644 index 0000000..5ab3254 --- /dev/null +++ b/examples/example_4.py @@ -0,0 +1,198 @@ +""" +Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter +This file is part of the Modular Assessment of Rainfall-Runoff Models +Toolbox (MARRMoT). +MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY +WARRANTY. See for details. + +Contact: l.trotter@unimelb.edu.au + +This example workflow contains an example application of calibration of a +single models to a single catchment. +It includes 7 steps: + +1. Data preparation +2. Model choice and setup +3. Model solver settings and time-stepping scheme +4. Calibration settings +5. Model calibration +6. Evaluation of calibration results +7. Output visualization + +NOTE: this example uses a custom function 'my_cmaes' to perform the +optimization, it is a wrapper around 'cmaes' to ensure inputs and outputs +are consistent with other MATLAB's optimization algorithms (e.g. +'fminsearch' or 'fminsearchbnd'). +While the wrapper is provided as part of MARRMoT, it requires the source +code to 'cmaes' to function, it is available at: +http://cma.gforge.inria.fr/cmaes.m +The wrapper is necessary for the optimizer to function within the +MARRMoT_model.calibrate method. +Alternatively any model can be calibrated using any optimization +algorithm using the MARRMoT_model.calc_par_fitness method which returns +the value of an objective function and can be used as input to an +optimizer. +""" + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +from datetime import datetime +from pymarrmot.models.models.m_29_hymod_5p_5s import m_29_hymod_5p_5s +import pymarrmot.functions.objective_functions as objective_functions + +# 1. Prepare data +df = pd.read_csv('c:/users/ssheeder/repos/pymarrmot/examples/Example_DataSet.csv') + +# Create a climatology data input structure +input_climatology = { + 'dates': df['Date'].to_numpy(), # Daily data: date in ??? format + 'precip': df['Precip'].to_numpy(), # Daily data: P rate [mm/d] + 'temp': df['Temp'].to_numpy(), # Daily data: mean T [degree C] + 'pet': df['PET'].to_numpy(), # Daily data: Ep rate [mm/d] +} + + + +# Extract observed streamflow +Q_obs = df['Q'].to_numpy() + +# 2. Define the model settings and create the model object +m = m_29_hymod_5p_5s() # Create an instance of the model object +parRanges = m.parRanges # Parameter ranges +numParams = m.numParams # Number of parameters +numStores = m.numStores # Number of stores +input_s0 = np.zeros(numStores) # Initial storages (see note in paragraph 5 on model warm-up) + +# 3. Define the solver settings +input_solver_opts = { + 'resnorm_tolerance': 0.1, # Root-finding convergence tolerance; + # users have reported differences in simulation accuracy (KGE scores) during calibration between Matlab and Octave for a given tolerance. + # In certain cases, Octave seems to require tighter tolerances to obtain the same KGE scores as Matlab does. + 'resnorm_maxiter': 6 # Maximum number of re-runs +} + +# 4. Define calibration settings +# Settings for 'my_cmaes' +# the opts struct is made up of two fields, the names are hardcoded, so +# they cannot be changed: +# .sigma0: initial value of sigma +# .cmaes_opts: struct of options for cmaes, see cmaes documentation +# or type cmaes to see list of options and default values + +# starting sigma +optim_opts = {'insigma': .3 * (parRanges[:, 1] - parRanges[:, 0])} # starting sigma (this is default, could have left it blank) + +# other options +optim_opts['LBounds'] = parRanges[:, 0] # lower bounds of parameters +optim_opts['UBounds'] = parRanges[:, 1] # upper bounds of parameters +optim_opts['PopSize'] = 4 + np.floor(3 * np.log(numParams)) # population size (default) +optim_opts['TolX'] = 1e-6 * np.min(optim_opts['insigma']) # stopping criterion on changes to parameters +optim_opts['TolFun'] = 1e-6 # stopping criterion on changes to fitness function +optim_opts['TolHistFun'] = 1e-6 # stopping criterion on changes to fitness function +optim_opts['SaveFilename'] = 'wf_ex_4_cmaesvars.txt' # output file of cmaes variables +optim_opts['LogFilenamePrefix'] = 'wf_ex_4_' # prefix for cmaes log-files + +# Other useful options +optim_opts['EvalParallel'] = False # change to true to run in parallel on a pool of CPUs (e.g. on a cluster) +# uncomment to restart r-times until the condition +# r = 2 +# optim_opts['Restarts'] = r +# optim_opts['RestartIf'] = 'fmin > -.8' # OF is inverted, so this restarts +# unless the OF (KGE here) > 0.8 + +# debugging options +#optim_opts['MaxIter'] = 5 # just do 5 iterations, to check if it works +#optim_opts['Seed'] = 1234 # for reproducibility + +# initial parameter set +par_ini = np.mean(parRanges, axis=1) # same as default value +m.theta = par_ini # set the initial parameter set + +# Choose the objective function +of_name = 'of_kge' # This function is provided as part of MARRMoT. See ./MARRMoT/Functions/Objective functions +weights = [1, 1, 1] # Weights for the three KGE components + +# Time periods for calibration. +# Indices of timestep to use for calibration, here we are using 1 year for +# warm-up, and 2 years for calibration, leaving the rest of the dataset for +# independent evaluation. +n = len(Q_obs) +warmup = 365 +cal_idx = list(range(warmup + 1, warmup + 365 * 2 + 1)) +eval_idx = list(range(max(cal_idx), n)) + +# 5. Calibrate the model +# MARRMoT model objects have a "calibrate" method that takes uses a chosen +# optimization algorithm and objective function to optimize the parameter +# set. See MARRMoT_model class for details. + +# first set up the model +m.delta_t = 1 +m.input_climate = input_climatology +m.solver_opts = input_solver_opts +m.S0 = input_s0 + +par_opt, of_cal, stopflag, output = m.calibrate( + Q_obs, # observed streamflow + cal_idx, # timesteps to use for model calibration + 'my_cmaes', # function to use for optimization (must have same structure as fminsearch) + par_ini, # initial parameter estimates + optim_opts, # options to optim_fun + of_name, # name of objective function to use + True, True, # should the OF be inversed? Should I display details about the calibration? + weights # additional arguments to of_name +) + +# 6. Evaluate the calibrated parameters on unseen data +# Run the model with calibrated parameters, get only the streamflow +Q_sim = m.get_streamflow(par_opt) + +# Compute evaluation performance +obj_func = getattr(objective_functions, of_name) # Objective function (here 'of_KGE') +of_eval = obj_func(Q_obs, # Observed flow during evaluation period + Q_sim, # Simulated flow during evaluation period, using calibrated parameters + eval_idx, # Indices of evaluation period + weights) # KGE component weights +of_eval = of_eval[0] # Extract the KGE value from the tuple + +# 7. Visualize the results +# Prepare a time vector +t = [datetime.fromtimestamp(t) for t in m.input_climate['dates']] + +# Compare simulated and observed streamflow + +plt.figure(figsize=(10, 6), facecolor='w') +plt.box(True) + +# Flows +plt.plot(t, Q_obs, 'k', label='Q_obs') +plt.plot(t, Q_sim, 'r', label='Q_sim') + +# Dividing line +plt.axvline(t[max(cal_idx)], color='b', linestyle='--', linewidth=2, label='Cal // Eval') +plt.axvline(t[warmup], color='g', linestyle='--', linewidth=2, label='warmup // Cal') + +# Legend & text +plt.legend(loc='upper left') +plt.title('Model calibration and evaluation results') +plt.ylabel('Streamflow [mm/d]') +plt.xlabel('Time [d]') + +max_obs = np.max(Q_obs) +max_sim = np.max(Q_sim) +max_val = 1.05 * max(max_obs, max_sim) +max_t = max(t) +min_t = min(t) +del_t = max(t) - min(t) + +plt_txt = 'Cal ' + of_name + ' = %.2f' % round(of_cal,2) +plt_txt += '\nEval ' + of_name + ' = %.2f' % round(of_eval,2) + +plt.text(min_t + del_t *.85, max_val * .85, plt_txt, fontsize = 10, bbox=dict(facecolor='white', edgecolor='red')) +plt.xticks(rotation=45) + +plt.ylim([0, max_val]) +plt.tick_params(axis='both', which='major', labelsize=12) +plt.tight_layout() +plt.show(block=True) \ No newline at end of file diff --git a/outcmaes/axlen.dat b/outcmaes/axlen.dat new file mode 100644 index 0000000..f264bf3 --- /dev/null +++ b/outcmaes/axlen.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, sigma, max axis length, min axis length, all principal axes lengths (sorted square roots of eigenvalues of C)", seed=581622, Thu May 23 13:25:19 2024 +1 8 0.914647981032298 1.0000400008000108 1.0 1.0 1.00001000005 1.0000200002000015 1.0000300004500045 1.0000400008000108 diff --git a/outcmaes/axlencorr.dat b/outcmaes/axlencorr.dat new file mode 100644 index 0000000..22fef07 --- /dev/null +++ b/outcmaes/axlencorr.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, min max(neg(.)) min(pos(.)) max correlation, correlation matrix principal axes lengths (sorted square roots of eigenvalues of correlation matrix)", seed=581622, Thu May 23 13:25:19 2024 +1 8 -0.10816836668420154 -0.027037173283884528 0.04466017330713045 0.08447965524701906 0.9154832908288328 0.9685323926232423 1.0016151054911655 1.0242153450485072 1.082398011856912 diff --git a/outcmaes/axlenprec.dat b/outcmaes/axlenprec.dat new file mode 100644 index 0000000..5dc8eb6 --- /dev/null +++ b/outcmaes/axlenprec.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, min max(neg(.)) min(pos(.)) max correlation, correlation matrix principal axes lengths (sorted square roots of eigenvalues of correlation matrix)", seed=581622, Thu May 23 13:25:19 2024 +1 8 -0.08454984751796392 -0.04660240315446251 0.0303716355834803 0.10588756675618491 0.9163263717187823 0.970760339100305 0.9957568698812483 1.0254102610826064 1.0839614370785302 diff --git a/outcmaes/fit.dat b/outcmaes/fit.dat new file mode 100644 index 0000000..6c554ec --- /dev/null +++ b/outcmaes/fit.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, sigma, axis ratio, bestever, best, median, worst objective function value, further objective values of best", seed=581622, Thu May 23 13:25:19 2024, {} +1 8 0.914647981032298 1.0000400008000108 -0.03268569325229487 -3.2685693252294867e-02 -0.03268569325229487 -0.03268569325229487 diff --git a/outcmaes/sigvec.dat b/outcmaes/sigvec.dat new file mode 100644 index 0000000..f816ba6 --- /dev/null +++ b/outcmaes/sigvec.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, sigma, beta, void, sigvec==sigma_vec.scaling factors from diagonal decoding", seed=581622, Thu May 23 13:25:19 2024, {} +1 8 0.914647981032298 1 0 599.6999999999999 3.0 0.3 0.3 0.3 diff --git a/outcmaes/stddev.dat b/outcmaes/stddev.dat new file mode 100644 index 0000000..74c2aa5 --- /dev/null +++ b/outcmaes/stddev.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, sigma, void, void, stds==sigma*sigma_vec.scaling*sqrt(diag(C))", seed=581622, Thu May 23 13:25:19 2024, {} +1 8 0.914647981032298 0 0 548.3099360911514 2.805167756156038 0.2614078771696393 0.2624990947176959 0.2713936714567164 diff --git a/outcmaes/xmean.dat b/outcmaes/xmean.dat new file mode 100644 index 0000000..8d84b5a --- /dev/null +++ b/outcmaes/xmean.dat @@ -0,0 +1,2 @@ +% # columns="iteration, evaluation, void, void, void, xmean", seed=581622, Thu May 23 13:25:19 2024, {} # scaling_of_variables: 1, typical_x: 0 +1 8 0 0 0 654.5122932599646 3.661683907932788 0.2777202634101034 0.439071381950522 0.6169564844827764 diff --git a/outcmaes/xrecentbest.dat b/outcmaes/xrecentbest.dat new file mode 100644 index 0000000..5f1c24d --- /dev/null +++ b/outcmaes/xrecentbest.dat @@ -0,0 +1,2 @@ +% # columns="iter, evals, sigma, 0, fitness, xbest" seed=581622, Thu May 23 13:25:19 2024, {} +1 8 0.914647981032298 0 -0.03268569325229487 266.7118468801151 6.042852169155995 0.13597395782092375 0.23869199917925976 0.7612633551792338 diff --git a/src/pymarrmot/Functions/Objective_functions/__init__.py b/src/pymarrmot/Functions/Objective_functions/__init__.py new file mode 100644 index 0000000..3382e1c --- /dev/null +++ b/src/pymarrmot/Functions/Objective_functions/__init__.py @@ -0,0 +1,13 @@ +from pymarrmot.functions.objective_functions.of_bias_penalised_log import of_bias_penalised_log as of_bias_penalised_log +from pymarrmot.functions.objective_functions.of_inverse_kge import of_inverse_kge as of_inverse_kge +from pymarrmot.functions.objective_functions.of_inverse_nse import of_inverse_nse as of_inverse_nse +from pymarrmot.functions.objective_functions.of_kge import of_kge as of_kge +from pymarrmot.functions.objective_functions.of_log_nse import of_log_nse as of_log_nse +from pymarrmot.functions.objective_functions.of_mare import of_mare as of_mare +from pymarrmot.functions.objective_functions.of_mean_hilo_kge import of_mean_hilo_kge as of_mean_hilo_kge +from pymarrmot.functions.objective_functions.of_mean_hilo_root5_kge import of_mean_hilo_root5_kge as of_mean_hilo_root5_kge +from pymarrmot.functions.objective_functions.of_nrmse import of_nrmse as of_nrmse +from pymarrmot.functions.objective_functions.of_nse import of_nse as of_nse +from pymarrmot.functions.objective_functions.of_pcmare import of_pcmare as of_pcmare +from pymarrmot.functions.objective_functions.of_rmse import of_rmse as of_rmse +from pymarrmot.functions.objective_functions.of_root5_kge import of_root5_kge as of_root5_kge \ No newline at end of file diff --git a/src/pymarrmot/Functions/Objective_functions/of_inverse_nse.py b/src/pymarrmot/Functions/Objective_functions/of_inverse_nse.py index 1092d5b..297d430 100644 --- a/src/pymarrmot/Functions/Objective_functions/of_inverse_nse.py +++ b/src/pymarrmot/Functions/Objective_functions/of_inverse_nse.py @@ -2,7 +2,7 @@ from pymarrmot.functions.objective_functions.check_and_select import check_and_select from typing import Tuple -def of_inverse_NSE(obs: np.array, sim: np.array, idx: np.array=None) -> Tuple[float, np.array]: +def of_inverse_nse(obs: np.array, sim: np.array, idx: np.array=None) -> Tuple[float, np.array]: """ Calculates the Nash-Sutcliffe Efficiency (Nash & Sutcliffe, 1970) of the log of simulated streamflow. Ignores time steps with negative flow values. Adds a constant e of 1/100 of mean(obs) to avoid issues with zero diff --git a/src/pymarrmot/Functions/Objective_functions/of_kge.py b/src/pymarrmot/Functions/Objective_functions/of_kge.py index 3349206..bb98448 100644 --- a/src/pymarrmot/Functions/Objective_functions/of_kge.py +++ b/src/pymarrmot/Functions/Objective_functions/of_kge.py @@ -55,8 +55,8 @@ def of_kge(obs: np.array, sim: np.array, idx: np.array=None, w=None): if w is None: w = w_default else: - if not (w.size == 3 and np.ndim(w) == 1): - raise ValueError('Weights should be a 3x1 or 1x3 vector.') + if not (len(w) == 3 and np.ndim(w) == 1): + raise ValueError('Weights should be a 1x3 vector.') # Calculate components c = np.array([np.corrcoef(obs, sim)[0, 1], np.std(sim) / np.std(obs), np.mean(sim) / np.mean(obs)]) diff --git a/src/pymarrmot/Functions/Objective_functions/of_mean_hilo_root5_kge.py b/src/pymarrmot/Functions/Objective_functions/of_mean_hilo_root5_kge.py index c1af037..813b999 100644 --- a/src/pymarrmot/Functions/Objective_functions/of_mean_hilo_root5_kge.py +++ b/src/pymarrmot/Functions/Objective_functions/of_mean_hilo_root5_kge.py @@ -1,10 +1,10 @@ import numpy as np from pymarrmot.functions.objective_functions.check_and_select import check_and_select from pymarrmot.functions.objective_functions.of_kge import of_kge -from pymarrmot.functions.objective_functions.of_mean_hilo_root5_kge import of_root5_kge +from pymarrmot.functions.objective_functions.of_root5_kge import of_root5_kge from typing import Tuple -def of_mean_hilo_root5_KGE(obs: np.array, sim: np.array, idx: np.array=None, w: np.array=None) -> Tuple[float, np.array, np.array, np.array]: +def of_mean_hilo_root5_kge(obs: np.array, sim: np.array, idx: np.array=None, w: np.array=None) -> Tuple[float, np.array, np.array, np.array]: """ Calculates the average Kling-Gupta Efficiency of the untransformed and fifth root of streamflow (Gupta et al., 2009) of the untransformed and fifth root (Chiew et al., 1993) of diff --git a/src/pymarrmot/Functions/solver_functions/NewtonRaphson.py b/src/pymarrmot/Functions/solver_functions/NewtonRaphson.py index fce8681..209620e 100644 --- a/src/pymarrmot/Functions/solver_functions/NewtonRaphson.py +++ b/src/pymarrmot/Functions/solver_functions/NewtonRaphson.py @@ -43,7 +43,6 @@ def NewtonRaphson(fun, x0, options=None): # Initialize x0 = np.array(x0).astype(float) x0 = np.array(x0).reshape(-1, 1) # Ensure column vector - #x0 = np.array(x0).reshape(-1, 1) # Ensure column vector defaultopt = {'TolX': 1e-12, 'TolFun': 1e-6, 'MaxIter': 1000} if options is None: options = defaultopt diff --git a/src/pymarrmot/Models/Flux/interflow/interflow_10.py b/src/pymarrmot/Models/Flux/interflow/interflow_10.py index 4d14883..567f136 100644 --- a/src/pymarrmot/Models/Flux/interflow/interflow_10.py +++ b/src/pymarrmot/Models/Flux/interflow/interflow_10.py @@ -1,22 +1,20 @@ - def interflow_10(S, p1, p2, p3): -%interflow_10 -% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter -% This file is part of the Modular Assessment of Rainfall-Runoff Models -% Toolbox (MARRMoT). -% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY -% WARRANTY. See for details. +# Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter +# This file is part of the Modular Assessment of Rainfall-Runoff Models +# Toolbox (MARRMoT). +# MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY +# WARRANTY. See for details. -% Flux function -% ------------------ -% Description: Scaled linear interflow if storage exceeds a threshold -% Constraints: f = 0, for S < p2 -% @(Inputs): p1 - time coefficient [d-1] -% p2 - threshold for flow generation [mm] -% p3 - linear scaling parameter [-] -% S - current storage [mm] +# Flux function +# ------------------ +# Description: Scaled linear interflow if storage exceeds a threshold +# Constraints: f = 0, for S < p2 +# @(Inputs): p1 - time coefficient [d-1] +# p2 - threshold for flow generation [mm] +# p3 - linear scaling parameter [-] +# S - current storage [mm]''' -out = p1*max(0,S-p2)/(p3); + out = p1*max(0,S-p2)/(p3) return out diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_2.py b/src/pymarrmot/Models/Flux/recharge/recharge_2.py index 68b7b13..31f24e2 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_2.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_2.py @@ -8,4 +8,6 @@ def recharge_2(p1, S, Smax, flux): # Smax - maximum contributing storage [mm] # flux - incoming flux [mm/d] - return flux * ((max(S, 0) / Smax) ** p1) + out = flux * ((max(S, 0) / Smax) ** p1) + + return out diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_3.py b/src/pymarrmot/Models/Flux/recharge/recharge_3.py index 2741b44..ecefa35 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_3.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_3.py @@ -1,20 +1,16 @@ def recharge_3(p1, S): -%recharge_3 + # Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter + # This file is part of the Modular Assessment of Rainfall-Runoff Models + # Toolbox (MARRMoT). + # MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY + # WARRANTY. See for details. -% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter -% This file is part of the Modular Assessment of Rainfall-Runoff Models -% Toolbox (MARRMoT). -% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY -% WARRANTY. See for details. - -% Flux function -% ------------------ -# Description: Linear recharge -# Constraints: - -# Inputs: p1 - time coefficient [d-1] -% S - current storage [mm] + # Flux function + # ------------------ + # Description: Linear recharge + # Constraints: - + # Inputs: p1 - time coefficient [d-1] + # S - current storage [mm] return p1 * S -end - diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_5.py b/src/pymarrmot/Models/Flux/recharge/recharge_5.py index 22d00c2..cfcc0a3 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_5.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_5.py @@ -8,4 +8,6 @@ def recharge_5(p1, p2, S1, S2): # S1 - current storage in S1 [mm] # S2 - current storage in S1 [mm] - return p1 * S1 * (1 - min(1, S2 / p2)) + out = p1 * S1 * (1 - min(1, S2 / p2)) + + return out \ No newline at end of file diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_6.py b/src/pymarrmot/Models/Flux/recharge/recharge_6.py index 056e195..b0ca981 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_6.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_6.py @@ -11,4 +11,5 @@ def recharge_6(p1, p2, S, dt): # dt - time step size [d] out = min(max(S/dt, 0), p1 * max(S, 0)**p2) + return out diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_7.py b/src/pymarrmot/Models/Flux/recharge/recharge_7.py index ca49089..b7770ac 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_7.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_7.py @@ -7,4 +7,5 @@ def recharge_7(p1, fin): :param fin: incoming flux [mm/d] :return: minimum of p1 and fin ''' - return min(p1, fin) + out = min(p1, fin) + return out diff --git a/src/pymarrmot/Models/Flux/recharge/recharge_8.py b/src/pymarrmot/Models/Flux/recharge/recharge_8.py index c66ce24..69d1706 100644 --- a/src/pymarrmot/Models/Flux/recharge/recharge_8.py +++ b/src/pymarrmot/Models/Flux/recharge/recharge_8.py @@ -9,5 +9,7 @@ def recharge_8(p1, S, Smax, p2, dt): # Smax - maximum contributing storage [mm] # p2 - maximum flux rate [mm/d] # dt - time step size [d] + out = min(p2*((max(S,0)/Smax)**p1),max(S/dt,0)) + return out diff --git a/src/pymarrmot/Models/Models/m_01_collie1_1p_1s.py b/src/pymarrmot/Models/Models/m_01_collie1_1p_1s.py index 00725f5..99bbb1e 100644 --- a/src/pymarrmot/Models/Models/m_01_collie1_1p_1s.py +++ b/src/pymarrmot/Models/Models/m_01_collie1_1p_1s.py @@ -1,8 +1,9 @@ +import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.evaporation import evap_7 from pymarrmot.models.flux.saturation import saturation_1 -class M_01_Collie1_1p_1s(MARRMoT_model): +class m_01_collie1_1p_1s(MARRMoT_model): """ Class for hydrologic conceptual model: Collie River 1 (traditional bucket model) Model reference: @@ -12,11 +13,16 @@ class M_01_Collie1_1p_1s(MARRMoT_model): 174–198. doi: 10.1016/S0022-1694(01)497-6. """ def __init__(self): + super().__init__() self.numStores = 1 # number of model stores self.numFluxes = 2 # number of model fluxes self.numParams = 1 self.JacobPattern = [1] # Jacobian matrix of model store ODEs - self.parRanges = [1, 2000] # Smax [mm] + #self.parRanges = [1, 2000] # Smax [mm] + + self.parRanges = np.array([[1, 2000]]) # Smax, Maximum soil moisture storage [mm] + + self.StoreNames = ["S1"] # Names for the stores self.FluxNames = ["ea", "qse"] # Names for the fluxes self.FluxGroups = {'Ea': 1, 'Q': 2} # Index or indices of fluxes to add to Actual ET and Streamflow @@ -47,7 +53,7 @@ def model_fun(self, S): dS1 = P - flux_ea - flux_qse # outputs - dS = [dS1] + dS = np.array([dS1]) fluxes = [flux_ea, flux_qse] return dS, fluxes diff --git a/src/pymarrmot/Models/Models/m_02_wetland_4p_1s.py b/src/pymarrmot/Models/Models/m_02_wetland_4p_1s.py index 5daca09..ebec95c 100644 --- a/src/pymarrmot/Models/Models/m_02_wetland_4p_1s.py +++ b/src/pymarrmot/Models/Models/m_02_wetland_4p_1s.py @@ -1,3 +1,4 @@ +import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.interception import interception_2 from pymarrmot.models.flux.evaporation import evap_1 @@ -11,6 +12,7 @@ class M_02_Wetland_4p_1s(MARRMoT_model): Savenije, H. H. G. (2010). “Topography driven conceptual modelling (FLEX-Topo).” Hydrology and Earth System Sciences, 14(12), 2681–2692. https://doi.org/10.5194/hess-14-2681-2010 """ def __init__(self): + super().__init__() self.num_stores = 1 # number of model stores self.num_fluxes = 5 # number of model fluxes self.num_params = 4 @@ -52,7 +54,7 @@ def model_fun(self, S): # stores ODEs dS1 = flux_pe - flux_ew - flux_qwsof - flux_qwgw # outputs - dS = [dS1] + dS = np.array([dS1]) fluxes = [flux_pe, flux_ei, flux_ew, flux_qwsof, flux_qwgw] return dS, fluxes diff --git a/src/pymarrmot/Models/Models/m_03_collie2_4p_1s.py b/src/pymarrmot/Models/Models/m_03_collie2_4p_1s.py index a103c65..502030a 100644 --- a/src/pymarrmot/Models/Models/m_03_collie2_4p_1s.py +++ b/src/pymarrmot/Models/Models/m_03_collie2_4p_1s.py @@ -1,3 +1,4 @@ +import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.evaporation import (evap_7, evap_3) from pymarrmot.models.flux.saturation import saturation_1 @@ -85,7 +86,7 @@ def model_fun(self, S): dS1 = P - flux_eb - flux_ev - flux_qse - flux_qss - dS = [dS1] + dS = np.array([dS1]) fluxes = [flux_eb, flux_ev, flux_qse, flux_qss] return dS, fluxes diff --git a/src/pymarrmot/Models/Models/m_04_newzealand1_6p_1s.py b/src/pymarrmot/Models/Models/m_04_newzealand1_6p_1s.py index 435cc6e..eabba41 100644 --- a/src/pymarrmot/Models/Models/m_04_newzealand1_6p_1s.py +++ b/src/pymarrmot/Models/Models/m_04_newzealand1_6p_1s.py @@ -1,3 +1,4 @@ +import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.evaporation import (evap_6, evap_5) from pymarrmot.models.flux.saturation import saturation_1 @@ -95,7 +96,7 @@ def model_fun(self, S): dS1 = P - flux_veg - flux_ebs - flux_qse - flux_qss - flux_qbf # outputs - dS = [dS1] + dS = np.array([dS1]) fluxes = [flux_veg, flux_ebs, flux_qse, flux_qss, flux_qbf] return dS, fluxes diff --git a/src/pymarrmot/Models/Models/m_05_ihacres_7p_1s.py b/src/pymarrmot/Models/Models/m_05_ihacres_7p_1s.py index e45c5e1..241a516 100644 --- a/src/pymarrmot/Models/Models/m_05_ihacres_7p_1s.py +++ b/src/pymarrmot/Models/Models/m_05_ihacres_7p_1s.py @@ -1,3 +1,4 @@ +import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.evaporation import evap_12 from pymarrmot.models.flux.split import split_1 @@ -121,7 +122,7 @@ def model_fun(self, S): dS1 = -P + flux_ea + flux_u # outputs - dS = dS1 + dS = np.array([dS1]) fluxes = [flux_ea, flux_u, flux_uq, flux_us, flux_xq, flux_xs, flux_xt] return dS, fluxes diff --git a/src/pymarrmot/Models/Models/m_06_alpine1_4p_2s.py b/src/pymarrmot/Models/Models/m_06_alpine1_4p_2s.py index 1a3de78..d4012fe 100644 --- a/src/pymarrmot/Models/Models/m_06_alpine1_4p_2s.py +++ b/src/pymarrmot/Models/Models/m_06_alpine1_4p_2s.py @@ -17,6 +17,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 6 # number of model fluxes self.numParams = 4 diff --git a/src/pymarrmot/Models/Models/m_07_gr4j_4p_2s.py b/src/pymarrmot/Models/Models/m_07_gr4j_4p_2s.py index 1f927cb..b9ceec5 100644 --- a/src/pymarrmot/Models/Models/m_07_gr4j_4p_2s.py +++ b/src/pymarrmot/Models/Models/m_07_gr4j_4p_2s.py @@ -38,6 +38,7 @@ def init(self): """ INITialisation function """ + super().__init__() theta = self.theta delta_t = self.delta_t x1 = theta[0] # Maximum soil moisture storage [mm] @@ -102,7 +103,7 @@ def model_fun(self, S): dS2 = flux_q9 + flux_fr - flux_qr # outputs - dS = [dS1, dS2] + dS = np.array([dS1, dS2]) fluxes = [flux_pn, flux_en, flux_ef, flux_ps, flux_es, flux_perc, flux_q9, flux_q1, flux_fr, flux_fq, flux_qr, flux_qt, flux_ex] diff --git a/src/pymarrmot/Models/Models/m_08_us1_5p_2s.py b/src/pymarrmot/Models/Models/m_08_us1_5p_2s.py index 768a483..f42b734 100644 --- a/src/pymarrmot/Models/Models/m_08_us1_5p_2s.py +++ b/src/pymarrmot/Models/Models/m_08_us1_5p_2s.py @@ -68,7 +68,7 @@ def model_fun(self, S): dS2 = flux_rg + flux_se - flux_esatveg - flux_esatbs - flux_qse - flux_qss # outputs - dS = [dS1, dS2] + dS = np.array([dS1, dS2]) fluxes = [flux_eusei, flux_eusveg, flux_eusbs, flux_esatveg, flux_esatbs, flux_rg, flux_se, flux_qse, flux_qss] diff --git a/src/pymarrmot/Models/Models/m_09_susannah1_6p_2s.py b/src/pymarrmot/Models/Models/m_09_susannah1_6p_2s.py index 1c74a55..cdf8d49 100644 --- a/src/pymarrmot/Models/Models/m_09_susannah1_6p_2s.py +++ b/src/pymarrmot/Models/Models/m_09_susannah1_6p_2s.py @@ -27,6 +27,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 7 # number of model fluxes self.numParams = 6 diff --git a/src/pymarrmot/Models/Models/m_10_susannah2_6p_2s.py b/src/pymarrmot/Models/Models/m_10_susannah2_6p_2s.py index 6afb64d..77ceb97 100644 --- a/src/pymarrmot/Models/Models/m_10_susannah2_6p_2s.py +++ b/src/pymarrmot/Models/Models/m_10_susannah2_6p_2s.py @@ -27,6 +27,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 8 # number of model fluxes self.numParams = 6 diff --git a/src/pymarrmot/Models/Models/m_11_collie3_6p_2s.py b/src/pymarrmot/Models/Models/m_11_collie3_6p_2s.py index 14637a8..018a1ba 100644 --- a/src/pymarrmot/Models/Models/m_11_collie3_6p_2s.py +++ b/src/pymarrmot/Models/Models/m_11_collie3_6p_2s.py @@ -29,6 +29,7 @@ def __init__(self, delta_t=None, theta=None): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 7 # number of model fluxes self.numParams = 6 diff --git a/src/pymarrmot/Models/Models/m_12_alpine2_6p_2s.py b/src/pymarrmot/Models/Models/m_12_alpine2_6p_2s.py index a1ba2a0..3389927 100644 --- a/src/pymarrmot/Models/Models/m_12_alpine2_6p_2s.py +++ b/src/pymarrmot/Models/Models/m_12_alpine2_6p_2s.py @@ -30,6 +30,7 @@ def __init__(self, delta_t=None, theta=None): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 7 # number of model fluxes self.numParams = 6 diff --git a/src/pymarrmot/Models/Models/m_13_hillslope_7p_2s.py b/src/pymarrmot/Models/Models/m_13_hillslope_7p_2s.py index c94e2ad..ada24ab 100644 --- a/src/pymarrmot/Models/Models/m_13_hillslope_7p_2s.py +++ b/src/pymarrmot/Models/Models/m_13_hillslope_7p_2s.py @@ -29,6 +29,7 @@ def __init__(self, delta_t=None, theta=None): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 10 # number of model fluxes self.numParams = 7 diff --git a/src/pymarrmot/Models/Models/m_14_topmodel_7p_2s.py b/src/pymarrmot/Models/Models/m_14_topmodel_7p_2s.py index 678fccc..ed03453 100644 --- a/src/pymarrmot/Models/Models/m_14_topmodel_7p_2s.py +++ b/src/pymarrmot/Models/Models/m_14_topmodel_7p_2s.py @@ -41,6 +41,7 @@ def __init__(self, delta_t=None, theta=None): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 6 # number of model fluxes self.numParams = 7 diff --git a/src/pymarrmot/Models/Models/m_15_plateau_8p_2s.py b/src/pymarrmot/Models/Models/m_15_plateau_8p_2s.py index 245f5c3..b810be4 100644 --- a/src/pymarrmot/Models/Models/m_15_plateau_8p_2s.py +++ b/src/pymarrmot/Models/Models/m_15_plateau_8p_2s.py @@ -56,6 +56,7 @@ def init(self): """ INITialisation function """ + super().__init__() # parameters theta = self.theta delta_t = self.delta_t diff --git a/src/pymarrmot/Models/Models/m_16_newzealand2_8p_2s.py b/src/pymarrmot/Models/Models/m_16_newzealand2_8p_2s.py index 84cd1d4..9bb0d86 100644 --- a/src/pymarrmot/Models/Models/m_16_newzealand2_8p_2s.py +++ b/src/pymarrmot/Models/Models/m_16_newzealand2_8p_2s.py @@ -29,6 +29,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 2 # number of model stores self.numFluxes = 8 # number of model fluxes self.numParams = 8 diff --git a/src/pymarrmot/Models/Models/m_17_penman_4p_3s.py b/src/pymarrmot/Models/Models/m_17_penman_4p_3s.py index 0b6294d..f13d0cc 100644 --- a/src/pymarrmot/Models/Models/m_17_penman_4p_3s.py +++ b/src/pymarrmot/Models/Models/m_17_penman_4p_3s.py @@ -31,6 +31,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 3 # number of model stores self.numFluxes = 7 # number of model fluxes self.numParams = 4 diff --git a/src/pymarrmot/Models/Models/m_18_simhyd_7p_3s.py b/src/pymarrmot/Models/Models/m_18_simhyd_7p_3s.py index 419a90f..e5850e3 100644 --- a/src/pymarrmot/Models/Models/m_18_simhyd_7p_3s.py +++ b/src/pymarrmot/Models/Models/m_18_simhyd_7p_3s.py @@ -30,6 +30,7 @@ def __init__(self): """ creator method """ + super().__init__() self.numStores = 3 # number of model stores self.numFluxes = 10 # number of model fluxes self.numParams = 7 # number of model parameters diff --git a/src/pymarrmot/Models/Models/m_19_australia_8p_3s.py b/src/pymarrmot/Models/Models/m_19_australia_8p_3s.py index c096198..c9e5122 100644 --- a/src/pymarrmot/Models/Models/m_19_australia_8p_3s.py +++ b/src/pymarrmot/Models/Models/m_19_australia_8p_3s.py @@ -23,6 +23,7 @@ def __init__(self): """ Creator method """ + super().__init__() self.numStores = 3 # number of model stores self.numFluxes = 8 # number of model fluxes self.numParams = 8 # number of model parameters diff --git a/src/pymarrmot/Models/Models/m_20_gsfb_8p_3s.py b/src/pymarrmot/Models/Models/m_20_gsfb_8p_3s.py index 3bfa807..7222ca4 100644 --- a/src/pymarrmot/Models/Models/m_20_gsfb_8p_3s.py +++ b/src/pymarrmot/Models/Models/m_20_gsfb_8p_3s.py @@ -24,6 +24,7 @@ def __init__(self): """ Creator method """ + super().__init__() self.numStores = 3 # number of model stores self.numFluxes = 6 # number of model fluxes self.numParams = 8 # number of model parameters diff --git a/src/pymarrmot/Models/Models/m_21_flexb_9p_3s.py b/src/pymarrmot/Models/Models/m_21_flexb_9p_3s.py index 8d176d3..5bcd455 100644 --- a/src/pymarrmot/Models/Models/m_21_flexb_9p_3s.py +++ b/src/pymarrmot/Models/Models/m_21_flexb_9p_3s.py @@ -23,6 +23,7 @@ def __init__(self): """ Creator method """ + super().__init__() self.numStores = 3 # number of model stores self.numFluxes = 9 # number of model fluxes self.numParams = 9 # number of model parameters diff --git a/src/pymarrmot/Models/Models/m_28_xinanjiang_12p_4s.py b/src/pymarrmot/Models/Models/m_28_xinanjiang_12p_4s.py index fd496df..65e66e1 100644 --- a/src/pymarrmot/Models/Models/m_28_xinanjiang_12p_4s.py +++ b/src/pymarrmot/Models/Models/m_28_xinanjiang_12p_4s.py @@ -17,6 +17,7 @@ def __init__(self): """ Creator method. """ + super().__init__() self.numStores = 4 # number of model stores self.numFluxes = 10 # number of model fluxes self.numParams = 12 diff --git a/src/pymarrmot/Models/Models/m_30_mopex2_7p_5s.py b/src/pymarrmot/Models/Models/m_30_mopex2_7p_5s.py index 0cd651d..6d50b5c 100644 --- a/src/pymarrmot/Models/Models/m_30_mopex2_7p_5s.py +++ b/src/pymarrmot/Models/Models/m_30_mopex2_7p_5s.py @@ -100,7 +100,7 @@ def model_fun(self, S): dS5 = flux_q2u - flux_qs # outputs - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_ps, flux_pr, flux_qn, flux_et1, flux_q1f, flux_qw, flux_et2, flux_q2u, flux_qf, flux_qs] diff --git a/src/pymarrmot/Models/Models/m_32_mopex4_10p_5s.py b/src/pymarrmot/Models/Models/m_32_mopex4_10p_5s.py index c7287a6..0a1e83f 100644 --- a/src/pymarrmot/Models/Models/m_32_mopex4_10p_5s.py +++ b/src/pymarrmot/Models/Models/m_32_mopex4_10p_5s.py @@ -1,3 +1,5 @@ +import numpy as np + from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.snowfall import snowfall_1 from pymarrmot.models.flux.rainfall import rainfall_1 @@ -112,7 +114,7 @@ def model_fun(self, S): dS5 = flux_q2u - flux_qs # Outputs - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_ps, flux_pr, flux_qn, flux_et1, flux_i, flux_q1f, flux_qw, flux_et2, flux_q2f, flux_q2u, flux_qf, flux_qs] diff --git a/src/pymarrmot/Models/Models/m_33_sacramento_11p_5s.py b/src/pymarrmot/Models/Models/m_33_sacramento_11p_5s.py index d9c699b..f6de8d7 100644 --- a/src/pymarrmot/Models/Models/m_33_sacramento_11p_5s.py +++ b/src/pymarrmot/Models/Models/m_33_sacramento_11p_5s.py @@ -32,6 +32,7 @@ class m_33_sacramento_11p_5s(MARRMoT_model): Beach, CA, (1), 103�106.""" def __init__(self): + super().__init__() self.theta_derived = None self.numStores = 5 @@ -139,7 +140,7 @@ def model_fun(self, S: list[float]) -> tuple[list[float], list[float]]: dS4 = flux_twexlp + flux_pcfwp - flux_rlp - flux_qbfp dS5 = flux_twexls + flux_pcfws - flux_rls - flux_qbfs - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_qdir, flux_peff, flux_ru, flux_euztw, flux_twexu, flux_qsur, flux_qint, flux_euzfw, flux_pc, flux_pctw, flux_elztw, flux_twexl, flux_twexlp, flux_twexls, flux_pcfwp, diff --git a/src/pymarrmot/Models/Models/m_36_modhydrolog_15p_5s.py b/src/pymarrmot/Models/Models/m_36_modhydrolog_15p_5s.py index 845ad98..7f91123 100644 --- a/src/pymarrmot/Models/Models/m_36_modhydrolog_15p_5s.py +++ b/src/pymarrmot/Models/Models/m_36_modhydrolog_15p_5s.py @@ -1,3 +1,5 @@ +import numpy as np + from pymarrmot.models.models.marrmot_model import MARRMoT_model from pymarrmot.models.flux.evaporation import evap_1, evap_2 from pymarrmot.models.flux.interception import interception_1 @@ -146,7 +148,7 @@ def model_fun(self, S): dS5 = flux_SRUN + flux_INT + flux_FLOW - flux_Q # outputs - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_Ei, flux_EXC, flux_INF, flux_INT, flux_REC, flux_SMF, flux_Et, flux_GWF, flux_TRAP, flux_Ed, flux_DINF, flux_SEEP, diff --git a/src/pymarrmot/Models/Models/m_37_hbv_15p_5s.py b/src/pymarrmot/Models/Models/m_37_hbv_15p_5s.py index 0e6daac..bcdfdfb 100644 --- a/src/pymarrmot/Models/Models/m_37_hbv_15p_5s.py +++ b/src/pymarrmot/Models/Models/m_37_hbv_15p_5s.py @@ -1,7 +1,8 @@ import numpy as np from pymarrmot.models.models.marrmot_model import MARRMoT_model -from models.unit_hydro import(uh_4_full, update_uh) +from models.unit_hydro.uh_4_full import(uh_4_full) +from models.unit_hydro.update_uh import(update_uh) class m_37_hbv_15p_5s(MARRMoT_model): """ @@ -24,6 +25,7 @@ def __init__(self): """ Creator method """ + super().__init__() self.num_stores = 5 # number of model stores self.num_fluxes = 13 # number of model fluxes self.num_params = 15 diff --git a/src/pymarrmot/Models/Models/m_38_tank2_16p_5s.py b/src/pymarrmot/Models/Models/m_38_tank2_16p_5s.py index 4cd6625..c779b2a 100644 --- a/src/pymarrmot/Models/Models/m_38_tank2_16p_5s.py +++ b/src/pymarrmot/Models/Models/m_38_tank2_16p_5s.py @@ -21,6 +21,7 @@ def __init__(self): """ creator method """ + super().__init__() self.num_stores = 5 # number of model stores self.num_fluxes = 14 # number of model fluxes self.num_params = 16 @@ -157,7 +158,7 @@ def model_fun(self, S): dS5 = flux_t2 # outputs - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_t1, flux_t2, flux_y1, flux_y2, flux_y3, flux_y4, flux_y5, flux_e1, flux_e2, flux_e3, flux_e4, diff --git a/src/pymarrmot/Models/Models/m_39_mcrm_16p_5s.py b/src/pymarrmot/Models/Models/m_39_mcrm_16p_5s.py index f95dc27..d110fe7 100644 --- a/src/pymarrmot/Models/Models/m_39_mcrm_16p_5s.py +++ b/src/pymarrmot/Models/Models/m_39_mcrm_16p_5s.py @@ -140,7 +140,7 @@ def model_fun(self, S): dS4 = flux_uib - flux_uob - flux_qic dS5 = flux_uob - flux_qoc - dS = [dS1, dS2, dS3, dS4, dS5] + dS = np.array([dS1, dS2, dS3, dS4, dS5]) fluxes = [flux_ec, flux_qt, flux_qr, flux_er, flux_qn, flux_qd, flux_qp, flux_qb, flux_uib, flux_uob, flux_qic, flux_qoc] diff --git a/src/pymarrmot/Models/Models/m_40_smar_8p_6s.py b/src/pymarrmot/Models/Models/m_40_smar_8p_6s.py index 9cdaadd..14c71bb 100644 --- a/src/pymarrmot/Models/Models/m_40_smar_8p_6s.py +++ b/src/pymarrmot/Models/Models/m_40_smar_8p_6s.py @@ -134,7 +134,7 @@ def model_fun(self, S): dS6 = flux_rg - flux_qg # outputs - dS = [dS1, dS2, dS3, dS4, dS5, dS6] + dS = np.arry([dS1, dS2, dS3, dS4, dS5, dS6]) fluxes = [flux_pstar, flux_estar, flux_evap, flux_r1, flux_i, flux_r2, flux_e1, flux_e2, flux_e3, flux_e4, flux_e5, flux_q1, flux_q2, flux_q3, flux_q4, diff --git a/src/pymarrmot/Models/Models/m_41_nam_10p_6s.py b/src/pymarrmot/Models/Models/m_41_nam_10p_6s.py index 58296d0..1bf8cec 100644 --- a/src/pymarrmot/Models/Models/m_41_nam_10p_6s.py +++ b/src/pymarrmot/Models/Models/m_41_nam_10p_6s.py @@ -31,6 +31,7 @@ def __init__(self): """ Creator method """ + super().__init__() self.numStores = 6 # number of model stores self.numFluxes = 14 # number of model fluxes self.numParams = 10 diff --git a/src/pymarrmot/Models/Models/m_42_hycymodel_12p_6s.py b/src/pymarrmot/Models/Models/m_42_hycymodel_12p_6s.py index f6c1efb..77711ee 100644 --- a/src/pymarrmot/Models/Models/m_42_hycymodel_12p_6s.py +++ b/src/pymarrmot/Models/Models/m_42_hycymodel_12p_6s.py @@ -29,6 +29,7 @@ def __init__(self): """ Creator method. """ + super().__init__() self.numStores = 6 # number of model stores self.numFluxes = 18 # number of model fluxes self.numParams = 12 @@ -152,7 +153,7 @@ def model_fun(self, S): dS6 = flux_rc - flux_qc # outputs - dS = [dS1, dS2, dS3, dS4, dS5, dS6] + dS = np.array([dS1, dS2, dS3, dS4, dS5, dS6]) fluxes = [flux_rc, flux_rg, flux_eic, flux_qie, flux_qis, flux_rt, flux_eis, flux_rs, flux_rn, flux_esu, flux_re, flux_qin, flux_esb, flux_qb, flux_qh, flux_qc, flux_ec, flux_qt] diff --git a/src/pymarrmot/Models/Models/m_44_echo_16p_6s.py b/src/pymarrmot/Models/Models/m_44_echo_16p_6s.py index f696c9c..9ca12fe 100644 --- a/src/pymarrmot/Models/Models/m_44_echo_16p_6s.py +++ b/src/pymarrmot/Models/Models/m_44_echo_16p_6s.py @@ -39,6 +39,7 @@ def __init__(self): """ creator method """ + super().__init__() self.aux_theta = None self.numStores = 6 @@ -150,7 +151,7 @@ def model_fun(self, S): dS4 = flux_fi - flux_et - flux_rd - flux_l dS5 = flux_lf - flux_rf dS6 = flux_ls - flux_rs - dS = [dS1, dS2, dS3, dS4, dS5, dS6] + dS = np.array([dS1, dS2, dS3, dS4, dS5, dS6]) fluxes = [flux_ei, flux_pn, flux_ps, flux_pr, flux_ms, flux_fs, flux_gs, flux_mw, flux_ew, flux_eq, flux_rh, flux_eps, flux_et, flux_fi, flux_rd, diff --git a/src/pymarrmot/Models/Models/m_45_prms_18p_7s.py b/src/pymarrmot/Models/Models/m_45_prms_18p_7s.py index 0cb3595..288c03d 100644 --- a/src/pymarrmot/Models/Models/m_45_prms_18p_7s.py +++ b/src/pymarrmot/Models/Models/m_45_prms_18p_7s.py @@ -174,11 +174,11 @@ def model_fun(self, S): dS7 = flux_sep + flux_gad - flux_bas - flux_snk # outputs dS = np.array([dS1, dS2, dS3, dS4, dS5, dS6, dS7]) - fluxes = np.array([flux_ps, flux_pr, flux_pim, flux_psm, flux_pby, - flux_pin, flux_ptf, flux_m, flux_mim, flux_msm, - flux_sas, flux_sro, flux_inf, flux_pc, flux_excs, - flux_qres, flux_sep, flux_gad, flux_ras, flux_bas, - flux_snk, flux_ein, flux_eim, flux_ea, flux_et]) + fluxes = [flux_ps, flux_pr, flux_pim, flux_psm, flux_pby, + flux_pin, flux_ptf, flux_m, flux_mim, flux_msm, + flux_sas, flux_sro, flux_inf, flux_pc, flux_excs, + flux_qres, flux_sep, flux_gad, flux_ras, flux_bas, + flux_snk, flux_ein, flux_eim, flux_ea, flux_et] return dS, fluxes def step(self): diff --git a/src/pymarrmot/Models/Models/m_46_classic_12p_8s.py b/src/pymarrmot/Models/Models/m_46_classic_12p_8s.py index 33d3c33..27da55e 100644 --- a/src/pymarrmot/Models/Models/m_46_classic_12p_8s.py +++ b/src/pymarrmot/Models/Models/m_46_classic_12p_8s.py @@ -163,7 +163,7 @@ def model_fun(self, S): dS8 = flux_pie - flux_u # Outputs - dS = [dS1, dS2, dS3, dS4, dS5, dS6, dS7, dS8] + dS = np.array([dS1, dS2, dS3, dS4, dS5, dS6, dS7, dS8]) fluxes = [flux_pp, flux_ps, flux_pi, flux_epx, flux_ppx, flux_epy, flux_ppe, flux_q, flux_psd, flux_psi, flux_esx, flux_psx, flux_esy, flux_pse, flux_psq, diff --git a/src/pymarrmot/Models/Models/marrmot_model.py b/src/pymarrmot/Models/Models/marrmot_model.py index 17abada..1eddcae 100644 --- a/src/pymarrmot/Models/Models/marrmot_model.py +++ b/src/pymarrmot/Models/Models/marrmot_model.py @@ -1,13 +1,13 @@ import numpy as np import numbers -from scipy import optimize -from scipy.sparse import csr_matrix + from scipy.optimize import fsolve, least_squares -from scipy.sparse.linalg import spsolve -from scipy.linalg import solve +import pymarrmot.functions.objective_functions as objective_functions from pymarrmot.functions.solver_functions import NewtonRaphson as nr +import cma as cma + class MARRMoT_model: """ Superclass for all MARRMoT models @@ -484,7 +484,9 @@ def get_streamflow(obj, *args): if args or not obj.status or obj.status == 0: obj.run(*args) - Q = np.sum(obj.fluxes[:, obj.FluxGroups['Q']], axis=1) + q = [x - 1 for x in obj.FluxGroups['Q']] + Q = np.sum(obj.fluxes[:,q], axis=1) + #Q = np.sum(obj.fluxes[:, obj.FluxGroups['Q']], axis=1) return Q def calibrate(obj, Q_obs, cal_idx, optim_fun, par_ini=None, optim_opts=None, of_name=None, @@ -527,11 +529,15 @@ def calibrate(obj, Q_obs, cal_idx, optim_fun, par_ini=None, optim_opts=None, of_ output : dict Output, see the documentation of the optimization function for detail. """ - if (not obj.input_climate) or (not obj.delta_t) or (not obj.S0) or (not obj.solver_opts): + + + + if (obj.input_climate is None) or (obj.delta_t is None) or (obj.S0 is None) or (obj.solver_opts is None): raise ValueError('input_climate, delta_t, S0, and solver_opts attributes must be specified before calling calibrate.') - if not cal_idx: - cal_idx = np.arange(len(Q_obs)) + 1 + # if the list of timesteps to use for calibration is empty, use all timesteps + if cal_idx is None: + cal_idx = np.arange(1,len(Q_obs)+1) if display is None: display = True @@ -541,13 +547,11 @@ def calibrate(obj, Q_obs, cal_idx, optim_fun, par_ini=None, optim_opts=None, of_ if isinstance(cal_idx, list): cal_idx = sorted(cal_idx) - cal_idx_str = ', '.join(map(str, cal_idx)) - if display: print('---') print(f'Starting calibration of model {type(obj)}.') print(f'Simulation will run for timesteps 1-{max(cal_idx)}.') - print(f'Objective function {of_name} will be calculated in time steps {cal_idx_str}.') + print(f'Objective function {of_name} will be calculated in time steps {cal_idx[0]} through {cal_idx[len(cal_idx)-1]}.') print(f'The optimiser {optim_fun} will be used to optimise the objective function.') print('Options passed to the optimiser:') print(optim_opts) @@ -555,23 +559,36 @@ def calibrate(obj, Q_obs, cal_idx, optim_fun, par_ini=None, optim_opts=None, of_ print('---') input_climate_all = obj.input_climate.copy() - obj.input_climate = input_climate_all[:max(cal_idx)] + + # Use the data from the start to the last value of cal_idx to run the model + obj.input_climate.clear() + obj.input_climate['dates'] = input_climate_all['dates'][0:max(cal_idx)] + obj.input_climate['precip'] = input_climate_all['precip'][0:max(cal_idx)] + obj.input_climate['temp'] = input_climate_all['temp'][0:max(cal_idx)] + obj.input_climate['pet'] = input_climate_all['pet'][0:max(cal_idx)] Q_obs = Q_obs[:max(cal_idx)] if par_ini is None: par_ini = np.mean(obj.parRanges, axis=1) + # Helper function to calculate fitness given a set of parameters def fitness_fun(par): - Q_sim = obj.get_streamflow(par=par) - return (-1) ** inverse_flag * of_name(Q_obs, Q_sim, cal_idx, *args) - - [par_opt, of_cal, stopflag, output] = optim_fun(fitness_fun, par_ini, **optim_opts) + Q_sim = obj.get_streamflow(par) + obj_func = getattr(objective_functions, of_name) + result = obj_func(Q_obs, Q_sim, cal_idx, *args) + if inverse_flag: + fitness = (-1) * result[0] + else: + fitness = result[0] + return fitness - of_cal = (-1) ** inverse_flag * of_cal + #Documentation states, to provide sigma values for each parameter, specify sigma = 1, then CMA_stds as multipliers for each parameter + xopt, es = cma.fmin2(fitness_fun, par_ini, 1, {'CMA_stds': optim_opts['insigma'], 'bounds': [optim_opts['LBounds'], optim_opts['UBounds']]}) - obj.input_climate = input_climate_all + obj.input_climate.clear() + obj.input_climate = input_climate_all.copy() - return par_opt, of_cal, stopflag, output + return es.best.x, es.best.f, es.stop(), es.result def default_solver_opts(obj): """ diff --git a/src/pymarrmot/Models/Models/temp.py b/src/pymarrmot/Models/Models/temp.py new file mode 100644 index 0000000..64d5b7a --- /dev/null +++ b/src/pymarrmot/Models/Models/temp.py @@ -0,0 +1,18 @@ +import matplotlib.pyplot as plt +import numpy as np + +x = np.arange(-10, 10, 0.01) +y = x**2 + +#adding text inside the plot +a = float(2.0987543023874523490875423) +a_txt = 'a = %.2f' % round(a, 2) +plt_text = 'line 1\n'+ a_txt +plt.text(-5, 40, plt_text, fontsize = 22) + +plt.plot(x, y, c='g') + +plt.xlabel("X-axis", fontsize = 15) +plt.ylabel("Y-axis",fontsize = 15) + +plt.show(block = True) \ No newline at end of file