Basic benchmark for FUSE functions in RHydro and fuse packages

FUSEbenchmark.R

if(!require(RHydro)) install.packages("RHydro",repos="http://R-Forge.R-project.org")
if(!require(fuse)) install.packages("fuse", repos="http://R-Forge.R-project.org")

library(RHydro)
library(fuse)

data(DATA)
data(modlist)

# Set parameter space
DefaultRanges <- data.frame(t(data.frame(c(fusesma.ranges(),fuserouting.ranges()))))
names(DefaultRanges) <- c("Min","Max")

# Sample from the above ranges using the Latin Hypercube Sampling method.
library(tgp)
nRuns <- 100 
parameters <- lhs( nRuns, as.matrix(DefaultRanges) )
parameters <- data.frame(parameters)
names(parameters) <- row.names(DefaultRanges)

counter <- 1
ParameterSet <- parameters[counter,]

# Set up the functions to benchmark
f <- function(DATA,ParameterSets) {
  
  myDELTIM <- 1
  myMID <- 60
  
  for (i in 1:100){
    
    ParameterSet <- ParameterSets[i,]
    
    U <- RHydro::fusesma.sim(DATA,
                             mid=myMID,
                             modlist,
                             deltim=myDELTIM,
                             states=FALSE, fluxes=FALSE, fracstate0=0.25,
                             ParameterSet$rferr_add, 
                             ParameterSet$rferr_mlt, 
                             ParameterSet$frchzne, 
                             ParameterSet$fracten, 
                             ParameterSet$maxwatr_1,
                             ParameterSet$percfrac,
                             ParameterSet$fprimqb, 
                             ParameterSet$qbrate_2a, 
                             ParameterSet$qbrate_2b, 
                             ParameterSet$qb_prms, 
                             ParameterSet$maxwatr_2, 
                             ParameterSet$baserte, 
                             ParameterSet$rtfrac1, 
                             ParameterSet$percrte, 
                             ParameterSet$percexp, 
                             ParameterSet$sacpmlt, 
                             ParameterSet$sacpexp, 
                             ParameterSet$iflwrte,  
                             ParameterSet$axv_bexp, 
                             ParameterSet$sareamax, 
                             ParameterSet$loglamb, 
                             ParameterSet$tishape,
                             ParameterSet$qb_powr)
    
    Q <- RHydro::fuserouting.sim(U, 
                                 mid=myMID, 
                                 modlist=modlist, 
                                 timedelay=ParameterSet$timedelay, 
                                 deltim=myDELTIM)
    
    
  }  
  
}

g <- function(DATA,ParameterSet) {
  
  myDELTIM <- 1
  myMID <- 60
  
  for (i in 1:100){
    
    ParameterSet <- ParameterSets[i,]
    
    U <- fuse::fusesma.sim(DATA,
                           mid=myMID,
                           deltim=myDELTIM,
                           ParameterSet$frchzne, 
                           ParameterSet$fracten, 
                           ParameterSet$maxwatr_1,
                           ParameterSet$percfrac,
                           ParameterSet$fprimqb, 
                           ParameterSet$qbrate_2a, 
                           ParameterSet$qbrate_2b, 
                           ParameterSet$qb_prms, 
                           ParameterSet$maxwatr_2, 
                           ParameterSet$baserte, 
                           ParameterSet$rtfrac1, 
                           ParameterSet$percrte, 
                           ParameterSet$percexp, 
                           ParameterSet$sacpmlt, 
                           ParameterSet$sacpexp, 
                           ParameterSet$iflwrte,  
                           ParameterSet$axv_bexp, 
                           ParameterSet$sareamax, 
                           ParameterSet$loglamb, 
                           ParameterSet$tishape,
                           ParameterSet$qb_powr)
    
    Q <- fuse::fuserouting.sim(U, 
                               mid=myMID, 
                               deltim=myDELTIM,
                               timedelay=ParameterSet$timedelay) 
    
  } 
  
}

# benchmark
library(microbenchmark)
compare <- microbenchmark(f(DATA, parameters), g(DATA, parameters), times = 10)

library(ggplot2)
autoplot(compare)


# Here are my results:
# compare
# Unit: seconds
#                expr        min         lq     median        uq        max neval
# f(DATA, parameters) 423.230827 433.070465 446.845983 451.28512 461.818262    10
# g(DATA, parameters)   2.893856   2.988898   3.076531   3.59736   3.713473    10

# Here is my session info:
# sessionInfo()
# R version 3.1.1 (2014-07-10)
# Platform: x86_64-pc-linux-gnu (64-bit)
#
# locale:
# [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8       
# [4] LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
# [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
#[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
#
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] ggplot2_1.0.0        microbenchmark_1.3-0 tgp_2.4-9            fuse_1.1.0           RHydro_2014-04.1    
# [6] qualV_0.3            KernSmooth_2.23-12   XML_3.98-1.1         deSolve_1.10-9       lhs_0.10            
#[11] sp_1.0-15            xts_0.9-7            zoo_1.7-11          
#
# loaded via a namespace (and not attached):
# [1] colorspace_1.2-4 digest_0.6.4     grid_3.1.1       gtable_0.1.2     lattice_0.20-29  MASS_7.3-33     
# [7] munsell_0.4.2    plyr_1.8.1       proto_0.3-10     Rcpp_0.11.2      reshape2_1.4     scales_0.2.4    
#[13] stringr_0.6.2    tools_3.1.1