Calling plot_results does not produce the sim-vs-obs ggof graph when fn!=hydromod because read_results can not obtain the best model output time series.

When fn=hydromod, fn yields a list containing 1) the simulated values and 2) the GoF value. For each particle, the simulated values are assigned and then written into Model_Out.txt. When calling plot_result, this file is then read in by read_out which selects the best model output time series by min/max GoF.

If fn!=hydromod, fn should "return a scalar result" which represents "both model output and its corresponding goodness-of-fit measure". In this case model output is set to GoF, and read_out can not obtain the best model output because it isn't saved for each particle.

Some ideas for allowing plot_results to work with output created by fn!=hydromod:

1. Read best-sim from file in read_out
   Analog to obs, read the file BestModel_out.txt (which is also created if fn=hydromod) in case Model_out.txt contains only GoF. This would also provide an easier/faster alternative to obtaining best-sim from ModelOut.txt.

2. Allow fn to return both sim and gof
   Analog to hydromod, save model output per particles. Would enable read_out to look up the best model output.

3. Add sim (and obs) as optional arguments to plot_results
   Simple and quick fix, as in 774d01d. It would then be up to the user to directly supply sim and obs.

1) and 2) would still require the user to call another custom function which runs the model with the best paramters returned from hydroPSO and writes out the best-sim and obs text files before being able to call plot_results. Doing this from hydroPSO directly is currently not possible because unlike with fn=hydromod, the obs parameter name is unknown (and fn currently does not return sim). I suppose it would be possible to pull obs from the ellipsis if it exists or create a stripped-down version of hydromod for running R code instead of external commands.

In general, I think it would be great to make saving model output per particle optional and implement 1) as an alternative way to supply the best model out. One would lose the ecdf plots, but improve performance by not writing the full sim to disk each model run and having to process the resulting ModelOut.txt which can easily grow to several gigabytes. This would possibly involve setting model.out to NA or similar instead of copying GoF, because at the moment plot_results exits with error once it hits the ecdf plots if model.values contains GoF instead of the full sim time series.

Hello,

In the lhoat function, the REPORT variable as set in the control list, seems to not change any of the output printed to the console:

library(hydroPSO)
#> (C) 2011-2020 M. Zambrano-Bigiarini and R. Rojas (GPL >=2 license)
#> Type 'citation('hydroPSO')' to see how to cite this package

# using example 1 in help(lhoat)
nparam <- 3
linear <- function(x) x[1] + x[2] + x[3]
lower <- c(0.5, 1.5, 4.5)
upper <- c(1.5, 4.5, 13.5)
names(lower) <- c("a","b","c")

# 200 parameter sets, print every 100 (default) but does not print any to console
set.seed(123) 
tmp <- lhoat(
  fn=linear, 
  lower=lower, 
  upper=upper,
  control=list(N=50, f=0.015, write2disk=FALSE, verbose=TRUE)
)
#> 
#> [ Number of strata for each parameter (N) : 50 ]
#> 
#> [ Number of Parameter Sets to be run      : 200 ]
#> 
#> ==============================================================
#> [  2)   Initial  LHS ...                                     ]
#> ==============================================================
#> 
#> ==============================================================
#> [  3)   Running LH-OAT ...                                   ]
#> ==============================================================
#>                              |                                
#>                              |
#> ==============================================================
#> [==================    LH-OAT finished !    =================]
#> ==============================================================

# 200 parameter sets, print every 10 but does not print any to console
set.seed(123) 
tmp <- lhoat(
  fn=linear, 
  lower=lower, 
  upper=upper,
  control=list(N=50, f=0.015, write2disk=FALSE, verbose=TRUE, REPORT = 10)
)
#> 
#> [ Number of strata for each parameter (N) : 50 ]
#> 
#> [ Number of Parameter Sets to be run      : 200 ]
#> 
#> ==============================================================
#> [  2)   Initial  LHS ...                                     ]
#> ==============================================================
#> 
#> ==============================================================
#> [  3)   Running LH-OAT ...                                   ]
#> ==============================================================
#>                              |                                
#>                              |
#> ==============================================================
#> [==================    LH-OAT finished !    =================]
#> ==============================================================

^{Created on 2021-01-13 by the reprex package (v0.3.0)}

