A fast (Rcpp) version of the Nash–Sutcliffe model efficiency coefficient

example.r

# Load the libraries
library(Rcpp)
library(hydroGOF)
# Load the source cpp file with a few gof values. 
sourceCpp("gof.cpp")
# Create some data to test and compare the two functions.
obs <- rnorm(10)
sim <- rnorm(10)
# Compare the two functions
NSE_fast(sim,obs)
hydroGOF::NSE(sim,obs)
## And now a more realistic situation with many simulations (50k).
sim_results <- matrix(rnorm(365*50000), 365, 50000)
obs <- rnorm(365)
## Wrap the functions in a system.time function to have a look at the time taken to evaluate each call.
## Using an apply statement to go across the columns of the simulated values.
system.time(nse_fast_values <- apply(sim_results, 2, NSE_fast, obs=obs))
## A comparison with the hydroGOF library
system.time(nse_hydroGOF_values <- apply(sim_results, 2, NSE, obs=obs))

gof.cpp

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
using namespace Rcpp;

/* These functions are Rcpp versions of the goodness of fit functions from
  the hydroGOF package see ?hydroGOF::gof for more info. 
*/

// This is r2 * b (slope)
// [[Rcpp::export]]
double br2_fast(arma::vec sim, arma::vec obs) {
   arma::uvec non_nas = arma::find_finite(sim);
   
   
   arma::colvec coef = arma::solve(obs.elem(non_nas), sim.elem(non_nas));
   double slope = coef(0);
   arma::mat temp = cor(sim.elem(non_nas), obs.elem(non_nas));
   double r2 = pow(temp(0,0),2);
   double res;
   if(std::abs(slope)<=1){
     res = std::abs(slope)*r2;
   } else {
     res = r2/std::abs(slope);
   }
   

   return res;
   
}

// This should be r2
// [[Rcpp::export]]
double r2_fast(arma::vec sim, arma::vec obs) {
  // this locates all real numbers (ie non NAs)
   arma::uvec non_nas = arma::find_finite(sim);
   
   // this finds the coefficients of the linear relationship (not realy needed)
   arma::colvec coef = arma::solve(obs.elem(non_nas), sim.elem(non_nas));
   // this gets the correlation between the two
   arma::mat temp = cor(sim.elem(non_nas), obs.elem(non_nas));
   // then a simply square the correlation to get the r2 value
   double r2 = pow(temp(0,0),2);
   
   return r2;
   
}


//Volumetric efficiency 
// [[Rcpp::export]]
double VE_fast(arma::vec sim, arma::vec obs) {
   arma::uvec non_nas = arma::find_finite(sim);
   
   return 1-accu(abs(obs.elem(non_nas)-sim.elem(non_nas)))/accu(obs.elem(non_nas));
   
}


double minWONA(NumericVector x) {
  int n = x.size();
  int count = 0;
  for(int i = 0; i<n; i++){
    if(NumericVector::is_na(x(i))){
      count++;
    }
  }
  if(count==n){
    return NA_REAL;
  } else {
   return(min(na_omit(x)));
  }
   
}


// [[Rcpp::export]]
double NSE_fast(arma::vec sim, arma::vec obs) {
   arma::uvec non_nas = arma::find_finite(sim);
   return 1 - (accu(pow(obs.elem(non_nas)-sim.elem(non_nas), 2.0)))/ accu(pow(obs.elem(non_nas) - mean(obs.elem(non_nas)),2.0));
   
}


// [[Rcpp::export]]
double tNSE(Rcpp::NumericVector sim, Rcpp::NumericVector obs, int w1= 5.0, int w2 = 5.0, int b=4){
  int n = sim.size();
  Rcpp::NumericMatrix res(n,n);
  Rcpp::NumericMatrix c1(n,n);
  Rcpp::NumericMatrix c2(n,n);
  
  std::fill(res.begin(), res.end(), NA_REAL);
  std::fill(c1.begin(), c1.end(), NA_REAL);
  std::fill(c2.begin(), c2.end(), NA_REAL);
  
  int k = 0;
  for(int j = 0;j<(w1+1); j++){
    c1(j,k) = pow(sim(j)-obs(k),2) + pow(b,2)*pow(j-k,2);
    res(j,k) = c1(j,k);
  }
  double cmin;
  Rcpp::NumericVector tempVect2(2);
  Rcpp::NumericVector tempVect4(4);
  
  for(int k = 1; k<n;k++){
    for(int j = (k-w2); j<= (k+w1); j++){
      
      if(j < 0 || j >= n){
        
      } else {
      
      c2(j,k) = pow(sim(j)-obs(k), 2) + pow(b, 2)* pow(j-k, 2) + c1(j,k-1);
        if(j==0){
          cmin = NA_REAL;
        }
        if(j>0){
          tempVect2(0) = c1(j-1,k-1);
          tempVect2(1) = c2(j-1,k-1);
          cmin = minWONA(tempVect2);
          
        }
        if(j>1){
          tempVect4(0) = c1(j-1,k-1);
          tempVect4(1) = c2(j-1,k-1);
          tempVect4(2) = c1(j-2,k-1);
          tempVect4(3) = c2(j-2,k-1);
          cmin = minWONA(tempVect4);
          
        }
        if(NumericVector::is_na(cmin)){
          c1(j,k) = NA_REAL;
        } else {
          c1(j,k) = pow(sim(j)-obs(k), 2) + pow(b, 2)* pow(j-k, 2) + cmin;
        }
        
      }
    }
  }
  Rcpp::NumericVector easy(2);
  easy(0) = minWONA(c1(_, n-1));
  easy(1) = minWONA(c2(_, n-1));
  double tNSE =  1-(min(easy)/sum(pow(obs-mean(obs),2)));
  
  return(tNSE);
  
}