From 41f4a94f2c5e011e69849a426ba35c93cf8bcad3 Mon Sep 17 00:00:00 2001 From: ThibHlln Date: Mon, 15 Oct 2018 18:45:40 +0100 Subject: [PATCH] moves obj fn into separate file also adds more citations to relevant work --- hydroeval/__init__.py | 4 +- hydroeval/hydroeval.py | 93 ----------------------- hydroeval/objective_functions.py | 125 +++++++++++++++++++++++++++++++ 3 files changed, 128 insertions(+), 94 deletions(-) create mode 100644 hydroeval/objective_functions.py diff --git a/hydroeval/__init__.py b/hydroeval/__init__.py index a43950c..fd13d2a 100644 --- a/hydroeval/__init__.py +++ b/hydroeval/__init__.py @@ -18,4 +18,6 @@ # You should have received a copy of the GNU General Public License # along with HYDEVA. If not, see . -from .hydroeval import evaluator, nse, nse_c2m, kge, kge_c2m, kgeprime, kgeprime_c2m, rmse, mare, pbias +from .hydroeval import evaluator + +from .objective_functions import nse, nse_c2m, kge, kge_c2m, kgeprime, kgeprime_c2m, rmse, mare, pbias diff --git a/hydroeval/hydroeval.py b/hydroeval/hydroeval.py index 55b94cb..3561037 100644 --- a/hydroeval/hydroeval.py +++ b/hydroeval/hydroeval.py @@ -74,96 +74,3 @@ def evaluator(func, simulation_s, evaluation, axis=1, transform=None, epsilon=No # calculate the requested function return func(my_simu, my_eval) - - -def nse(simulation_s, evaluation): - # calculate Nash-Sutcliffe Efficiency - nse_ = 1 - (np.sum((evaluation - simulation_s) ** 2, axis=1, dtype=np.float64) / - np.sum((evaluation - np.mean(evaluation)) ** 2)) - - return nse_ - - -def nse_c2m(simulation_s, evaluation): - # calculate bounded formulation of NSE following C2M transformation after Mathevet et al. (2006) - nse_ = nse(simulation_s, evaluation) - nse_c2m_ = nse_ / (2 - nse_) - - return nse_c2m_ - - -def kge(simulation_s, evaluation): - # calculate correlation coefficient - sim_mean = np.reshape(np.mean(simulation_s, axis=1), (simulation_s.shape[0], 1)) - obs_mean = np.mean(evaluation) - r = np.sum((simulation_s - sim_mean) * (evaluation - obs_mean), axis=1, dtype=np.float64) / \ - np.sqrt(np.sum((simulation_s - sim_mean) ** 2, axis=1, dtype=np.float64) * - np.sum((evaluation - obs_mean) ** 2, dtype=np.float64)) - # calculate alpha - alpha = np.reshape(np.std(simulation_s, axis=1), (simulation_s.shape[0], 1)) / \ - np.std(evaluation) - # calculate beta - beta = np.reshape(np.sum(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / \ - np.sum(evaluation, dtype=np.float64) - # calculate the Kling-Gupta Efficiency KGE - kge_ = 1 - np.sqrt(np.sum((np.reshape(r, (simulation_s.shape[0], 1)) - 1) ** 2 + - (alpha - 1) ** 2 + (beta - 1) ** 2, axis=1)) - - return np.vstack((kge_, r, alpha[:, 0], beta[:, 0])).T - - -def kge_c2m(simulation_s, evaluation): - # calculate bounded formulation of KGE following C2M transformation after Mathevet et al. (2006) - kge_ = kge(simulation_s, evaluation)[0] - kge_c2m_ = kge_ / (2 - kge_) - - return kge_c2m_ - - -def kgeprime(simulation_s, evaluation): - # calculate correlation coefficient - sim_mean = np.reshape(np.mean(simulation_s, axis=1), (simulation_s.shape[0], 1)) - obs_mean = np.mean(evaluation) - r = np.sum((simulation_s - sim_mean) * (evaluation - obs_mean), axis=1, dtype=np.float64) / \ - np.sqrt(np.sum((simulation_s - sim_mean) ** 2, axis=1, dtype=np.float64) * - np.sum((evaluation - obs_mean) ** 2, dtype=np.float64)) - # calculate gamma - gamma = (np.reshape(np.std(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / sim_mean) / \ - (np.std(evaluation, dtype=np.float64) / obs_mean) - # calculate beta - beta = np.reshape(np.sum(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / \ - np.sum(evaluation, dtype=np.float64) - # calculate the modified Kling-Gupta Efficiency KGE' - kge_ = 1 - np.sqrt(np.sum((np.reshape(r, (simulation_s.shape[0], 1)) - 1) ** 2 + - (gamma - 1) ** 2 + (beta - 1) ** 2, axis=1)) - - return np.vstack((kge_, r, gamma[:, 0], beta[:, 0])).T - - -def kgeprime_c2m(simulation_s, evaluation): - # calculate bounded formulation of KGE' following C2M transformation after Mathevet et al. (2006) - kgeprime_ = kgeprime(simulation_s, evaluation)[0] - kgeprime_c2m_ = kgeprime_ / (2 - kgeprime_) - - return kgeprime_c2m_ - - -def rmse(simulation_s, evaluation): - # calculate root mean square error - rmse_ = np.sqrt(np.mean((evaluation - simulation_s) ** 2, axis=1, dtype=np.float64)) - - return rmse_ - - -def mare(simulation_s, evaluation): - # calculate mean absolute relative error (MARE) - mare_ = np.sum(np.abs(evaluation - simulation_s), axis=1, dtype=np.float64) / np.sum(evaluation) - - return mare_ - - -def pbias(simulation_s, evaluation): - # calculate percent bias - pbias_ = 100 * np.sum(evaluation - simulation_s, axis=1, dtype=np.float64) / np.sum(evaluation) - - return pbias_ diff --git a/hydroeval/objective_functions.py b/hydroeval/objective_functions.py new file mode 100644 index 0000000..37f765f --- /dev/null +++ b/hydroeval/objective_functions.py @@ -0,0 +1,125 @@ +# -*- coding: utf-8 -*- + +# This file is part of HydroEval: An Evaluator for Hydrological Time Series +# Copyright (C) 2018 Thibault Hallouin (1) +# +# (1) Dooge Centre for Water Resources Research, University College Dublin, Ireland +# +# HydroEval is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# HydroEval is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with HydroEval. If not, see . + +import numpy as np + + +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# OBJECTIVE FUNCTIONS +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# Nash-Sutcliffe Efficiency (Nash and Sutcliffe 1970 - https://doi.org/10.1016/0022-1694(70)90255-6) +def nse(simulation_s, evaluation): + nse_ = 1 - (np.sum((evaluation - simulation_s) ** 2, axis=1, dtype=np.float64) / + np.sum((evaluation - np.mean(evaluation)) ** 2)) + + return nse_ + + +# Original Kling-Gupta Efficiency (Gupta et al. 2009 - https://doi.org/10.1016/j.jhydrol.2009.08.003) +def kge(simulation_s, evaluation): + # calculate correlation coefficient + sim_mean = np.reshape(np.mean(simulation_s, axis=1), (simulation_s.shape[0], 1)) + obs_mean = np.mean(evaluation) + r = np.sum((simulation_s - sim_mean) * (evaluation - obs_mean), axis=1, dtype=np.float64) / \ + np.sqrt(np.sum((simulation_s - sim_mean) ** 2, axis=1, dtype=np.float64) * + np.sum((evaluation - obs_mean) ** 2, dtype=np.float64)) + # calculate alpha + alpha = np.reshape(np.std(simulation_s, axis=1), (simulation_s.shape[0], 1)) / \ + np.std(evaluation) + # calculate beta + beta = np.reshape(np.sum(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / \ + np.sum(evaluation, dtype=np.float64) + # calculate the Kling-Gupta Efficiency KGE + kge_ = 1 - np.sqrt(np.sum((np.reshape(r, (simulation_s.shape[0], 1)) - 1) ** 2 + + (alpha - 1) ** 2 + (beta - 1) ** 2, axis=1)) + + return np.vstack((kge_, r, alpha[:, 0], beta[:, 0])).T + + +# Modified kling-Gupta Efficiency (kling et al. 2012 - https://doi.org/10.1016/j.jhydrol.2012.01.011) +def kgeprime(simulation_s, evaluation): + # calculate correlation coefficient + sim_mean = np.reshape(np.mean(simulation_s, axis=1), (simulation_s.shape[0], 1)) + obs_mean = np.mean(evaluation) + r = np.sum((simulation_s - sim_mean) * (evaluation - obs_mean), axis=1, dtype=np.float64) / \ + np.sqrt(np.sum((simulation_s - sim_mean) ** 2, axis=1, dtype=np.float64) * + np.sum((evaluation - obs_mean) ** 2, dtype=np.float64)) + # calculate gamma + gamma = (np.reshape(np.std(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / sim_mean) / \ + (np.std(evaluation, dtype=np.float64) / obs_mean) + # calculate beta + beta = np.reshape(np.sum(simulation_s, axis=1, dtype=np.float64), (simulation_s.shape[0], 1)) / \ + np.sum(evaluation, dtype=np.float64) + # calculate the modified Kling-Gupta Efficiency KGE' + kge_ = 1 - np.sqrt(np.sum((np.reshape(r, (simulation_s.shape[0], 1)) - 1) ** 2 + + (gamma - 1) ** 2 + (beta - 1) ** 2, axis=1)) + + return np.vstack((kge_, r, gamma[:, 0], beta[:, 0])).T + + +# Root Mean Square Error +def rmse(simulation_s, evaluation): + rmse_ = np.sqrt(np.mean((evaluation - simulation_s) ** 2, axis=1, dtype=np.float64)) + + return rmse_ + + +# Mean Absolute Relative Error +def mare(simulation_s, evaluation): + mare_ = np.sum(np.abs(evaluation - simulation_s), axis=1, dtype=np.float64) / np.sum(evaluation) + + return mare_ + + +# Percent Bias +def pbias(simulation_s, evaluation): + pbias_ = 100 * np.sum(evaluation - simulation_s, axis=1, dtype=np.float64) / np.sum(evaluation) + + return pbias_ + +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# BOUNDED VERSIONS OF SOME OBJECTIVE FUNCTIONS +# (After Mathevet et al. 2006 - https://iahs.info/uploads/dms/13614.21--211-219-41-MATHEVET.pdf) +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +# Bounded Version of the Nash-Sutcliffe Efficiency +def nse_c2m(simulation_s, evaluation): + nse_ = nse(simulation_s, evaluation) + nse_c2m_ = nse_ / (2 - nse_) + + return nse_c2m_ + + +# Bounded Version of the Original Kling-Gupta Efficiency +def kge_c2m(simulation_s, evaluation): + kge_ = kge(simulation_s, evaluation)[0] + kge_c2m_ = kge_ / (2 - kge_) + + return kge_c2m_ + + +# Bounded Version of the Modified Kling-Gupta Efficiency +def kgeprime_c2m(simulation_s, evaluation): + kgeprime_ = kgeprime(simulation_s, evaluation)[0] + kgeprime_c2m_ = kgeprime_ / (2 - kgeprime_) + + return kgeprime_c2m_