devtools::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value                       
#>  version  R version 4.0.3 (2020-10-10)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  Dutch_Belgium.1252          
#>  ctype    Dutch_Belgium.1252          
#>  tz       Europe/Paris                
#>  date     2021-01-13                  
#> 
#> - Packages -------------------------------------------------------------------
#>  package       * version date       lib source        
#>  assertthat      0.2.1   2019-03-21 [1] CRAN (R 4.0.2)
#>  automap         1.0-14  2013-08-29 [1] CRAN (R 4.0.2)
#>  backports       1.2.1   2020-12-09 [1] CRAN (R 4.0.3)
#>  base64enc       0.1-3   2015-07-28 [1] CRAN (R 4.0.0)
#>  callr           3.5.1   2020-10-13 [1] CRAN (R 4.0.3)
#>  checkmate       2.0.0   2020-02-06 [1] CRAN (R 4.0.2)
#>  class           7.3-17  2020-04-26 [2] CRAN (R 4.0.3)
#>  cli             2.2.0   2020-11-20 [1] CRAN (R 4.0.3)
#>  cluster         2.1.0   2019-06-19 [2] CRAN (R 4.0.3)
#>  colorspace      2.0-0   2020-11-11 [1] CRAN (R 4.0.3)
#>  crayon          1.3.4   2017-09-16 [1] CRAN (R 4.0.2)
#>  data.table      1.13.6  2020-12-30 [1] CRAN (R 4.0.3)
#>  desc            1.2.0   2018-05-01 [1] CRAN (R 4.0.2)
#>  devtools        2.3.2   2020-09-18 [1] CRAN (R 4.0.3)
#>  digest          0.6.27  2020-10-24 [1] CRAN (R 4.0.3)
#>  dplyr           1.0.2   2020-08-18 [1] CRAN (R 4.0.3)
#>  e1071           1.7-4   2020-10-14 [1] CRAN (R 4.0.3)
#>  ellipsis        0.3.1   2020-05-15 [1] CRAN (R 4.0.2)
#>  evaluate        0.14    2019-05-28 [1] CRAN (R 4.0.2)
#>  fansi           0.4.1   2020-01-08 [1] CRAN (R 4.0.2)
#>  FNN             1.1.3   2019-02-15 [1] CRAN (R 4.0.2)
#>  foreign         0.8-80  2020-05-24 [2] CRAN (R 4.0.3)
#>  Formula         1.2-4   2020-10-16 [1] CRAN (R 4.0.3)
#>  fs              1.5.0   2020-07-31 [1] CRAN (R 4.0.3)
#>  generics        0.1.0   2020-10-31 [1] CRAN (R 4.0.3)
#>  ggplot2         3.3.3   2020-12-30 [1] CRAN (R 4.0.3)
#>  glue            1.4.2   2020-08-27 [1] CRAN (R 4.0.3)
#>  gridExtra       2.3     2017-09-09 [1] CRAN (R 4.0.2)
#>  gstat           2.0-6   2020-05-18 [1] CRAN (R 4.0.2)
#>  gtable          0.3.0   2019-03-25 [1] CRAN (R 4.0.2)
#>  highr           0.8     2019-03-20 [1] CRAN (R 4.0.2)
#>  Hmisc           4.4-2   2020-11-29 [1] CRAN (R 4.0.3)
#>  htmlTable       2.1.0   2020-09-16 [1] CRAN (R 4.0.3)
#>  htmltools       0.5.0   2020-06-16 [1] CRAN (R 4.0.2)
#>  htmlwidgets     1.5.3   2020-12-10 [1] CRAN (R 4.0.3)
#>  hydroGOF        0.4-0   2020-03-12 [1] CRAN (R 4.0.2)
#>  hydroPSO      * 0.5-1   2020-04-29 [1] CRAN (R 4.0.2)
#>  hydroTSM        0.6-0   2020-03-11 [1] CRAN (R 4.0.2)
#>  intervals       0.15.2  2020-04-04 [1] CRAN (R 4.0.0)
#>  jpeg            0.1-8.1 2019-10-24 [1] CRAN (R 4.0.0)
#>  knitr           1.30    2020-09-22 [1] CRAN (R 4.0.3)
#>  lattice         0.20-41 2020-04-02 [2] CRAN (R 4.0.3)
#>  latticeExtra    0.6-29  2019-12-19 [1] CRAN (R 4.0.2)
#>  lhs             1.1.1   2020-10-05 [1] CRAN (R 4.0.3)
#>  lifecycle       0.2.0   2020-03-06 [1] CRAN (R 4.0.2)
#>  magrittr        2.0.1   2020-11-17 [1] CRAN (R 4.0.3)
#>  maptools        1.0-2   2020-08-24 [1] CRAN (R 4.0.3)
#>  Matrix          1.2-18  2019-11-27 [2] CRAN (R 4.0.3)
#>  memoise         1.1.0   2017-04-21 [1] CRAN (R 4.0.2)
#>  munsell         0.5.0   2018-06-12 [1] CRAN (R 4.0.2)
#>  nnet            7.3-14  2020-04-26 [2] CRAN (R 4.0.3)
#>  pillar          1.4.7   2020-11-20 [1] CRAN (R 4.0.3)
#>  pkgbuild        1.2.0   2020-12-15 [1] CRAN (R 4.0.3)
#>  pkgconfig       2.0.3   2019-09-22 [1] CRAN (R 4.0.2)
#>  pkgload         1.1.0   2020-05-29 [1] CRAN (R 4.0.2)
#>  plyr            1.8.6   2020-03-03 [1] CRAN (R 4.0.2)
#>  png             0.1-7   2013-12-03 [1] CRAN (R 4.0.0)
#>  prettyunits     1.1.1   2020-01-24 [1] CRAN (R 4.0.2)
#>  processx        3.4.5   2020-11-30 [1] CRAN (R 4.0.3)
#>  ps              1.5.0   2020-12-05 [1] CRAN (R 4.0.3)
#>  purrr           0.3.4   2020-04-17 [1] CRAN (R 4.0.2)
#>  R6              2.5.0   2020-10-28 [1] CRAN (R 4.0.3)
#>  RColorBrewer    1.1-2   2014-12-07 [1] CRAN (R 4.0.0)
#>  Rcpp            1.0.5   2020-07-06 [1] CRAN (R 4.0.2)
#>  remotes         2.2.0   2020-07-21 [1] CRAN (R 4.0.3)
#>  reshape         0.8.8   2018-10-23 [1] CRAN (R 4.0.2)
#>  rlang           0.4.10  2020-12-30 [1] CRAN (R 4.0.3)
#>  rmarkdown       2.6     2020-12-14 [1] CRAN (R 4.0.3)
#>  rpart           4.1-15  2019-04-12 [2] CRAN (R 4.0.3)
#>  rprojroot       2.0.2   2020-11-15 [1] CRAN (R 4.0.3)
#>  rstudioapi      0.13    2020-11-12 [1] CRAN (R 4.0.3)
#>  scales          1.1.1   2020-05-11 [1] CRAN (R 4.0.2)
#>  scatterplot3d   0.3-41  2018-03-14 [1] CRAN (R 4.0.0)
#>  sessioninfo     1.1.1   2018-11-05 [1] CRAN (R 4.0.2)
#>  sm              2.2-5.6 2018-09-27 [1] CRAN (R 4.0.2)
#>  sp              1.4-4   2020-10-07 [1] CRAN (R 4.0.3)
#>  spacetime       1.2-3   2020-01-21 [1] CRAN (R 4.0.2)
#>  stringi         1.5.3   2020-09-09 [1] CRAN (R 4.0.3)
#>  stringr         1.4.0   2019-02-10 [1] CRAN (R 4.0.2)
#>  survival        3.2-7   2020-09-28 [2] CRAN (R 4.0.3)
#>  testthat        3.0.1   2020-12-17 [1] CRAN (R 4.0.3)
#>  tibble          3.0.4   2020-10-12 [1] CRAN (R 4.0.3)
#>  tidyselect      1.1.0   2020-05-11 [1] CRAN (R 4.0.2)
#>  usethis         2.0.0   2020-12-10 [1] CRAN (R 4.0.3)
#>  vctrs           0.3.6   2020-12-17 [1] CRAN (R 4.0.3)
#>  vioplot         0.3.5   2020-06-15 [1] CRAN (R 4.0.2)
#>  withr           2.3.0   2020-09-22 [1] CRAN (R 4.0.3)
#>  xfun            0.20    2021-01-06 [1] CRAN (R 4.0.3)
#>  xts             0.12.1  2020-09-09 [1] CRAN (R 4.0.3)
#>  yaml            2.2.1   2020-02-01 [1] CRAN (R 4.0.0)
#>  zoo             1.8-8   2020-05-02 [1] CRAN (R 4.0.2)
#> 
#> [1] C:/Users/casne/OneDrive/Documenten/R/library_cas
#> [2] C:/Users/casne/OneDrive/Documenten/R/R-4.0.3/library

I have following function to maximize with 18<=x1<=20; 4000<=x2<=4800; 5<=x3<=13
#######################################################
f<-function(x1,x2,x3){
-53.06453125+5.146767857x1+0.001306786x2+0.25028125x3-0.00019125x1x2-
0.00005594x2x3-0.111714286x1^2+0.00000032*x2^2
}
################################################################
I am running following code
lower<-c(18,4000,5)
upper<-c(20,4800,13)
out<-hydroPSO(fn=f, lower=lower, upper=upper, control=list(write2disk=FALSE) )

But i am getting the following Error
Error in FUN(newX[, i], ...) : argument "x2" is missing, with no default
Any help would be Appriciated

Raven is a good hydrological modelling platform (http://raven.uwaterloo.ca/) and is still in fast development. I would suggest you to provide some examples about using hydroPSO with Raven, so that people will have more alternatives for auto calibration.

Kangsheng

I have tried the example of hydroPSO-MF2005 but I always got an error message as below. The version of R I used is 3.4.4. I have tried on both Win 7 and Win 10 but same error remains. Can anyone suggest what I should do to find out the issue? Thanks in advance for help.

Error in file(con, "r") : cannot open the connection
In addition: There were 28 warnings (use warnings() to see them)

There are additional warning messages as below.

1: running command '"run_me.bat"' had status 1
2: In readLines(fname) : line 22 appears to contain an embedded nul
3: In readLines(fname) : line 58 appears to contain an embedded nul
4: In readLines(fname) : line 90 appears to contain an embedded nul
5: In readLines(fname) : line 124 appears to contain an embedded nul
6: In readLines(fname) : line 22 appears to contain an embedded nul
7: In readLines(fname) : line 58 appears to contain an embedded nul
8: In readLines(fname) : line 90 appears to contain an embedded nul
9: In readLines(fname) : line 124 appears to contain an embedded nul
10: running command '"run_me.bat"' had status 1
11: In readLines(fname) : line 22 appears to contain an embedded nul
12: In readLines(fname) : line 58 appears to contain an embedded nul
13: In readLines(fname) : line 90 appears to contain an embedded nul
14: In readLines(fname) : line 124 appears to contain an embedded nul
15: In readLines(fname) : line 22 appears to contain an embedded nul
16: In readLines(fname) : line 58 appears to contain an embedded nul
17: In readLines(fname) : line 90 appears to contain an embedded nul
18: In readLines(fname) : line 124 appears to contain an embedded nul
19: running command '"run_me.bat"' had status 1
20: In readLines(fname) : line 22 appears to contain an embedded nul
21: In readLines(fname) : line 58 appears to contain an embedded nul
22: In readLines(fname) : line 90 appears to contain an embedded nul
23: In readLines(fname) : line 124 appears to contain an embedded nul
24: In readLines(fname) : line 22 appears to contain an embedded nul
25: In readLines(fname) : line 58 appears to contain an embedded nul
26: In readLines(fname) : line 90 appears to contain an embedded nul
27: In readLines(fname) : line 124 appears to contain an embedded nul
28: In file(con, "r") : cannot open file 'pre_param.txt': Permission denied
There were 14 warnings (use warnings() to see them)

warnings()
警告信息：
1: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 23 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/LocalBestPerIter.txt)
2: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 22 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/PbestPerIter.txt)
3: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 21 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/BestParamPerIter.txt)
4: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 20 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/ConvergenceMeasures.txt)
5: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 19 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Velocities.txt)
6: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 18 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Particles.txt)
7: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 17 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Model_out.txt)
8: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 9 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/LocalBestPerIter.txt)
9: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 8 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/PbestPerIter.txt)
10: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 7 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/BestParamPerIter.txt)
11: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 6 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/ConvergenceMeasures.txt)
12: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 5 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Velocities.txt)
13: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 4 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Particles.txt)
14: In tryCatchOne(expr, names, parentenv, handlers[[1L]]) :
closing unused connection 3 (D:/mynutstore/Reading/hydropso/MF2005/PSO.70p3000i.rand_TVc1_TVlambda.out/Model_out.txt)

The hydroPSO function seems to be broken for hydroPSO 0.5-0 as installed from CRAN.

See example from hydroPSO docs:

library(hydroPSO)
#> (C) 2011-2020 M. Zambrano-Bigiarini and R. Rojas (GPL >=2 license)
#> Type 'citation('hydroPSO')' to see how to cite this package

# Number of dimensions of the optimisation problem (for all the examples)
D <- 5

# Boundaries of the search space (Rastrigin function)
lower <- rep(-5.12, D)
upper <- rep(5.12, D)

## Not run: 
# Setting the home directory of the user as working directory 
setwd("~")

################################ 
# Example 1. Basic use         #
################################ 

# Setting the seed (for reproducible results)         
set.seed(100)

# Basic use 1. Rastrigin function (non-linear and multi-modal with many local minima)
# Results are not saved to the hard disk, for faster execution ('write2disk=FALSE')
hydroPSO(fn=rastrigin, lower=lower, upper=upper, control=list(write2disk=FALSE) )
#> 
#> ================================================================================
#> [                                Initialising  ...                             ]
#> ================================================================================
#> 
#> [npart=40 ; maxit=1000 ; method=spso2011 ; topology=random ; boundary.wall=absorbing2011]
#> 
#> [ user-definitions in control: write2disk=FALSE ]
#> 
#> 
#> ================================================================================
#> [                                 Running  PSO ...                             ]
#> ================================================================================
#> 
#> Error in hydroPSO(fn = rastrigin, lower = lower, upper = upper, control = list(write2disk = FALSE)): object 'GoF' not found

^{Created on 2020-04-23 by the reprex package (v0.3.0)}

devtools::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.6.3 (2020-02-29)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  Dutch_Belgium.1252          
#>  ctype    Dutch_Belgium.1252          
#>  tz       Europe/Paris                
#>  date     2020-04-23                  
#> 
#> - Packages -------------------------------------------------------------------
#>  package       * version date       lib source        
#>  acepack         1.4.1   2016-10-29 [1] CRAN (R 3.6.1)
#>  assertthat      0.2.1   2019-03-21 [1] CRAN (R 3.6.1)
#>  automap         1.0-14  2013-08-29 [1] CRAN (R 3.6.1)
#>  backports       1.1.5   2019-10-02 [1] CRAN (R 3.6.1)
#>  base64enc       0.1-3   2015-07-28 [1] CRAN (R 3.6.0)
#>  callr           3.4.2   2020-02-12 [1] CRAN (R 3.6.3)
#>  checkmate       2.0.0   2020-02-06 [1] CRAN (R 3.6.3)
#>  class           7.3-15  2019-01-01 [2] CRAN (R 3.6.3)
#>  cli             2.0.2   2020-02-28 [1] CRAN (R 3.6.3)
#>  cluster         2.1.0   2019-06-19 [2] CRAN (R 3.6.3)
#>  colorspace      1.4-1   2019-03-18 [1] CRAN (R 3.6.1)
#>  crayon          1.3.4   2017-09-16 [1] CRAN (R 3.6.1)
#>  data.table      1.12.8  2019-12-09 [1] CRAN (R 3.6.3)
#>  desc            1.2.0   2018-05-01 [1] CRAN (R 3.6.1)
#>  devtools        2.2.2   2020-02-17 [1] CRAN (R 3.6.3)
#>  digest          0.6.25  2020-02-23 [1] CRAN (R 3.6.3)
#>  dplyr           0.8.5   2020-03-07 [1] CRAN (R 3.6.3)
#>  e1071           1.7-3   2019-11-26 [1] CRAN (R 3.6.3)
#>  ellipsis        0.3.0   2019-09-20 [1] CRAN (R 3.6.1)
#>  evaluate        0.14    2019-05-28 [1] CRAN (R 3.6.1)
#>  fansi           0.4.1   2020-01-08 [1] CRAN (R 3.6.3)
#>  FNN             1.1.3   2019-02-15 [1] CRAN (R 3.6.1)
#>  foreign         0.8-75  2020-01-20 [2] CRAN (R 3.6.3)
#>  Formula         1.2-3   2018-05-03 [1] CRAN (R 3.6.0)
#>  fs              1.3.2   2020-03-05 [1] CRAN (R 3.6.3)
#>  ggplot2         3.3.0   2020-03-05 [1] CRAN (R 3.6.3)
#>  glue            1.3.2   2020-03-12 [1] CRAN (R 3.6.3)
#>  gridExtra       2.3     2017-09-09 [1] CRAN (R 3.6.1)
#>  gstat           2.0-4   2020-01-21 [1] CRAN (R 3.6.3)
#>  gtable          0.3.0   2019-03-25 [1] CRAN (R 3.6.1)
#>  highr           0.8     2019-03-20 [1] CRAN (R 3.6.1)
#>  Hmisc           4.3-1   2020-02-07 [1] CRAN (R 3.6.3)
#>  htmlTable       1.13.3  2019-12-04 [1] CRAN (R 3.6.3)
#>  htmltools       0.4.0   2019-10-04 [1] CRAN (R 3.6.1)
#>  htmlwidgets     1.5.1   2019-10-08 [1] CRAN (R 3.6.1)
#>  hydroGOF        0.4-0   2020-03-12 [1] CRAN (R 3.6.3)
#>  hydroPSO      * 0.5-0   2020-03-18 [1] CRAN (R 3.6.3)
#>  hydroTSM        0.6-0   2020-03-11 [1] CRAN (R 3.6.3)
#>  intervals       0.15.1  2015-08-27 [1] CRAN (R 3.6.0)
#>  jpeg            0.1-8.1 2019-10-24 [1] CRAN (R 3.6.1)
#>  knitr           1.28    2020-02-06 [1] CRAN (R 3.6.3)
#>  lattice         0.20-38 2018-11-04 [2] CRAN (R 3.6.3)
#>  latticeExtra    0.6-29  2019-12-19 [1] CRAN (R 3.6.3)
#>  lhs             1.0.1   2019-02-03 [1] CRAN (R 3.6.1)
#>  lifecycle       0.2.0   2020-03-06 [1] CRAN (R 3.6.3)
#>  magrittr        1.5     2014-11-22 [1] CRAN (R 3.6.1)
#>  maptools        0.9-9   2019-12-01 [1] CRAN (R 3.6.3)
#>  Matrix          1.2-18  2019-11-27 [2] CRAN (R 3.6.3)
#>  memoise         1.1.0   2017-04-21 [1] CRAN (R 3.6.1)
#>  munsell         0.5.0   2018-06-12 [1] CRAN (R 3.6.1)
#>  nnet            7.3-12  2016-02-02 [2] CRAN (R 3.6.3)
#>  pillar          1.4.3   2019-12-20 [1] CRAN (R 3.6.3)
#>  pkgbuild        1.0.6   2019-10-09 [1] CRAN (R 3.6.1)
#>  pkgconfig       2.0.3   2019-09-22 [1] CRAN (R 3.6.1)
#>  pkgload         1.0.2   2018-10-29 [1] CRAN (R 3.6.1)
#>  plyr            1.8.6   2020-03-03 [1] CRAN (R 3.6.3)
#>  png             0.1-7   2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits     1.1.1   2020-01-24 [1] CRAN (R 3.6.3)
#>  processx        3.4.2   2020-02-09 [1] CRAN (R 3.6.3)
#>  ps              1.3.2   2020-02-13 [1] CRAN (R 3.6.3)
#>  purrr           0.3.3   2019-10-18 [1] CRAN (R 3.6.1)
#>  R6              2.4.1   2019-11-12 [1] CRAN (R 3.6.3)
#>  RColorBrewer    1.1-2   2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp            1.0.4   2020-03-17 [1] CRAN (R 3.6.3)
#>  remotes         2.1.1   2020-02-15 [1] CRAN (R 3.6.3)
#>  reshape         0.8.8   2018-10-23 [1] CRAN (R 3.6.1)
#>  rlang           0.4.5   2020-03-01 [1] CRAN (R 3.6.3)
#>  rmarkdown       1.16    2019-10-01 [1] CRAN (R 3.6.1)
#>  rpart           4.1-15  2019-04-12 [2] CRAN (R 3.6.3)
#>  rprojroot       1.3-2   2018-01-03 [1] CRAN (R 3.6.1)
#>  rstudioapi      0.11    2020-02-07 [1] CRAN (R 3.6.3)
#>  scales          1.1.0   2019-11-18 [1] CRAN (R 3.6.3)
#>  scatterplot3d   0.3-41  2018-03-14 [1] CRAN (R 3.6.0)
#>  sessioninfo     1.1.1   2018-11-05 [1] CRAN (R 3.6.1)
#>  sm              2.2-5.6 2018-09-27 [1] CRAN (R 3.6.1)
#>  sp              1.4-1   2020-02-28 [1] CRAN (R 3.6.3)
#>  spacetime       1.2-3   2020-01-21 [1] CRAN (R 3.6.3)
#>  stringi         1.4.6   2020-02-17 [1] CRAN (R 3.6.2)
#>  stringr         1.4.0   2019-02-10 [1] CRAN (R 3.6.1)
#>  survival        3.1-8   2019-12-03 [2] CRAN (R 3.6.3)
#>  testthat        2.3.2   2020-03-02 [1] CRAN (R 3.6.3)
#>  tibble          2.1.3   2019-06-06 [1] CRAN (R 3.6.1)
#>  tidyselect      1.0.0   2020-01-27 [1] CRAN (R 3.6.3)
#>  usethis         1.5.1   2019-07-04 [1] CRAN (R 3.6.1)
#>  vioplot         0.3.4   2019-11-29 [1] CRAN (R 3.6.3)
#>  withr           2.1.2   2018-03-15 [1] CRAN (R 3.6.1)
#>  xfun            0.12    2020-01-13 [1] CRAN (R 3.6.3)
#>  xts             0.12-0  2020-01-19 [1] CRAN (R 3.6.3)
#>  yaml            2.2.1   2020-02-01 [1] CRAN (R 3.6.2)
#>  zoo             1.8-7   2020-01-10 [1] CRAN (R 3.6.3)
#> 
#> [1] C:/Users/casne/Documents/R/library_cas
#> [2] C:/Users/casne/Documents/R/R-3.6.3/library

When reinstalling 0.4-1 this example does work:

library(hydroPSO)
#> (C) 2011-2018 M. Zambrano-Bigiarini and R. Rojas (GPL >=2 license)
#> Type 'citation('hydroPSO')' to see how to cite this package

# Number of dimensions of the optimisation problem (for all the examples)
D <- 5

# Boundaries of the search space (Rastrigin function)
lower <- rep(-5.12, D)
upper <- rep(5.12, D)

## Not run: 
# Setting the home directory of the user as working directory 
setwd("~")

################################ 
# Example 1. Basic use         #
################################ 

# Setting the seed (for reproducible results)         
set.seed(100)

# Basic use 1. Rastrigin function (non-linear and multi-modal with many local minima)
# Results are not saved to the hard disk, for faster execution ('write2disk=FALSE')
hydroPSO(fn=rastrigin, lower=lower, upper=upper, control=list(write2disk=FALSE) )
#> 
#> ================================================================================
#> [                                Initialising  ...                             ]
#> ================================================================================
#> 
#> [npart=40 ; maxit=1000 ; method=spso2011 ; topology=random ; boundary.wall=absorbing2011]
#> 
#> [ user-definitions in control: write2disk=FALSE ]
#> 
#> 
#> ================================================================================
#> [                                 Running  PSO ...                             ]
#> ================================================================================
#> 
#> iter: 100  Gbest: 2.992E+00  Gbest_rate:  0.08%  Iter_best_fit: 2.992E+00  nSwarm_Radius: 2.96E-02  |g-mean(p)|/mean(p): 69.36%
#>                                     |
#> ================================================================================
#> [                          Creating the R output ...                           ]
#> ================================================================================
#> $par
#>       Param1       Param2       Param3       Param4       Param5 
#> 4.382399e-05 9.947626e-01 9.950853e-01 9.948407e-01 5.619973e-05 
#> 
#> $value
#> [1] 2.984892
#> 
#> $best.particle
#> [1] 5
#> 
#> $counts
#> function.calls     iterations    regroupings 
#>           5920            148              0 
#> 
#> $convergence
#> [1] 1
#> 
#> $message
#> [1] "Converged ('reltol' criterion)"

^{Created on 2020-04-23 by the reprex package (v0.3.0)}

Even with control=list(verbose=FALSE), hydroPSO produces several lines of output. A (nonsensical example to replicate the issue is)

hydroPSO(fn=function(x){x^4},lower=1,upper=1,control=list(verbose=FALSE))
>[npart=40 ; maxit=1000 ; method=spso2011 ; topology=random ; boundary.wall=absorbing2011]
>       
> [ user-definitions in control: verbose=FALSE ]

On top of this, the whole text appears in red, which makes it harder to spot error messages from other parts of my code. Is there a way this can be fully eliminated?

I have edited the ModifyInputFile to allow for different options for parameter changes. However the parameter changes are not affecting the model output.
ModifyInputFile.txt

The argument method="wfips" appears to not work or the documentation / error messages regarding missing parameters/misuse of the method are unclear. Consider:

temp <- function(x) {x[1]+x[2]}
hydroPSO(fn=temp, lower=c(1000,1000), upper=c(20000,20000), method="wfips") 

Which results in an error message that varies from run to run but fits the general pattern:

Error in is.data.frame(x) : 
  dims [product 2] do not match the length of object [40]

... where the number following product varies.

�The code that triggers the error appears to be:

V[j, ] <- compute.veloc(x = X[j, ], v = V[j, ], w = w, c1 = c1, 
    c2 = c2, CF = CF, Pbest = X.best.part, part.index = j, gbest = X.best.part[gbest.pos, 
        ], topology = ltopology, method = method, MinMax = MinMax, 
    neighs.index = X.neighbours[j, ], localBest = X.best.part[LocalBest.pos[j], 
        ], localBest.pos = LocalBest.pos[j], ngbest.fit = ngbest.fit, 
    ngbest = X.best.part[ngbest.pos, ], lpbest.fit = pbest.fit[X.neighbours[j, 
        ]])

I am working on inverting SWAP model (http://www.swap.alterra.nl) using satellite soil moisture to upscale soil hydraulic parameters and want to use hydroPSO for model inversion.

I have developed a wrapper function in R for updating and running the SWAP parameters and use the "hydromodInR" option to link with hydroPSO. However, when I use the option "write2disk=TRUE" the optimization stops at the next iteration and I receive the following error (no parameters are saved on disk).

Error in writeLines(as.character(c(iter, j, formatC(GoF, format = "E", : invalid connection
The code otherwise works fine if I use "write2disk=FALSE".
How can I save the optimized parameters for each iteration to extract 95% CI of the optimized params?

Hi,
In the paper A model-independent Particle Swarm Optimisation software for model calibration it is mentioned that multi-objective optimization is under development (year 2013).

Has this functionality been implemented? Is there a work around?

Cheers

How to determine npart--number of particles in the swarm? it seems that npart =40 is used as default.

Does npart depend on number of parameters used in calibration or other factors? For example, I am calibrating VIC model, which is a grid based model. I have 6 parameters to be calibrated, but I need to repeat 100 times for each parameter because I have 100 grids. Each grid is independent and is calibrated independently, although I use the same value for each grid at each calibration run.

So what npart should I use? 5, 40, or 100 in such case?

So far, I have not found a way to suggest the positions of particles.
Would this make sense anyway?
I intended to start particles at points, where I assume a good fit close by and by that shorten the calculation time.
Would it then be possible, to suggest some positions and generate others similar to a default run?

the example code source: https://zenodo.org/record/1044951
program version: latest
error info:
Error in checkForRemoteErrors(val) : 10 nodes produced errors; first error: line 8 did not have 7 elements Calls: hydroPSO ... <Anonymous> -> staticClusterApply -> checkForRemoteErrors

could someone help about this problem? I have got into this problem for days.

According to to the documentation, there is no way to pass arguments to model.FUN when fn is set to "hydromodInR":

From how I understand the intention behind "hydromodInR" (arg checking etc) I believe model.FUN.args should be used here instead.

If you like I can make a PR if you let me know how you want this resolved.

BasinObs<-read.csv("Monthly Katar.csv",header=TRUE)
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR2M, DatesR = as.POSIXlt(BasinObs$DatesR),
Precip = BasinObs$P, PotEvap = BasinObs$E)
Ind_Run <- seq(which(format(as.POSIXlt(BasinObs$DatesR), format = "%m/%Y")=="01/1986"),
which(format(as.POSIXlt(BasinObs$DatesR), format = "%m/%Y")=="11/1989"))

RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)

InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])

CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR2M, FUN_CALIB = Calibration_Michel)

OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR2M, FUN_CRIT = ErrorCrit_NSE)
OptimGR2M <- function(ParamOptim) {

Transformation of the parameter set to real space

RawParamOptim <- airGR::TransfoParam_GR2M(ParamIn = ParamOptim,
Direction = "TR")

Simulation given a parameter set

OutputsModel <- airGR::RunModel_GR2M(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = RawParamOptim)

Computation of the value of the performance criteria

OutputsCrit <- airGR::ErrorCrit_RMSE(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
verbose = FALSE)
return(OutputsCrit$CritValue)
}
lowerGR2M <- rep(0, times = 2)
upperGR2M <- rep(+999.99, times = 2)
startGR2M <- c(500, 2)
optPORT <- stats::nlminb(start = startGR2M,
objective = OptimGR2M,
lower = lowerGR2M, upper = upperGR2M,
control = list(trace = 1))
startGR2M <- expand.grid(data.frame(CalibOptions$StartParamDistrib))
optPORT_ <- function(x) {
opt <- stats::nlminb(start = x,
objective = OptimGR2M,
lower = lowerGR2M, upper = upperGR2M,
control = list(trace = 1))
}
listOptPORT <- apply(startGR2M, MARGIN = 1, FUN = optPORT_)
parPORT <- t(sapply(listOptPORT, function(x) x$par))
objPORT <- sapply(listOptPORT, function(x) x$objective)
resPORT <- data.frame(parPORT, RMSE = objPORT)
summary(resPORT)
library(airGR)
library(DEoptim)
library(hydroPSO)
library(Rmalschains)
optDE <- DEoptim::DEoptim(fn = OptimGR2M,
lower = lowerGR2M, upper = upperGR2M,
control = DEoptim::DEoptim.control(NP = 40, trace = 10))
optPSO <- hydroPSO::hydroPSO(fn = OptimGR2M,
lower = lowerGR2M, upper = upperGR2M,
control = list(write2disk = FALSE, verbose = FALSE))
optMALS <- Rmalschains::malschains(fn = OptimGR2M,
lower = lowerGR2M, upper = upperGR2M,
maxEvals = 2000)
resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
round(rbind(
OutputsCalib$ParamFinalR
airGR::TransfoParam_GR2M(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
airGR::TransfoParam_GR2M(ParamIn = as.numeric(optPSO$par) , Direction = "TR"),
airGR::TransfoParam_GR2M(ParamIn = optMALS$sol , Direction = "TR")),
digits = 3))
resGLOB

But I keep getting this error

Error in if (any(ProdStore < 0) | any(RoutStore < 0) | any(UH1 < 0) | :
missing value where TRUE/FALSE needed
Called from: CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L],
ExpStore = NULL, UH1 = NULL, UH2 = NULL, GCemaNeigeLayers = NULL,
eTGCemaNeigeLayers = NULL, verbose = FALSE)
Browse[1]> Q

Can you please help?

out <- hydroPSO(fn="BGCM.basic",
lower=lower,
upper=upper,
method="spso2011",
control=list(write2disk=TRUE, MinMax="max", npart=10,
maxit=50, normalise=TRUE, REPORT=1, #parallel="parallelWin",
reltol=1E-10)
)

I write a function, run the model, and compare the results. The whole structure is similar to the example in tutorial for using hydropso to calibrate tuwmodel. However, when 10-npart done.

================================================================================
[ Running PSO ... ]

[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
[1] "Done"
Error in Xt.fitness[iter, ] : incorrect number of dimensions

Hello all,
I am unsure of the correct format for inputting the subdaily scale data in the hydromod as it seems hydromod is suitable for daily scale computation (as can be seen in argument "date.format", and there is no option available for hourly data input).

Any help in this regards will be highly appreciated.

With regards,
Ashu

Using R version 3.2.2
hydroPSO version 0.3-4
Windows 7

When trying to optimize a single parameter, I received an error that subsequently did not occur when I optimized 2 parameters. The message was:

[ 4) Computing the goodness of fit ... ]

[NSE= -3.152]

[Iter: 1/10. Particle: 4/4. Finished !. GoF: -3.152331e+00]

                                 |  
                                 |  

Error in X[j, ] <- out[["x.new"]] :
number of items to replace is not a multiple of replacement length

I tracked this error to the position.update.and.boundary.treatment function. I don’t believe it is indexing the vectors correctly in this function. For instance, when optimizing only 1 parameter, my x.new was a vector of length 4, and x.min was a vector of length 1. Byd.min.pos was a vector of length 1, and its value was 2. Then, when boundary.wall == ‘absorbing2011’), x.new[2] <- x.min[2], but x.min is only of length 1, so it creates the error copied above.

As I pointed out, I could probably get around the problem by always optimizing at least 2 parameters, but I’m afraid this error may cause the program to execute in a way it was not intended to.

Here are the parameter settings and function call that I used:

model.drty = sim.wd
setwd(model.drty)
model.FUN.args=list(
model.drty=getwd(),
param.files='ParamFiles.txt',
exe.fname=exe.fname,
verbose=TRUE, #FALSE,
stderr="",
out.FUN=out.FUN,
out.FUN.args=out.FUN.args,
gof.FUN=gof.FUN,
gof.FUN.args=gof.FUN.args,
obs=obs
)

maxit=10
par.type <- 'none' #"parallelWin" #"none"
par.nnodes <- max(detectCores()-2,2) #4
par.pkgs <- c("hydroPSO","lhs","plyr","ggplot2","gridExtra","zoo","parallel","LaF")

drty.out = "PSO_Hays"

drty.out = ""

hydroPSO(fn="hydromod"
, method="spso2011"
, model.FUN="hydromod"
, model.FUN.args=model.FUN.args
, control=list(
drty.in=model.drty
, drty.out=drty.out
, MinMax="max" #"min" #max"
, maxit=maxit #50 #25 #100
, npart=4 #should be more than 4. Can set to NA for the default value of 40
#, c1=2.05
#, c2=2.05
#, use.IW=FALSE
#, use.CF=TRUE
, abstol=0.95 #0.001 #0.95
, reltol=1E-8 #1E-5 #1E-8 #0.0005 #1e-20
#, Xini.type="lhs"
#, Vini.type="lhs2011"
#, best.update="sync"
#, boundary.wall="absorbing2011"
#, topology="random"
#, K=11
#, use.TVlambda=TRUE
#, TVlambda.rng=c(1, 0.5)
#, plot=FALSE
#, write2disk=TRUE
, REPORT=1 #25
, parallel=par.type
, par.nnodes=par.nnodes
, par.pkgs=par.pkgs
)
)

I cannot install hydropso package using conda. Please let me know how I should proceed (I was able to install it on my own computer directly through rstudio using install.packages("hydropso") but I wasn't able to install it on a linux virtual environment on the university server using conda.)

My calibration is running in windows system. It is fast at the beginning, then it gets slow and slow. Is there any way to fix this, such as remove some temporary files during the calibration?

Kangsheng

When I use hydroPSO, I found that I cannot control digit numbers after the decimal point when writing a new value for the calibration parameter. For example, hydroPSO would write one value as 32.045634, while I only want 32.04. How to fix this, Mauricio?

Kangsheng

It might be convenient to include a progress bar to tell the progress. For example I have 10,000 calibration runs, it goes for several days with my current laptop. Built-in progress bar would be a nice to have.

Hello.
I've implemented the VIC hydrologic model, now I want to calibrate the model.
Can I use your package?

Hi Mauricio,

I would like to use parallel function with hydroPSO. I think our machine here is Window based. I read your instruction, it seems that there are not many examples to show how to use such function, or maybe I just missed that part. Please provide some examples on how to use parallel in Window based machine.
Thanks so much
Kangsheng

Can someone help me with vic calibration

In R 4.x "array" was added to the class of matrix objects in addition to the "matrix" class. The checks in the hydroPSO function for vector vs matrix/data.frame throw warnings in R 4.x.

tidymodels/recipes#422
https://cran.r-project.org/doc/manuals/r-release/NEWS.html

00check.log
00install.out.txt
is a recent check log from running under _SP_EVOLUTION_STATUS_=2 without retiring r-spatial packages (suggested in this package) on the library path. This is caused by hydroTSM depending strongly on maptools to use maptools::readShapePoly, which is has been deprecated for years. maptools will retire in October 2023. Follow up on hzambran/hydroTSM#19 (see https://r-spatial.org/r/2022/04/12/evolution.html, https://r-spatial.org/r/2022/12/14/evolution2.html).

The 'multicore' core package is no longer maintained by CRAN, and all its functionality is inherited in the 'parallel' package:
setting the control=list(parallel='multicore') in the function call yields the warning:
[ Package 'multicore' is not installed => parallel='none' ]

To still be able to use this functionality though, simply setting control=list(parallel='parallel') calls the right package