.Rbuildignore

@@ -2,3 +2,14 @@
 ^\.Rproj\.user$
 ^\.Rprofile$
 ^packrat/
+^tests/tmp/
+^\.gitlab-ci.yml$
+^\.regressionignore$
+^\.vignettechunkignore$
+^\.gitlab-ci\.yml$
+^\.vscode$
+^Rplots\.pdf$
+^ci$
+^data-raw$
+^revdep$
+^\.devcontainer$

.devcontainer/devcontainer.json

+{
+	"image": "rocker/geospatial:devel",
+	"customizations": {
+		"vscode": {
+			"extensions": [
+				"eamodio.gitlens",
+				"REditorSupport.r"
+			]
+		}
+	},
+	// Use 'postCreateCommand' to run commands after the container is created.
+	"postCreateCommand": "R -q -e 'install.packages(\"languageserver\");remotes::install_deps(dep = TRUE)'",
+	"postStartCommand": "R -q -e 'devtools::install()'"
+}

.gitignore

-.Rproj.user
+# Specific files for airGR
+packrat/lib*/
+
+# Compiled files
+/src/*.o
+/src/*.so
+/src/*.dll
+/src-*
+
+# Test temporary files
+/tests/tmp/
+*.pdf
+!man/figures/*.pdf
+
+# revdep
+/revdep/
+
+######################################################################################################
+### Generic .gitignore for R (source: https://github.com/github/gitignore/blob/master/R.gitignore) ###
+######################################################################################################
+
+# History files
 .Rhistory
+.Rapp.history
+
+# Session Data files
 .RData
-airGR.Rproj
-packrat/lib*/
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron
+.Rprofile
+
+# pkgdown site
+docs/
+
+# vscode IDE
+.vscode/*
+*.code-workspace
+.history/

.gitlab-ci.yml

+stages:
+  - check
+  - scheduled_tests
+  - revdepcheck
+
+variables:
+  R_LIBS_USER: "$CI_PROJECT_DIR/ci/lib"
+
+default:
+  tags: [docker]
+
+cache:
+  key: "airGR-${R_VERSION}"
+  paths:
+    - $R_LIBS_USER
+
+before_script:
+  - mkdir -p $R_LIBS_USER
+  - echo "R_LIBS='$R_LIBS_USER'" > .Renviron
+  - echo 'options(repos = c(CRAN = "https://packagemanager.posit.co/cran/__linux__/focal/latest"))' > .Rprofile
+  - apt-get update && apt-get install -y libstdc++6
+  - R -q -e 'remotes::install_deps(dependencies = TRUE)'
+
+.R-devel:
+  image: rocker/geospatial:devel
+  variables:
+    R_VERSION: devel
+
+.R-patched:
+  only:
+    refs:
+      - tags
+      - schedules
+  image: rocker/geospatial:latest
+  variables:
+    R_VERSION: patched
+
+.R-oldrel:
+  only:
+    refs:
+      - tags
+      - schedules
+  image: rocker/geospatial:4.1
+  variables:
+    R_VERSION: oldrel
+
+.scheduled_tests:
+  only:
+    refs:
+      - tags
+      - schedules
+      - master
+  stage: scheduled_tests
+  script:
+    - R -q -e 'devtools::install()'
+    - Rscript tests/scheduled_tests/scheduled.R
+    - Rscript tests/scheduled_tests/regression.R stable
+    - R CMD INSTALL .
+    - Rscript tests/scheduled_tests/regression.R dev
+    - Rscript tests/scheduled_tests/regression.R compare
+
+.check:
+  stage: check
+  script:
+  - R -q -e 'tinytex::tlmgr("option repository https://ftp.tu-chemnitz.de/pub/tug/historic/systems/texlive/2021/tlnet-final")'
+  - tlmgr update --self && tlmgr install ec epstopdf-pkg amsmath
+  - R -q -e 'remotes::update_packages("rcmdcheck")'
+  - R -q -e 'rcmdcheck::rcmdcheck(args = c("--as-cran"), error_on = "warning")'
+
+.test_all:
+  stage: check
+  script:
+  - R -q -e 'testthat::test_local()'
+
+benchmark_devel:
+  extends: .R-devel
+  stage: check
+  allow_failure: true
+  script:
+  - R -q -e 'remotes::update_packages("microbenchmark", repos = "http://cran.r-project.org")'
+  - R -q -e 'install.packages("airGR", repos = "http://cran.r-project.org")'
+  - Rscript tests/scheduled_tests/benchmarkRunModel.R
+  - R CMD INSTALL .
+  - Rscript tests/scheduled_tests/benchmarkRunModel.R
+  artifacts:
+    paths:
+      - tests/tmp/benchmark.tsv
+      - tests/tmp/mean_execution_time.tsv
+
+scheduled_tests_patched:
+  extends:
+    - .scheduled_tests
+    - .R-patched
+
+scheduled_tests_devel:
+  extends:
+    - .scheduled_tests
+    - .R-devel
+
+scheduled_tests_oldrel:
+  extends:
+    - .scheduled_tests
+    - .R-oldrel
+
+check_patched:
+  extends:
+    - .check
+    - .R-patched
+
+check_devel:
+  extends:
+    - .check
+    - .R-devel
+
+test_all_patched:
+  extends:
+    - .test_all
+    - .R-patched
+
+test_all_devel:
+  extends:
+    - .test_all
+    - .R-devel
+
+test_all_oldrel:
+  extends:
+    - .test_all
+    - .R-oldrel
+
+revdepcheck_devel:
+  stage: revdepcheck
+  only:
+    refs:
+      - tags
+      - schedule
+      - master
+  extends: .R-devel
+  script:
+    - R -q -e 'remotes::install_github("https://github.com/r-lib/revdepcheck")'
+    - R -q -e 'revdepcheck::revdep_check(timeout = as.difftime(20, units = "mins"))'
+    - R -q -e 'stopifnot(all("+" == sapply(revdepcheck::revdep_summary(), "[[", "status")))'
+    - R -q -e 'if (any(sapply(revdepcheck::revdep_summary(), function(x) {any(x$cmp$change == 1)}))) stop()'
+  artifacts:
+    when: on_failure
+    paths:
+      - revdep/README.md
+      - revdep/problems.md
+      - revdep/failures.md
+      - revdep/cran.md
+      - revdep/checks/*.log

.regressionignore

+# .test-regression.ignore contains the list of topic/variables produces by
+# documentation examples that should be ignore in the regression test
+# The format of this file is: 5 lines of comments followed by one line by
+# ignored variable : [Topic]<SPACE>[Variable] or *<SPACE>[Variable] for every variable whatever the topic
+# Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel

.vignettechunkignore

+# This file is used by the script tests/testthat/test-vignettes which test all
+# chunks including those with `eval=FALSE`
+# It serves to ignore chunks that should not be tested anyway
+# Format: `vignette file name`[space]`id of the chunk`
+V02.1_param_optim.Rmd hydroPSO1
+V02.1_param_optim.Rmd hydroPSO2
+V02.1_param_optim.Rmd resGLOB

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.0.9.55
-Date: 2017-10-26
+Version: 1.7.6.9000
+Date: 2023-10-25
 Authors@R: c(
-  person("Laurent", "Coron", role = c("aut", "trl")),
-  person("Charles", "Perrin", role = c("aut", "ths")),
+  person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
+  person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
+  person("Guillaume", "Thirel", role = c("aut", "ths"), comment = c(ORCID = "0000-0002-1444-1830")),
+  person("David", "Dorchies", role = c("aut"), comment = c(ORCID = "0000-0002-6595-7984")),
+  person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
   person("Claude", "Michel", role = c("aut", "ths")),
-  person("Vazken", "Andréassian", role = c("ctb", "ths")),
-  person("Pierre", "Brigode", role = c("ctb")),
-  person("Olivier", "Delaigue", role = c("cre", "ctb"), email = "airGR@irstea.fr"),
+  person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
+  person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260")),
+  person("Pierre", "Brigode", role = c("ctb"), comment = c(ORCID = "0000-0001-8257-0741")),
   person("Nicolas", "Le Moine", role = c("ctb")),
-  person("Thibaut", "Mathevet", role = c("ctb")),
+  person("Thibaut", "Mathevet", role = c("ctb"), comment = c(ORCID = "0000-0002-4142-4454")),
   person("Safouane", "Mouelhi", role = c("ctb")),
-  person("Ludovic", "Oudin", role = c("ctb")),
+  person("Ludovic", "Oudin", role = c("ctb"), comment = c(ORCID = "0000-0002-3712-0933")),
   person("Raji", "Pushpalatha", role = c("ctb")),
-  person("Guillaume", "Thirel", role = c("ctb")),
   person("Audrey", "Valéry", role = c("ctb"))
   )
-Depends: R (>= 3.0.1)
-Description: Hydrological modelling tools developed at Irstea-Antony (HBAN Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR4J, GR5J, GR6J, GR2M, GR1A), a snow accumulation and melt model (CemaNeige) and the associated functions for their calibration and evaluation. Use help(airGR) for package description.
+Depends: R (>= 3.1.0)
+Imports:
+  graphics,
+  grDevices,
+  stats,
+  utils
+Suggests:
+  knitr, markdown, rmarkdown,
+  caRamel, coda, DEoptim, FME, ggmcmc, Rmalschains,
+  GGally, ggplot2,
+  testthat
+Description: Hydrological modelling tools developed at INRAE-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR5H, GR4J, GR5J, GR6J, GR2M, GR1A) that can be applied either on a lumped or semi-distributed way. A snow accumulation and melt model (CemaNeige) and the associated functions for the calibration and evaluation of models are also included. Use help(airGR) for package description and references.
 License: GPL-2
-URL: https://webgr.irstea.fr/en/airGR/
+URL: https://hydrogr.github.io/airGR/
+BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues
 NeedsCompilation: yes
 Encoding: UTF-8
 VignetteBuilder: knitr
-Suggests: knitr
+RoxygenNote: 7.1.1

NAMESPACE

@@ -8,7 +8,13 @@ useDynLib(airGR, .registration = TRUE)
 #####################################
 ##            S3 methods           ##
 #####################################
-S3method("plot", "OutputsModel")
+S3method('[', InputsModel)
+#S3method('[', OutputsModel) ### to add in version 2.0
+S3method(plot, OutputsModel)
+S3method(SeriesAggreg, data.frame)
+S3method(SeriesAggreg, list)
+S3method(SeriesAggreg, InputsModel)
+S3method(SeriesAggreg, OutputsModel)
 
 
 
@@ -18,8 +24,10 @@ S3method("plot", "OutputsModel")
 export(Calibration)
 export(Calibration_Michel)
 export(CreateCalibOptions)
+export(CreateErrorCrit_GAPX)
 export(CreateIniStates)
 export(CreateInputsCrit)
+export(CreateInputsCrit_Lavenne)
 export(CreateInputsModel)
 export(CreateRunOptions)
 export(DataAltiExtrapolation_Valery)
@@ -28,30 +36,38 @@ export(ErrorCrit_KGE)
 export(ErrorCrit_KGE2)
 export(ErrorCrit_NSE)
 export(ErrorCrit_RMSE)
-export(PEdaily_Oudin)
+export(Imax)
+export(PE_Oudin)
+export(plot.OutputsModel) ### to remove from version 2.0
 export(RunModel)
 export(RunModel_CemaNeige)
+export(RunModel_CemaNeigeGR4H)
+export(RunModel_CemaNeigeGR5H)
 export(RunModel_CemaNeigeGR4J)
 export(RunModel_CemaNeigeGR5J)
 export(RunModel_CemaNeigeGR6J)
 export(RunModel_GR1A)
 export(RunModel_GR2M)
 export(RunModel_GR4H)
+export(RunModel_GR5H)
 export(RunModel_GR4J)
 export(RunModel_GR5J)
 export(RunModel_GR6J)
+export(RunModel_Lag)
 export(SeriesAggreg)
 export(TransfoParam)
 export(TransfoParam_CemaNeige)
+export(TransfoParam_CemaNeigeHyst)
 export(TransfoParam_GR1A)
 export(TransfoParam_GR2M)
 export(TransfoParam_GR4H)
+export(TransfoParam_GR5H)
 export(TransfoParam_GR4J)
 export(TransfoParam_GR5J)
 export(TransfoParam_GR6J)
-export(plot.OutputsModel)
-export(plot_OutputsModel)
-
+export(TransfoParam_Lag)
+#export(.ErrorCrit)
+#export(.FeatModels)
 
 
 #####################################

NEWS.rmd

----
-title: "Release History of the airGR Package"
-output: 
-  html_document:
-    toc: true
-    toc_float: true
-    depth: 3  # upto three depths of headings (specified by #, ### and ###)
-    number_sections: false  ### if you want number sections at each table header
-    theme: united  # many options for theme, this one is my favorite.
-    highlight: tango  # specifies the syntax highlighting style  
-    keep_md: true
----
-
-
-
-
-### 1.0.9.55 Release Notes (2017-10-26) 
-
-#### New features
-
-- An article describing the airGR package has been published. Its reference has been added and will be displayed with <code>citation("airGR")</code>. 
-
-- Added <code>CreateIniStates()</code> function in order to help user to format <code>IniStates</code> argument for <code>CreateRunOptions()</code>.
-
-- Added (<code>Param_Sets_GR4J</code>) dataset. It contains generalist parameter sets for the GR4J model.
-
-
-#### Bug fixes
-
-- Fixed bug in <code>RunModel_GR4H()</code>: in <code>frun_GR4H</code> Fortran subroutine, <code>St(2)</code> is now set to 0 (and not <code>St(1)</code>) when <code>St(2) < 0</code>. 
-
-- Fixed bug in <code>plot.OutputsModel()</code> for the regime plot when the period is less than 1 year.
-
-- Fixed bug in <code>plot.OutputsModel()</code> when there is no common data to plot the cumulative frequency or the correlation QQ.
-
-- Fixed bug in <code>plot.OutputsModel()</code> for y-axis labelling of flows time series when <code>log_scale = TRUE</code> and <code>BasinArea</code> used.
-
-
-#### Deprectated and defunct
-
-- The <code>RunSnowModule</code> argument is now deprecated in <code>CreateRunOptions()</code>.
-
-
-#### Major user-visible changes
-
-- <code>RunModel_GR4J()</code>, <code>RunModel_GR5J()</code> and <code>RunModel_GR6J()</code> (and <code>CemaNeige_GR&#42;J()</code>) now return Ps, Pn and actual exchanges. See the model Fortran codes for more details about the calculation of these variables.
-
-- <code>CreateInputsModel()</code> now returns an error when <code>DatesR</code> contains duplicated values.
-
-- <code>RunModel_GR5J</code> now returns <code>StateEnd</code> in the same order as the other models.
-
-
-#### Minor user-visible changes
-
-- <code>plot.OutputsModel()</code> now returns a warning message when the length of Qobs is different from the length of Qsim. 
-
-- The X1, X3 (and X6) parameters from GR4H, GR4J, GR2M, GR5J (and GR6J) are now set to 1e-2 when they are fixed to lower values. <code>RunModel&#42;()</code> functions now return a warning message. <code>RunModel&#42;()</code> functions now return a warning when X4 < 0.5 and set it to 0.5.
-
-- The commands <code>?L0123001</code>, <code>?L0123002</code> and <code>?L0123003</code> now return the documentation page related to <code>BasinObs</code>.
-
-- Many functions of the package were cleant or slightly modified, with no effect on their outputs. 
-
-- The documentation and help of several functions were improved.
-
-
-#### CRAN-comparibility updates
-
-- "airGR.c" file registers native routines.
-
-____________________________________________________________________________________
-
-
-### 1.0.5.12 Release Notes (2017-01-23) 
-
-
-#### New features
-
-- <code>DataAltiExtrapolation_Valery()</code> and <code>CreateInputsModel()</code> now present a PrecipScale argument which allows rescaling precipitation when it is interpolated on the elevation layers when CemaNeige is used.
-
-
-#### Bug fixes
-
-- Fixed bug in <code>DataAltiExtrapolation_Valery()</code>. The elevation gradients for air temperature returned by <code>CreateInputsModel()</code> are improved.
-
-
-#### User-visible changes
-
-- <code>DataAltiExtrapolation_Valery()</code> has been improved. <code>DataAltiExtrapolation_Valery()</code> now runs faster (and by consequence <code>CreateInputsModel()</code> too, when CemaNeige is used).
-
-____________________________________________________________________________________
-
-
-### 1.0.4 Release Notes (2017-01-18) 
-
-#### New features
-
-- <code>RunModel_CemaNeige()</code>, <code>RunModel_CemaNeigeGR4J()</code>, <code>RunModel_CemaNeigeGR5J()</code> and <code>RunModel_CemaNeigeGR6J()</code> now return air temperature for each elevation layer. 
-
-
-#### Deprectated and defunct
-
-- S3 plot method defined for <code>OutputsModel</code> objects. It means that the <code>plot_OutputsModel()</code> function is deprecated and his use has been replaced by the use of <code>plot.OutputsModel()</code> or <code>plot()</code>.
-
-- In <code>plot.OutputsModel()</code> the <code>PlotChoice</code> argument is deprecated and has been renamed <code>which</code>.
-
-
-#### User-visible changes
-
-- <code>plot.OutputsModel()</code> displays air temperature time series for each layer when <code>CemaNeige</code> is used (argument <code>which = "Temp"</code> or <code>"all"</code>).
-
-____________________________________________________________________________________
-
-
-### 1.0.3 Release Notes (2016-12-09) 
-
-#### New features
-
-- <code>ErrorCrit&#42;()</code> functions gain a warnings argument to replace the verbose action and the verbose argument now prints the criterion value(s).
-
-
-#### Bug fixes
-
-- Fixed bug in <code>CreateCalibOptions()</code> when <code>StartParamList</code> or <code>StartParamDistrib</code> arguments are used.
-
-
-#### User-visible changes
-
-- <code>CreateInputsModel()</code> now returns an error if <code>NLayers <= 0</code> when <code>CemaNeige</code> is used.
-
-- <code>plot_OutputsModel()</code> now displays raw values on the y-axis when the discharge time series is represented with log scale (formerly, log values of discharges were displayed on the y-axis).
-
-
-____________________________________________________________________________________
-
-
-### 1.0.2 Release Notes (2016-11-03) 
-
-
-#### New features
-
-- <code>SeriesAggreg()</code> gains a TimeLag argument that corresponds to a numeric value indicating a time lag (in seconds) for the time series aggregation (useful to aggregate hourly time series to the daily time step for instance).
- In addition, the function now accepts input dates in both POSIXt formats (POSIXct and POSIXlt). The output is in POSIXct format.
- 
-- <code>plot_OutputsModel()</code> gains a <code>log_scale</code> argument in order to plot the flow with a log scale.
-
-- A tutorial is available online on the following link: from http://webgr.irstea.fr/airGR.
- It can also be displayed with the <code>vignette("airGR")</code> command
- 
-
-#### Bug fixes
-
-- The value <code>sort</code> for the <code>transfo</code> argument of <code>CreateInputsCrit()</code> was not taken into account. It is now fixed. 
-
-
-#### Deprectated and defunct
-
-- <code>CreateCalibOptions()</code> loses the OptimParam argument that was redundant with the <code>FixedParam</code> argument. The <code>Calibration_Michel()</code> was modified to take into account this change by using directly <code>FixedParam</code>, but this is transparent to the user. 
-
-- <code>CreateCalibOptions()</code> loses the StartParam argument that was not used. 
-
- 
-#### Major user-visible changes
-
-- The <code>RunModel_GR6J()</code> and <code>RunModel_CemaNeigeGR6J()</code> models were modified back to versions previous to 1.0.1 to prevent from unwanted efficiency criteria deterioration related to the calibration with <code>Calibration_Michel()</code>. 
- The actual model codes were not modified but the <code>TransfoParam_GR6J()</code> and <code>CreateCalibOptions()</code> functions were modified regarding the X5 parameter. 
- It is strongly advised to use airGR 1.0.2 for the <code>RunModel_GR6J()</code> and <code>RunModel_CemaNeigeGR6J()</code> functions if you are using <code>Calibration_Michel()</code>, as they are much more efficient. 
- In case you were using your own calibration algorithm, you will not notice any difference. 
-
-
-#### Minor user-visible changes
-
-- <code>CreateInputsModel()</code> and <code>DataAltiExtrapolation_Valery()</code> functions now allow both POSIXt formats (POSIXct and POSIXlt).
-
-
-____________________________________________________________________________________
-
-
-### 1.0.1 Release Notes (2016-04-21) 
-
-#### Deprectated and defunct
-
-- The <code>Calibration_HBAN()</code> and <code>DataAltiExtrapolation_HBAN()</code> functions have respectively been renamed as <code>Calibration_Michel()</code> and <code>DataAltiExtrapolation_Valery()</code> after the names of their creators.
-
-- The <code>Calibration_optim()</code> function has been removed from the package.
-
-- The silent mode is now defined by the <code>verbose = TRUE</code> argument (formerly <code>quiet = FALSE</code>) in the following functions:
-<code>Calibration()</code>, <code>Calibration_Michel()</code>, <code>CreateInputsModel()</code>, <code>CreateRunOptions()</code>, <code>DataAltiExtrapolation_Valery()</code>, <code>ErrorCrit()</code>, <code>ErrorCrit_KGE()</code>, <code>ErrorCrit_KGE2()</code>, <code>ErrorCrit_NSE()</code>, <code>ErrorCrit_RMSE()</code>, <code>plot_OutputsModel()</code>, <code>SeriesAggreg()</code>.
-
-
-#### Major user-visible changes
-- The GR5J model has been modified: previously, two unit hydrographs were used, now only one is remaining.
- As a consequence, simulations from the GR5J (<code>RunModel_GR5J()</code> function) and CemaNeige (<code>RunModel_CemaNeigeGR5J()</code> function) models will be different.
-
-- An important proportion of the transformations of the parameters have been modified (<code>TransfoParam&#42;()</code> functions). Since this modifies the local search, calibration results will be different .
-
-- The quantiles of the parameters have been recalculated with the new transformations (<code>CreateCalibOptions()</code> function). Since these quantiles constitute the starting point of the calibration algorithm, calibration results will be different.
-
-
-#### Minor user-visible changes
-
-- The Fortran model core codes have been modified:
-	- optimisation of the codes for fastening of computation;
-	- simplification of the internal variables for easier reading and understanding.
-	
-- The list of the contributors and authors is now full.
-
-- The references of the package has been updated; they are returned by the following R-command <code>citation("airGR")</code>.
-
-____________________________________________________________________________________
-
-
-### 0.8.1.2 Release Notes (2015-08-21) 
-
-#### Bug fixes
-
-- Bug fxed in <code>CreateInputsModel()</code> related to the handling of missing values.
-
-- Bug fxed in <code>CreateRunOptions()</code> preventing the correct use of the <code>IniResLevels</code> argument (to manually set the filling rate of the production and routing stores).
-
-
-#### Minor user-visible changes
-
-- Removal of an unnecessary warning when <code>IndPeriod_WarmUp = 0</code>.
-
-
-#### CRAN-comparibility updates
-
-- Modification of namespace file to ensure proper use under linux whithout compilation issues.
-
-
-____________________________________________________________________________________
-
-
-### 0.8.0.2 Release Notes (2015-04-15) 
-
-
-#### New features
-
-- Three new hydrological models: <code>RunModel_GR4H() function for </code> GR4H (hourly), <code>RunModel_GR2M()</code> function for GR2M (monthly) and <code>RunModel_GR1A()</code> function for GR1A (yearly).
-
-- New function <code>SeriesAggreg()</code> to easily aggreg timesteps.
-
-
-#### Bug fixes
-
-- Fixed bug in <code>ErrorCrit_RMSE()</code> which led to incorrect calibration (the criterion was maximised instead of minimised).
-
-
-#### Major user-visible changes
-
-- Update of the functions <code>CreateRunOptions()</code>, <code>CreateCalibOptions()</code> and <code>plot_OutputsModel()</code> to handle the new models.
-
-- Modification of CemaNeige Fortran code to add an update of Gratio after the SnowPack update (no impact on snow simulation).
-
-
-#### Minor user-visible changes
-
-- Improvement of the <code>plot_OutputsModel()</code> function to allow a selection among available plots.
-
-- Minor update in <code>ErrorCrit_KGE()</code> and <code>ErrorCrit_KGE2()</code> to handle case when only one values in not NA.
-
-- Update of the scripts in airGR-advanced-example to match the structures of the <code>BasinData</code> objects.
-
-- Correction of formatting issue in airGR-advanced-example regarding the "List_HypsoData.txt" file.
-
-
-____________________________________________________________________________________
-
-
-### 0.7.4 Release Notes (2014-11-01) 
-
-#### User-visible changes
-
-- Stable version stemming from version 0.7.3.
-
-____________________________________________________________________________________
-
-
-### 0.7.3 Release Notes (2014-XX-XX) 
-
-#### Minor user-visible changes
-
-- Improvements allowing the arrival of new models.
-
-- Improvements of the argument verifications in <code>CreateInputsModel()</code>, <code>CreateRunOptions()</code>, <code>CreateInputsCrit()</code>, <code>CreateCalibOptions()</code>.
-
-- Improvements of all the <code>ErrorCrit()</code> functions to better account for the cases with constant flow values or local zeros.
-
-____________________________________________________________________________________
-
-
-### 0.7.2 Release Notes (2014-07-14) 
-
-#### New features
-
-- New argument in many functions (<code>quiet = TRUE</code> or <code>FALSE</code>) to choose if the warnings should be suppressed or not.
-
-
-#### Bug fixes
-
-- Fixed bug in <code>CreateCalibOptions()</code> to handle models with only one parameter.
-
-
-#### Major user-visible changes
-
-- CemaNeige users must now specify one <code>MeanAnSolidPrecip</code> for each elevation layer. The <code>CreateRunOptions()</code> function is impacted.
-
-- CemaNeige users can now specify the mean elevation of the input series (before it was always considered equal to the catchment median elevation). 
- The impacted functions are <code>CreateInputsModel()</code> and <code>DataAltiExtrapolation_HBAN()</code>.
-
-
-#### Minor user-visible changes
-
-- Improvement of the <code>plot_OutputsModel</code> function (to handle 0 in Qobs and Qsim).
-
-- Improved documentation.
-
-____________________________________________________________________________________
-
-
-### 0.7.1 Release Notes (2014-07-14) 
-
-#### Bug fixes
-
-- Fixed bug in <code>Calibration_HBAN()</code>. The function was not working properly with models having only one parameter.
-
-
-#### Deprecated and defunct
-
-- The <code>CalibrationAlgo_&#42;()</code> functions were renamed into<code>Calibration_&#42;()</code>.
-
-
-#### Major user-visible changes
-
-- New architecture with better format verification procedure (using classes) and simpler setting of default configuration.
-
-- New architecture where the model, calibration and error functions are in the arguments of the functions
- (the exotic use of "generic function" created by the users has been removed).
-
-- Improved documentation and examples.
-
-
-
-____________________________________________________________________________________
-
-
-### 0.7.0 Release Notes (2014-06-06) 
-
-#### Experimental features
-
-- unfinished version used for development purpose.
-
-____________________________________________________________________________________
-
-
-### 0.6.2 Release Notes (2014-02-12) 
-
-#### Bug fixes
-
-- RC11 bug correction: the automatic selection of the warm-up period was not working properly when no data was available from warm-up (i.e. when the user had set the run to start at the very first index).
-
-- RC10 bug correction: the <code>CalibrationAlgo_HBAN()</code> function was not working in the very rare case when the diagonal search was activated and lead to a set outside the authorised range.
-
-- RC9 bug correction: the <code>CalibrationAlgo_HBAN()</code> function was not working properly with models having only one parameter.
-
-- RC8 bug correction of the <code>ModelDefaultIniOptions()</code> function (this bug was introduced in the RC7 and caused an error when <code>IndPeriod_WarmUp = NULL</code>
-
-- RC7 bug correction of the <code>ModelDefaultIniOptions()</code> function (the automatic selection of one year for warm-up was not handling properly missing data).
-
-- RC6 correction of the help files (the description of CemaNeige parameters were inverted).
-
-- RC5 differs from previous releases in the way the data are read and stored (in a list instead of individual vectors).
- The package is similar, only the examples of Main and the files in MyScriptBlocks have changed.
- All basin data are now stored inside a list named <code>BasinData</code>. This will greatly ease the future use of Rdata files (instead of txt files) as storage format for the time series of observation. 
-
-____________________________________________________________________________________
-
- 
-### 0.6.1 Release Notes (2014-02-12) 
-
-
-#### New features
-
-- Additional functions for results plotting (the <code>{zoo}</code> package is required for some of them).
-
-- Add multi-objective calibration using <code>nsga2()</code> (the <code>{mco}</code> package is required).
-
-#### Major user-visible changes
-
-- The definition of the generic function is now made in a much simpler way (e.g. see <code>DefineFunctions_Model()</code> or <code>DefineFunctions_ErrorCrit()</code>).
-
-
-#### Minor user-visible changes
-
-- Code improvements to reduce the computation time.
-
-- Clearer instructions for the adding and modification of a model.
-
-- Improvements of the documentation.
-
-____________________________________________________________________________________
-
- 
-### 0.6.0 Release Notes (2014-02-12) 
-
-#### New features
-
-- The field Multiplier has been added in the ErrorCrit() outputs, to indicate whether the criterion is an error (to minimise) or and efficiency (to maximise).
- This allows to provide real efficiency values in the outputs e.g. NSE[Q] instead of (-1) &times; NSE[Q].
-
-
-#### Deprecated and defunct
-
-- <code>EfficiencyCrit()</code> have been replaced by <code>ErrorCrit()</code> to avoid misunderstanding (by default, the algorithms minimise the error criterion).
-
-____________________________________________________________________________________
-
-
-### 0.5.2 Release Notes (2014-02-05) 
-
-
-#### Deprectated and defunct
-
-- The <code>SelectPer</code> arguments are replaced by <code>IndPeriod</code> to ease understanding.
-
-- The <code>PE</code> arguments are replaced by <code>PotEvap()</code> to ease understanding.
-
-- The <code>Fsol</code> arguments are replaced by <code>FracSolidPrecip</code> to ease understanding.
-
-
-#### Major user-visible changes
-
-- R 2.15 in not supported by default.
-
-- The check that <code>SelectPer_Run()</code> is continuous is now made in the <code>CheckArg()</code> functions.
-
-- Check of the model functioning time step.
-
-- Name of the calibration criterion provided in <code>OutputsAlgo()</code>.
-
-
-#### Minor user-visible changes
-
-- Missing values in Fortran are now -999.999 instead of -9.999.
-
-
-____________________________________________________________________________________
-
-
-### 0.5.1 Release Notes (2014-01-27) 
-
-
-#### New features
-
-- New <code>EfficiencyCrit_NSE_sqrtQ()</code> function to compute NSE criterion on sqrt flows
-
-
-#### Bug fixes
-
-- Incorrect arguments in the call to <code>RunModelAndCrit</code> from <code>CalibrationAlgo_optim_stats</code> and <code>CalibrationAlgo_nlminb_stats</code>.
-
-- <code>CalibrationAlgo_nlminb_stats</code> argument was wrongly defined in <code>DefineFunctions_CalibrationAlgo()</code> (<code>optim</code> instead of <code>nlminb</code>).
-
-- Format checking for <code>RunOptions</code> was incorrectly made in <code>CheckArg()</code> function.
-
-

R/BasinData.R

-#' @name BasinInfo
-#' @docType data
-#' @title Data sample: characteristics of a fictional catchment (L0123001, L0123002 or L0123003)
-#' @description
-#' R-object containing the code, station's name, area and hypsometric curve of the catchment.
-#' @encoding UTF-8
-#' @format
-#' List named 'BasinInfo' containing
-#' \itemize{
-#' \item two strings: catchment's code and station's name
-#' \item one float: catchment's area in km2
-#' \item one numeric vector: catchment's hypsometric curve (min, quantiles 01 to 99 and max) in metres
-#' }
-#' @examples 
-#'    require(airGR)
-#'    data(L0123001)
-#'    str(BasinInfo)
-
-NULL
-
-
-#' @name BasinObs
-#' @docType data
-#' @title Data sample: time series of observations of a fictional catchment (L0123001, L0123002 or L0123003)
-#' @description
-#' R-object containing the times series of precipitation, temperature, potential evapotranspiration and discharges. \cr
-#' Times series for L0123001 or L0123002 are at the daily time-step for use with daily models such as GR4J, GR5J, GR6J, CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J.
-#' Times series for L0123003 are at the hourly time-step for use with hourly models such as GR4H.
-#' @encoding UTF-8
-#' @format
-#' Data frame named 'BasinObs' containing
-#' \itemize{
-#' \item one POSIXlt vector: time series dates in the POSIXlt format
-#' \item five numeric vectors: time series of catchment average precipitation [mm], catchment average air temperature [degC], catchment average potential evapotranspiration [mm], outlet discharge [l/s], outlet discharge [mm]
-#' }
-#' @examples 
-#'    require(airGR)
-#'    data(L0123001)
-#'    str(BasinObs)
-
-NULL
-

R/Calibration.R

-Calibration <- function(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, FUN_CALIB = Calibration_Michel, FUN_TRANSFO = NULL, verbose = TRUE) {
-  return(FUN_CALIB(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, FUN_TRANSFO, verbose = verbose))
+Calibration <- function(InputsModel,
+                        RunOptions,
+                        InputsCrit,
+                        CalibOptions,
+                        FUN_MOD,
+                        FUN_CRIT,                      # deprecated
+                        FUN_CALIB = Calibration_Michel,
+                        FUN_TRANSFO = NULL,
+                        verbose = TRUE,
+                        ...) {
+
+  FUN_MOD   <- match.fun(FUN_MOD)
+
+  if (!missing(FUN_CRIT)) {
+    FUN_CRIT <- match.fun(FUN_CRIT)
+  }
+
+  FUN_CALIB <- match.fun(FUN_CALIB)
+
+  if (!is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
+  }
+
+  return(FUN_CALIB(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit,
+                   CalibOptions = CalibOptions,
+                   FUN_MOD = FUN_MOD, FUN_TRANSFO = FUN_TRANSFO,
+                   verbose = verbose, ...))
 }
 

R/CreateCalibOptions.R

-CreateCalibOptions <-
-  function(FUN_MOD,
-           FUN_CALIB = Calibration_Michel,
-           FUN_TRANSFO = NULL,
-           FixedParam = NULL,
-           SearchRanges = NULL,
-           StartParamList = NULL,
-           StartParamDistrib = NULL) {
-    
-    ObjectClass <- NULL
-    
-    
-    ##check_FUN_MOD
-    BOOL <- FALSE
-    
-    if (identical(FUN_MOD, RunModel_GR4H)) {
-      ObjectClass <- c(ObjectClass, "GR4H")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_GR4J)) {
-      ObjectClass <- c(ObjectClass, "GR4J")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_GR5J)) {
-      ObjectClass <- c(ObjectClass, "GR5J")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_GR6J)) {
-      ObjectClass <- c(ObjectClass, "GR6J")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_GR2M)) {
-      ObjectClass <- c(ObjectClass, "GR2M")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_GR1A)) {
-      ObjectClass <- c(ObjectClass, "GR1A")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeige)) {
-      ObjectClass <- c(ObjectClass, "CemaNeige")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
-      ObjectClass <- c(ObjectClass, "CemaNeigeGR4J")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
-      ObjectClass <- c(ObjectClass, "CemaNeigeGR5J")
-      BOOL <- TRUE
-    }
-    if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-      ObjectClass <- c(ObjectClass, "CemaNeigeGR6J")
-      BOOL <- TRUE
-    }
-    if (!BOOL) {
-      stop("incorrect FUN_MOD for use in CreateCalibOptions \n")
-      return(NULL)
-    }
-    
-    ##check_FUN_CALIB
-    BOOL <- FALSE
-    
-    if (identical(FUN_CALIB, Calibration_Michel)) {
-      ObjectClass <- c(ObjectClass, "HBAN")
-      BOOL <- TRUE
-    }
-    if (!BOOL) {
-      stop("incorrect FUN_CALIB for use in CreateCalibOptions \n")
-      return(NULL)
-      
-    }
-    
-    ##check_FUN_TRANSFO
-    if (is.null(FUN_TRANSFO)) {
-      ##_set_FUN1
-      if (identical(FUN_MOD, RunModel_GR4H)) {
-        FUN1 <- TransfoParam_GR4H
-      }
-      if (identical(FUN_MOD, RunModel_GR4J) |
-          identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
-        FUN1 <- TransfoParam_GR4J
-      }
-      if (identical(FUN_MOD, RunModel_GR5J) |
-          identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
-        FUN1 <- TransfoParam_GR5J
-      }
-      if (identical(FUN_MOD, RunModel_GR6J) |
-          identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-        FUN1 <- TransfoParam_GR6J
-      }
-      if (identical(FUN_MOD, RunModel_GR2M)) {
-        FUN1 <- TransfoParam_GR2M
-      }
-      if (identical(FUN_MOD, RunModel_GR1A)) {
-        FUN1 <- TransfoParam_GR1A
-      }
-      if (identical(FUN_MOD, RunModel_CemaNeige)) {
-        FUN1 <- TransfoParam_CemaNeige
-      }
-      if (is.null(FUN1)) {
-        stop("FUN1 was not found \n")
-        return(NULL)
-      }
-      ##_set_FUN2
-      FUN2 <- TransfoParam_CemaNeige
-      
-      ##_set_FUN_TRANSFO
-      if (sum(ObjectClass %in% c("GR4H", "GR4J", "GR5J", "GR6J", "GR2M", "GR1A", "CemaNeige")) > 0) {
-        FUN_TRANSFO <- FUN1
-      } else {
-        FUN_TRANSFO <- function(ParamIn, Direction) {
-          Bool <- is.matrix(ParamIn)
-          if (Bool == FALSE) {
-            ParamIn <- rbind(ParamIn)
-          }
-          ParamOut <- NA * ParamIn
-          NParam   <- ncol(ParamIn)
-          if (NParam <= 3) {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
-          } else {
-            ParamOut[, 1:(NParam - 2)] <- FUN1(ParamIn[, 1:(NParam - 2)], Direction)
-          }
-          ParamOut[, (NParam - 1):NParam] <- FUN2(ParamIn[, (NParam - 1):NParam], Direction)
-          if (Bool == FALSE) {
-            ParamOut <- ParamOut[1, ]
-          }
-          return(ParamOut)
-        }
-      }
-    }
-    if (is.null(FUN_TRANSFO)) {
-      stop("FUN_TRANSFO was not found \n")
-      return(NULL)
-    }
-    
-    ##NParam
-    if ("GR4H" %in% ObjectClass) {
-      NParam <- 4
-    }
-    if ("GR4J" %in% ObjectClass) {
-      NParam <- 4
-    }
-    if ("GR5J" %in% ObjectClass) {
-      NParam <- 5
-    }
-    if ("GR6J" %in% ObjectClass) {
-      NParam <- 6
-    }
-    if ("GR2M" %in% ObjectClass) {
-      NParam <- 2
-    }
-    if ("GR1A" %in% ObjectClass) {
-      NParam <- 1
-    }
-    if ("CemaNeige" %in% ObjectClass) {
-      NParam <- 2
-    }
-    if ("CemaNeigeGR4J" %in% ObjectClass) {
-      NParam <- 6
-    }
-    if ("CemaNeigeGR5J" %in% ObjectClass) {
-      NParam <- 7
-    }
-    if ("CemaNeigeGR6J" %in% ObjectClass) {
-      NParam <- 8
-    }
-    
-    ##check_FixedParam
-    if (is.null(FixedParam)) {
-      FixedParam <- rep(NA, NParam)
-    } else {
-      if (!is.vector(FixedParam)) {
-        stop("FixedParam must be a vector \n")
-        return(NULL)
-      }
-      if (length(FixedParam) != NParam) {
-        stop("Incompatibility between FixedParam length and FUN_MOD \n")
-        return(NULL)
-      }
-    }
-    
-    ##check_SearchRanges
-    if (is.null(SearchRanges)) {
-      ParamT <-  matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
-                        ncol = NParam, byrow = TRUE)
-      
-      SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
-      
-    } else {
-      if (!is.matrix(SearchRanges)) {
-        stop("SearchRanges must be a matrix \n")
-        return(NULL)
-      }
-      if (!is.numeric(SearchRanges)) {
-        stop("SearchRanges must be a matrix of numeric values \n")
-        return(NULL)
-      }
-      if (sum(is.na(SearchRanges)) != 0) {
-        stop("SearchRanges must not include NA values \n")
-        return(NULL)
-      }
-      if (nrow(SearchRanges) != 2) {
-        stop("SearchRanges must have 2 rows \n")
-        return(NULL)
-      }
-      if (ncol(SearchRanges) != NParam) {
-        stop("Incompatibility between SearchRanges ncol and FUN_MOD \n")
-        return(NULL)
-      }
-    }
-    
-    ##check_StartParamList_and_StartParamDistrib__default_values
-    if (("HBAN"  %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
-      if ("GR4H" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
-                           +5.58, -0.85, +4.74, -9.47,
-                           +6.01, -0.50, +5.14, -8.87), ncol = NParam,  byrow = TRUE)
-      }
-      if ("GR4J" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
-                           +5.51, -0.61, +3.74, -8.51,
-                           +6.07, -0.02, +4.42, -8.06),  ncol = NParam, byrow = TRUE)
-      }
-      if ("GR5J" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
-                           +5.55, -0.46, +3.75, -9.09, -4.69,
-                           +6.10, -0.11, +4.43, -8.60, -0.66), ncol = NParam, byrow = TRUE)
-        
-      }
-      if ("GR6J" %in% ObjectClass) {
-        ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
-                           +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
-                           +4.50, +0.50, +5.00, -8.10, +1.10, +5.00), ncol = NParam, byrow = TRUE)
-      }
-      if ("GR2M" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.03, -7.15,
-                           +5.22, -6.74,
-                           +5.85, -6.37), ncol = NParam, byrow = TRUE)
-      }
-      if ("GR1A" %in% ObjectClass) {
-        ParamT <- matrix(c(-1.69,
-                           -0.38,
-                           +1.39), ncol = NParam, byrow = TRUE)
-      }
-      if ("CemaNeige" %in% ObjectClass) {
-        ParamT <- matrix(c(-9.96, +6.63,
-                           -9.14, +6.90,
-                           +4.10, +7.21), ncol = NParam, byrow = TRUE)
-      }
-      if ("CemaNeigeGR4J" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
-                           +5.51, -0.61, +3.74, -8.51, -9.14, +6.90,
-                           +6.07, -0.02, +4.42, -8.06, +4.10, +7.21), ncol = NParam, byrow = TRUE)
-      }
-      if ("CemaNeigeGR5J" %in% ObjectClass) {
-        ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45, -9.96, +6.63,
-                           +5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
-                           +6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21), ncol = NParam, byrow = TRUE)
-      }
-      if ("CemaNeigeGR6J" %in% ObjectClass) {
-        ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -9.96, +6.63,
-                           +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
-                           +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = NParam, byrow = TRUE)
-      }
-      
-      StartParamList    <- NULL
-      StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
-      
-    }
-    
-    ##check_StartParamList_and_StartParamDistrib__format
-    if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
-      if (!is.matrix(StartParamList)) {
-        stop("StartParamList must be a matrix \n")
-        return(NULL)
-      }
-      if (!is.numeric(StartParamList)) {
-        stop("StartParamList must be a matrix of numeric values \n")
-        return(NULL)
-      }
-      if (sum(is.na(StartParamList)) != 0) {
-        stop("StartParamList must not include NA values \n")
-        return(NULL)
-      }
-      if (ncol(StartParamList) != NParam) {
-        stop("Incompatibility between StartParamList ncol and FUN_MOD \n")
-        return(NULL)
-      }
-    }
-    if ("HBAN" %in% ObjectClass & !is.null(StartParamDistrib)) {
-      if (!is.matrix(StartParamDistrib)) {
-        stop("StartParamDistrib must be a matrix \n")
-        return(NULL)
-      }
-      if (!is.numeric(StartParamDistrib[1, ])) {
-        stop("StartParamDistrib must be a matrix of numeric values \n")
-        return(NULL)
-      }
-      if (sum(is.na(StartParamDistrib[1, ])) != 0) {
-        stop("StartParamDistrib must not include NA values on the first line \n")
-        return(NULL)
-      }
-      if (ncol(StartParamDistrib) != NParam) {
-        stop("Incompatibility between StartParamDistrib ncol and FUN_MOD \n")
-        return(NULL)
-      }
-    }
-    
-    
-    ##Create_CalibOptions
-    CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges)
-    
-    if (!is.null(StartParamList)) {
-      CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
-    }
-    if (!is.null(StartParamDistrib)) {
-      CalibOptions <- c(CalibOptions, list(StartParamDistrib = StartParamDistrib))
-    }
-    class(CalibOptions) <- c("CalibOptions", ObjectClass)
-    
-    return(CalibOptions)
-    
+CreateCalibOptions <- function(FUN_MOD,
+                               FUN_CALIB = Calibration_Michel,
+                               FUN_TRANSFO = NULL,
+                               IsHyst = FALSE,
+                               IsSD = FALSE,
+                               FixedParam = NULL,
+                               SearchRanges = NULL,
+                               StartParamList = NULL,
+                               StartParamDistrib = NULL) {
+
+  ObjectClass <- NULL
+
+  FUN_MOD     <- match.fun(FUN_MOD)
+  FUN_CALIB   <- match.fun(FUN_CALIB)
+  if (!is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
   }
+  if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+    stop("'IsHyst' must be a logical of length 1")
+  }
+  if (!is.logical(IsSD) | length(IsSD) != 1L) {
+    stop("'IsSD' must be a logical of length 1")
+  }
+
+  ## check FUN_MOD
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD)
+  FeatFUN_MOD$IsHyst <- IsHyst
+  FeatFUN_MOD$IsSD <- IsSD
+  ObjectClass <- FeatFUN_MOD$Class
+
+  if (identical(FUN_MOD, RunModel_Lag) && IsSD) {
+      stop("RunModel_Lag should not be used with 'IsSD=TRUE'")
+  }
+  if (IsHyst) {
+    ObjectClass <- c(ObjectClass, "hysteresis")
+  }
+  if (IsSD) {
+    ObjectClass <- c(ObjectClass, "SD")
+  }
+
+  ## check FUN_CALIB
+  BOOL <- FALSE
+
+  if (identical(FUN_CALIB, Calibration_Michel)) {
+    ObjectClass <- c(ObjectClass, "HBAN")
+    BOOL <- TRUE
+  }
+  if (!BOOL) {
+    stop("incorrect 'FUN_CALIB' for use in 'CreateCalibOptions'")
+    return(NULL)
+
+  }
+
+  ## check FUN_TRANSFO
+  if (is.null(FUN_TRANSFO)) {
+    FUN_TRANSFO <- .FunTransfo(FeatFUN_MOD)
+  }
+
+  ## NParam
+  NParam <- FeatFUN_MOD$NbParam
+
+  if (IsHyst) {
+    NParam <- NParam + 2
+  }
+  if (IsSD) {
+    NParam <- NParam + 1
+  }
+
+  ## check FixedParam
+  if (is.null(FixedParam)) {
+    FixedParam <- rep(NA, NParam)
+  } else {
+    if (!is.vector(FixedParam)) {
+      stop("FixedParam must be a vector")
+    }
+    if (length(FixedParam) != NParam) {
+      stop("Incompatibility between 'FixedParam' length and 'FUN_MOD'")
+    }
+    if (all(!is.na(FixedParam))) {
+      stop("At least one parameter must be not set (NA)")
+    }
+    if (all(is.na(FixedParam))) {
+      warning("You have not set any parameter in 'FixedParam'")
+    }
+  }
+
+  ## check SearchRanges
+  if (is.null(SearchRanges)) {
+    ParamT <- matrix(c(rep(-9.99, NParam), rep(+9.99, NParam)),
+                     ncol = NParam, byrow = TRUE)
+    SearchRanges <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
+
+  } else {
+    if (!is.matrix(SearchRanges)) {
+      stop("'SearchRanges' must be a matrix")
+    }
+    if (!is.numeric(SearchRanges)) {
+      stop("'SearchRanges' must be a matrix of numeric values")
+    }
+    if (sum(is.na(SearchRanges)) != 0) {
+      stop("'SearchRanges' must not include NA values")
+    }
+    if (nrow(SearchRanges) != 2) {
+      stop("'SearchRanges' must have 2 rows")
+    }
+    if (ncol(SearchRanges) != NParam) {
+      stop("Incompatibility between 'SearchRanges' ncol and 'FUN_MOD'")
+    }
+  }
+
+  ## check StartParamList and StartParamDistrib default values
+  if (("HBAN"  %in% ObjectClass & is.null(StartParamList) & is.null(StartParamDistrib))) {
+    if ("GR4H" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69,
+                         +5.58, -0.85, +4.74, -9.47,
+                         +6.01, -0.50, +5.14, -8.87), ncol = 4, byrow = TRUE)
+    }
+    if (("GR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34,
+                         +3.74, -0.41, +4.78, -8.94, -3.33,
+                         +4.29, +0.16, +5.39, -7.39, +3.33), ncol = 5, byrow = TRUE)
+    }
+    if (("GR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49,
+                         +3.62, -0.19, +4.80, -9.00, -6.31,
+                         +4.01, -0.04, +5.43, -7.53, -5.33), ncol = 5, byrow = TRUE)
+    }
+    if ("GR4J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05,
+                         +5.51, -0.61, +3.74, -8.51,
+                         +6.07, -0.02, +4.42, -8.06), ncol = 4, byrow = TRUE)
+    }
+    if ("GR5J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45,
+                         +5.55, -0.46, +3.75, -9.09, -4.69,
+                         +6.10, -0.11, +4.43, -8.60, -0.66), ncol = 5, byrow = TRUE)
+
+    }
+    if ("GR6J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
+                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
+                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00), ncol = 6, byrow = TRUE)
+    }
+    if ("GR2M" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.03, -7.15,
+                         +5.22, -6.74,
+                         +5.85, -6.37), ncol = 2, byrow = TRUE)
+    }
+    if ("GR1A" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(-1.69,
+                         -0.38,
+                         +1.39), ncol = 1, byrow = TRUE)
+    }
+
+
+    if ("CemaNeige" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(-9.96, +6.63,
+                         -9.14, +6.90,
+                         +4.10, +7.21), ncol = 2, byrow = TRUE)
+    }
+    if ("CemaNeigeGR4H" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.12, -1.18, +4.34, -9.69, -9.96, +6.63,
+                         +5.58, -0.85, +4.74, -9.47, -9.14, +6.90,
+                         +6.01, -0.50, +5.14, -8.87, +4.10, +7.21), ncol = 6, byrow = TRUE)
+    }
+    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & ("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.46, -1.25, +4.04, -9.53, -9.34, -9.96, +6.63,
+                         +3.74, -0.41, +4.78, -8.94, -3.33, -9.14, +6.90,
+                         +4.29, +0.16, +5.39, -7.39, +3.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if (("CemaNeigeGR5H" == FeatFUN_MOD$CodeMod) & !("interception" %in% ObjectClass)) {
+      ParamT <- matrix(c(+3.28, -0.39, +4.14, -9.54, -7.49, -9.96, +6.63,
+                         +3.62, -0.19, +4.80, -9.00, -6.31, -9.14, +6.90,
+                         +4.01, -0.04, +5.43, -7.53, -5.33, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if ("CemaNeigeGR4J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.13, -1.60, +3.03, -9.05, -9.96, +6.63,
+                         +5.51, -0.61, +3.74, -8.51, -9.14, +6.90,
+                         +6.07, -0.02, +4.42, -8.06, +4.10, +7.21), ncol = 6, byrow = TRUE)
+    }
+    if ("CemaNeigeGR5J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+5.17, -1.13, +3.08, -9.37, -7.45, -9.96, +6.63,
+                         +5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
+                         +6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21), ncol = 7, byrow = TRUE)
+    }
+    if ("CemaNeigeGR6J" == FeatFUN_MOD$CodeMod) {
+      ParamT <- matrix(c(+3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -9.96, +6.63,
+                         +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -9.14, +6.90,
+                         +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.10, +7.21), ncol = 8, byrow = TRUE)
+    }
+
+    if (IsHyst) {
+      ParamTHyst <- matrix(c(-7.00, -7.00,
+                             -0.00, -0.00,
+                             +7.00, +7.00), ncol = 2, byrow = TRUE)
+      ParamT <- cbind(ParamT, ParamTHyst)
+    }
+    if (IsSD) {
+      ParamTSD <- matrix(c(-8.75,
+                           -7.50,
+                           -5.00), ncol = 1, byrow = TRUE)
+      ParamT <- cbind(ParamTSD, ParamT)
+    }
+
+    StartParamList    <- NULL
+    StartParamDistrib <- TransfoParam(ParamIn = ParamT, Direction = "TR", FUN_TRANSFO = FUN_TRANSFO)
+
+  }
+
+  ## check StartParamList and StartParamDistrib format
+  if ("HBAN" %in% ObjectClass & !is.null(StartParamList)) {
+    if (!is.matrix(StartParamList)) {
+      stop("'StartParamList' must be a matrix")
+    }
+    if (!is.numeric(StartParamList)) {
+      stop("'StartParamList' must be a matrix of numeric values")
+    }
+    if (sum(is.na(StartParamList)) != 0) {
+      stop("'StartParamList' must not include NA values")
+    }
+    if (ncol(StartParamList) != NParam) {
+      stop("Incompatibility between 'StartParamList' ncol and 'FUN_MOD'")
+    }
+  }
+  if ("HBAN" %in% ObjectClass & !is.null(StartParamDistrib)) {
+    if (!is.matrix(StartParamDistrib)) {
+      stop("'StartParamDistrib' must be a matrix")
+    }
+    if (!is.numeric(StartParamDistrib[1, ])) {
+      stop("'StartParamDistrib' must be a matrix of numeric values")
+    }
+    if (sum(is.na(StartParamDistrib[1, ])) != 0) {
+      stop("'StartParamDistrib' must not include NA values on the first line")
+    }
+    if (ncol(StartParamDistrib) != NParam) {
+      stop("Incompatibility between 'StartParamDistrib' ncol and 'FUN_MOD'")
+    }
+  }
+
+
+  ##Create_CalibOptions
+  CalibOptions <- list(FixedParam = FixedParam, SearchRanges = SearchRanges, FUN_TRANSFO = FUN_TRANSFO)
+
+  if (!is.null(StartParamList)) {
+    CalibOptions <- c(CalibOptions, list(StartParamList = StartParamList))
+  }
+  if (!is.null(StartParamDistrib)) {
+    CalibOptions <- c(CalibOptions, list(StartParamDistrib = StartParamDistrib))
+  }
+  class(CalibOptions) <- c("CalibOptions", ObjectClass)
+
+  return(CalibOptions)
+
+}

R/CreateErrorCrit_GAPX.R

+CreateErrorCrit_GAPX <- function(FUN_TRANSFO) {
+
+  FUN_CRIT <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
+    ## Arguments check
+    if (!inherits(OutputsModel, "OutputsModel")) {
+      stop("'OutputsModel' must be of class 'OutputsModel'")
+    }
+
+    OutputsModel$RunOptions$ParamT <- FUN_TRANSFO(OutputsModel$RunOptions$Param, "RT")
+
+    EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "GAPX", OutputsModel = OutputsModel, warnings = warnings)
+
+    CritValue <- NA
+
+    if (EC$CritCompute) {
+
+      ParamApr <- EC$VarObs[!EC$TS_ignore]
+      ParamOpt <- EC$VarSim[!EC$TS_ignore]
+
+      ## ErrorCrit
+      Crit <- 1 - sum(((ParamApr - ParamOpt) / 20)^2)^0.5
+
+      if (is.numeric(Crit) & is.finite(Crit)) {
+        CritValue <- Crit
+      }
+
+      ## Verbose
+      if (verbose) {
+        message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
+      }
+    }
+
+    ## Output
+    OutputsCrit <- list(CritValue       = CritValue,
+                        CritName        = EC$CritName,
+                        CritBestValue   = EC$CritBestValue,
+                        Multiplier      = EC$Multiplier,
+                        Ind_notcomputed = EC$Ind_TS_ignore)
+
+    class(OutputsCrit) <- c("GAPX", "ErrorCrit")
+    return(OutputsCrit)
+  }
+
+  class(FUN_CRIT) <- c("FUN_CRIT", class(FUN_CRIT))
+
+  return(FUN_CRIT)
+}

R/CreateIniStates.R

-CreateIniStates <- function(FUN_MOD, InputsModel,
-                            ProdStore = 350, RoutStore = 90, ExpStore = NULL,
+CreateIniStates <- function(FUN_MOD, InputsModel, IsHyst = FALSE, IsIntStore = FALSE,
+                            ProdStore = 350, RoutStore = 90, ExpStore = NULL, IntStore = NULL,
                             UH1 = NULL, UH2 = NULL,
                             GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                            GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
+                            SD = NULL,
                             verbose = TRUE) {
-  
-  
+
+
   ObjectClass <- NULL
-  
+
   UH1n <- 20L
   UH2n <- UH1n * 2L
-  
-  
-  ## check FUN_MOD
-  BOOL <- FALSE
-  if (identical(FUN_MOD, RunModel_GR4H)) {
-    ObjectClass <- c(ObjectClass, "GR", "hourly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR4J) |
-      identical(FUN_MOD, RunModel_GR5J) |
-      identical(FUN_MOD, RunModel_GR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "daily")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR2M)) {
-    ObjectClass <- c(ObjectClass, "GR", "monthly")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_GR1A)) {
-    stop("'RunModel_GR1A' does not require 'IniStates' object")
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeige)) {
-    ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
-    BOOL <- TRUE
-  }
-  if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
-      identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
-      identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
-    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
-    BOOL <- TRUE
+
+  FUN_MOD <- match.fun(FUN_MOD)
+
+  FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = InputsModel$DatesR)
+  ObjectClass <- FeatFUN_MOD$Class
+
+  if (!"CemaNeige" %in% ObjectClass & IsHyst) {
+    stop("'IsHyst' cannot be TRUE if CemaNeige is not used in 'FUN_MOD'")
   }
-  if (!BOOL) {
-    stop("Incorrect 'FUN_MOD' for use in 'CreateIniStates'")
-    return(NULL)
+  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & IsIntStore) {
+    stop("'IsIntStore' cannot be TRUE if GR5H is not used in 'FUN_MOD'")
   }
-  
+
   ## check InputsModel
   if (!inherits(InputsModel, "InputsModel")) {
     stop("'InputsModel' must be of class 'InputsModel'")
-    return(NULL)
   }
   if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
     stop("'InputsModel' must be of class 'GR'")
-    return(NULL)
   }
   if ("CemaNeige" %in% ObjectClass &
       !inherits(InputsModel, "CemaNeige")) {
     stop("'InputsModel' must be of class 'CemaNeige'")
-    return(NULL)
   }
-  
-  
+
+
   ## check states
   if (any(eTGCemaNeigeLayers > 0)) {
     stop("Positive values are not allowed for 'eTGCemaNeigeLayers'")
-  }  
-  
+  }
+
   if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
     if (is.null(ExpStore)) {
       stop("'RunModel_*GR6J' need an 'ExpStore' value")
-      return(NULL)
     }
   } else if (!is.null(ExpStore)) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", as.character(substitute(FUN_MOD))))
+      warning(sprintf("'%s' does not require 'ExpStore'. Value set to NA", FeatFUN_MOD$NameFunMod))
     }
     ExpStore <- Inf
   }
-  
+
   if (identical(FUN_MOD, RunModel_GR2M)) {
     if (!is.null(UH1)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
       UH1 <- rep(Inf, UH1n)
     }
     if (!is.null(UH2)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
       UH2 <- rep(Inf, UH2n)
     }
   }
-  
+
   if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1)) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'UH1'. Values set to NA", as.character(substitute(FUN_MOD))))
+      warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
     }
     UH1 <- rep(Inf, UH1n)
   }
- 
+  if (!(identical(FUN_MOD, RunModel_GR5H) | identical(FUN_MOD, RunModel_CemaNeigeGR5H)) & !is.null(IntStore)) {
+    if (verbose) {
+      warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
+    }
+    IntStore <- Inf
+  }
+
   if ("CemaNeige" %in% ObjectClass & ! "GR" %in% ObjectClass) {
     if (!is.null(ProdStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'ProdStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     ProdStore <- Inf
     if (!is.null(RoutStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'RoutStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     RoutStore <- Inf
     if (!is.null(ExpStore)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'ExpStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     ExpStore <- Inf
+    if (!is.null(IntStore)) {
+      if (verbose) {
+        warning(sprintf("'%s' does not require 'IntStore'. Values set to NA", FeatFUN_MOD$NameFunMod))
+      }
+    }
+    IntStore <- Inf
     if (!is.null(UH1)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'UH1'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     UH1 <- rep(Inf, UH1n)
     if (!is.null(UH2)) {
       if (verbose) {
-        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", as.character(substitute(FUN_MOD))))
+        warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
       }
     }
     UH2 <- rep(Inf, UH2n)
   }
-  if("CemaNeige" %in% ObjectClass &
+  if (IsIntStore & is.null(IntStore)) {
+      stop(sprintf("'%s' need values for 'IntStore'", FeatFUN_MOD$NameFunMod))
+  }
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
      (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
-    stop("'RunModel_CemaNeigeGR*' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'")
-    return(NULL)
+      stop(sprintf("'%s' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
+  }
+  if ("CemaNeige" %in% ObjectClass & IsHyst &
+     (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers) |
+      is.null(GthrCemaNeigeLayers) | is.null(GlocmaxCemaNeigeLayers))) {
+    stop(sprintf("'%s' need values for 'GCemaNeigeLayers', 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'", FeatFUN_MOD$NameFunMod))
+  }
+  if ("CemaNeige" %in% ObjectClass & !IsHyst &
+     (!is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
+    if (verbose) {
+      warning(sprintf("'%s' does not require 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
+    }
+    GthrCemaNeigeLayers    <- Inf
+    GlocmaxCemaNeigeLayers <- Inf
   }
-  if(!"CemaNeige" %in% ObjectClass &
-     (!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers))) {
+  if (!"CemaNeige" %in% ObjectClass &
+     (!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
     if (verbose) {
-      warning(sprintf("'%s' does not require 'GCemaNeigeLayers' and 'GCemaNeigeLayers'. Values set to NA", as.character(substitute(FUN_MOD))))
+      warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
     }
-    GCemaNeigeLayers   <- Inf
-    eTGCemaNeigeLayers <- Inf
+    GCemaNeigeLayers       <- Inf
+    eTGCemaNeigeLayers     <- Inf
+    GthrCemaNeigeLayers    <- Inf
+    GlocmaxCemaNeigeLayers <- Inf
   }
-  
-  
+
+
   ## set states
-  if("CemaNeige" %in% ObjectClass) {
+  if ("CemaNeige" %in% ObjectClass) {
     NLayers <- length(InputsModel$LayerPrecip)
   } else {
     NLayers <- 1
   }
-  
-  
+
+
   ## manage NULL values
   if (is.null(ExpStore)) {
-    ExpStore <- Inf 
+    ExpStore <- Inf
+  }
+  if (is.null(IntStore)) {
+    IntStore <- Inf
   }
   if (is.null(UH1)) {
     if ("hourly"  %in% ObjectClass) {
@@ -165,7 +172,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
     } else {
       k <- 1
     }
-    UH1 <- rep(Inf, 20 * k)
+    UH1 <- rep(Inf, UH1n * k)
   }
   if (is.null(UH2)) {
     if ("hourly"  %in% ObjectClass) {
@@ -173,7 +180,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
     } else {
       k <- 1
     }
-    UH2 <- rep(Inf, 40 * k)
+    UH2 <- rep(Inf, UH2n * k)
   }
   if (is.null(GCemaNeigeLayers)) {
     GCemaNeigeLayers <- rep(Inf, NLayers)
@@ -181,19 +188,29 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
   if (is.null(eTGCemaNeigeLayers)) {
     eTGCemaNeigeLayers <- rep(Inf, NLayers)
   }
-  
-  
+  if (is.null(GthrCemaNeigeLayers)) {
+    GthrCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  if (any(is.infinite(GthrCemaNeigeLayers))) {
+    GthrCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  if (is.null(GlocmaxCemaNeigeLayers)) {
+    GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  if (any(is.infinite(GlocmaxCemaNeigeLayers))) {
+    GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+
   # check negative values
-  if (any(ProdStore < 0) | any(RoutStore < 0) |
+  if (any(ProdStore < 0) | any(RoutStore < 0) | any(IntStore < 0) |
       any(UH1 < 0) | any(UH2 < 0) |
       any(GCemaNeigeLayers < 0)) {
-    stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
+    stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'IntStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
   }
-  
-  
+
+
   ## check length
   if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
-    print(ProdStore)
     stop("'ProdStore' must be numeric of length one")
   }
   if (!is.numeric(RoutStore) || length(RoutStore) != 1L) {
@@ -202,16 +219,19 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
   if (!is.numeric(ExpStore) || length(ExpStore) != 1L) {
     stop("'ExpStore' must be numeric of length one")
   }
-  if ( "hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 480L)) {
+  if (!is.numeric(IntStore) || length(IntStore) != 1L) {
+    stop("'IntStore' must be numeric of length one")
+  }
+  if ( "hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != UH1n * 24)) {
     stop(sprintf("'UH1' must be numeric of length 480 (%i * 24)", UH1n))
   }
-  if (!"hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 20L)) {
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != UH1n)) {
     stop(sprintf("'UH1' must be numeric of length %i", UH1n))
   }
-  if ( "hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 960L)) {
+  if ( "hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != UH2n * 24)) {
     stop(sprintf("'UH2' must be numeric of length 960 (%i * 24)", UH2n))
   }
-  if (!"hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 40L)) {
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != UH2n)) {
     stop(sprintf("'UH2' must be numeric of length %i (2 * %i)", UH2n, UH1n))
   }
   if (!is.numeric(GCemaNeigeLayers) || length(GCemaNeigeLayers) != NLayers) {
@@ -220,18 +240,54 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
   if (!is.numeric(eTGCemaNeigeLayers) || length(eTGCemaNeigeLayers) != NLayers) {
     stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
   }
-  
+  if (IsHyst) {
+    if (!is.numeric(GthrCemaNeigeLayers) || length(GthrCemaNeigeLayers) != NLayers) {
+      stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
+    }
+    if (!is.numeric(GlocmaxCemaNeigeLayers) || length(GlocmaxCemaNeigeLayers) != NLayers) {
+      stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
+    }
+  }
+
+  # SD model state handling
+  if (!is.null(SD)) {
+    if (!inherits(InputsModel, "SD")) {
+      stop("'SD' argument provided and 'InputsModel' is not of class 'SD'")
+    }
+    if (!is.list(SD)) {
+      stop("'SD' argument must be a list")
+    }
+    lapply(SD, function(x) {
+      if (!is.numeric(x)) stop("Each item of 'SD' list argument must be numeric")
+    })
+    if (length(SD) != length(InputsModel$LengthHydro)) {
+      stop("Number of items of 'SD' list argument must be the same as the number of upstream connections",
+           sprintf(" (%i required, found %i)", length(InputsModel$LengthHydro), length(SD)))
+    }
+  }
 
   ## format output
-  IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore),
+  IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore, Int = IntStore),
                     UH = list(UH1 = UH1, UH2 = UH2),
-                    CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers))
+                    CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers,
+                                           Gthr = GthrCemaNeigeLayers, Glocmax = GlocmaxCemaNeigeLayers))
   IniStatesNA <- unlist(IniStates)
   IniStatesNA[is.infinite(IniStatesNA)] <- NA
   IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
-  
+
+  if (!is.null(SD)) {
+    IniStatesNA$SD <- SD
+  }
+
   class(IniStatesNA) <- c("IniStates", ObjectClass)
+  if (IsHyst) {
+    class(IniStatesNA) <- c(class(IniStatesNA), "hysteresis")
+  }
+  if (IsIntStore) {
+    class(IniStatesNA) <- c(class(IniStatesNA), "interception")
+  }
+
   return(IniStatesNA)
-  
-  
+
+
 }

R/CreateInputsCrit.R

-CreateInputsCrit <-
-  function(FUN_CRIT,
-           InputsModel,
-           RunOptions,
-           Qobs,
-           BoolCrit = NULL,
-           transfo = "",
-           Ind_zeroes = NULL,
-           epsilon = NULL) {
-    ObjectClass <- NULL
-    
-    
-    ##check_FUN_CRIT
-    BOOL <- FALSE
-    
-    if (identical(FUN_CRIT, ErrorCrit_NSE) | identical(FUN_CRIT, ErrorCrit_KGE) |
-        identical(FUN_CRIT, ErrorCrit_KGE2) | identical(FUN_CRIT, ErrorCrit_RMSE)) {
-      BOOL <- TRUE
-    }
-    if (!BOOL) {
-      stop("incorrect FUN_CRIT for use in CreateInputsCrit \n")
-      return(NULL)
-    }
-    
-    ##check_arguments
-    if (inherits(InputsModel, "InputsModel") == FALSE) {
-      stop("InputsModel must be of class 'InputsModel' \n")
-      return(NULL)
-    }
-    if (inherits(RunOptions , "RunOptions") == FALSE) {
-      stop("RunOptions must be of class 'RunOptions' \n")
-      return(NULL)
-    }
-    
-    LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
-    
-    if (is.null(Qobs)) {
-      stop("Qobs is missing \n")
-      return(NULL)
-    }
-    if (!is.vector(Qobs)) {
-      stop(paste("Qobs must be a vector of numeric values \n", sep = ""))
-      return(NULL)
-    }
-    if (!is.numeric(Qobs)) {
-      stop(paste("Qobs must be a vector of numeric values \n", sep = ""))
-      return(NULL)
-    }
-    if (length(Qobs) != LLL) {
-      stop("Qobs and InputsModel series must have the same length \n")
-      return(NULL)
-    }
-    if (is.null(BoolCrit)) {
-      BoolCrit <- rep(TRUE, length(Qobs))
-    }
-    if (!is.logical(BoolCrit)) {
-      stop("BoolCrit must be a vector of boolean \n")
-      return(NULL)
-    }
-    if (length(BoolCrit) != LLL) {
-      stop("BoolCrit and InputsModel series must have the same length \n")
-      return(NULL)
-    }
-    if (is.null(transfo)) {
-      stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' or 'sort' \n")
-      return(NULL)
-    }
-    if (!is.vector(transfo)) {
-      stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' or 'sort' \n")
-      return(NULL)
-    }
-    if (length(transfo) != 1) {
-      stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' or 'sort' \n")
-      return(NULL)
-    }
-    if (!is.character(transfo)) {
-      stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' or 'sort' \n")
-      return(NULL)
-    }
-    if (transfo %in% c("", "sqrt", "log", "inv", "sort") == FALSE) {
-      stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' or 'sort' \n")
-      return(NULL)
-    }
-    
-    if (!is.null(Ind_zeroes)) {
-      if (!is.vector(Ind_zeroes)) {
-        stop("Ind_zeroes must be a vector of integers \n")
-        return(NULL)
-      }
-      if (!is.integer(Ind_zeroes)) {
-        stop("Ind_zeroes must be a vector of integers \n")
-        return(NULL)
-      }
-    }
-    if (!is.null(epsilon)) {
-      if (!is.vector(epsilon) |
-          length(epsilon) != 1 | !is.numeric(epsilon)) {
-        stop("epsilon must be single numeric value \n")
-        return(NULL)
-        
-      }
-      epsilon = as.double(epsilon)
-    }
-    
-    ##Create_InputsCrit
-    InputsCrit <- list(BoolCrit = BoolCrit,
-                       Qobs     = Qobs,
-                       transfo  = transfo,
-                       Ind_zeroes = Ind_zeroes,
-                       epsilon = epsilon)
-    
-    class(InputsCrit) <- c("InputsCrit", ObjectClass)
-    
-    return(InputsCrit)
-    
-    
-    
+CreateInputsCrit <- function(FUN_CRIT,
+                             InputsModel,
+                             RunOptions,
+                             Obs,
+                             VarObs = "Q",
+                             BoolCrit = NULL,
+                             transfo = "",
+                             Weights = NULL,
+                             epsilon = NULL,
+                             warnings = TRUE) {
+
+
+  ObjectClass <- NULL
+
+
+  ## ---------- check arguments
+
+  ## check 'InputsModel'
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
   }
+
+
+  ## length of index of period to be used for the model run
+  LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
+
+
+  ## check 'Obs' and definition of idLayer
+  if (!is.numeric(unlist(Obs))) {
+    stop("'Obs' must be a (list of) vector(s) of numeric values")
+  }
+  Obs2 <- Obs
+  if ("ParamT" %in% VarObs) {
+    if (is.list(Obs2)) {
+      Obs2[[which(VarObs == "ParamT")]] <- NULL
+    } else {
+      Obs2 <- NULL
+    }
+  }
+  if (!is.null(Obs2)) {
+    vecObs <- unlist(Obs2)
+    if (length(vecObs) %% LLL != 0) {
+      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", LLL), call. = FALSE)
+    }
+  }
+  if (!is.list(Obs)) {
+    idLayer <- list(1L)
+    Obs <- list(Obs)
+  } else {
+    idLayer <- lapply(Obs, function(i) {
+      if (is.list(i)) {
+        length(i)
+      } else {
+        1L
+      }
+    })
+    Obs <- lapply(Obs, function(x) rowMeans(as.data.frame(x)))
+  }
+
+
+  ## create list of arguments
+  listArgs <- list(FUN_CRIT   = FUN_CRIT,
+                   Obs        = Obs,
+                   VarObs     = VarObs,
+                   BoolCrit   = BoolCrit,
+                   idLayer    = idLayer,
+                   transfo    = as.character(transfo),
+                   Weights    = Weights,
+                   epsilon    = epsilon)
+
+
+  ## check lists lengths
+  for (iArgs in names(listArgs)) {
+    if (iArgs %in% c("Weights", "BoolCrit", "epsilon")) {
+      if (any(is.null(listArgs[[iArgs]]))) {
+        listArgs[[iArgs]] <- lapply(seq_along(listArgs$FUN_CRIT), function(x) NULL)
+      }
+    }
+    if (iArgs %in% c("FUN_CRIT", "VarObs", "transfo", "Weights") & length(listArgs[[iArgs]]) > 1L) {
+      listArgs[[iArgs]] <- as.list(listArgs[[iArgs]])
+    }
+    if (!is.list(listArgs[[iArgs]])) {
+      listArgs[[iArgs]] <- list(listArgs[[iArgs]])
+    }
+  }
+
+  ## check 'FUN_CRIT'
+  listArgs$FUN_CRIT <- lapply(listArgs$FUN_CRIT, FUN = match.fun)
+
+
+  ## check 'VarObs'
+  if (missing(VarObs)) {
+    listArgs$VarObs <- as.list(rep("Q", times = length(listArgs$Obs)))
+    # if (warnings) {
+    #   warning("'VarObs' automatically set to \"Q\"")
+    # }
+  }
+
+
+  ## check 'VarObs' + 'RunOptions'
+  if ("Q" %in% VarObs & !inherits(RunOptions, "GR")) {
+    stop("'VarObs' cannot contain Q if a GR rainfall-runoff model is not used")
+  }
+  if (any(c("SCA", "SWE") %in% VarObs) & !inherits(RunOptions, "CemaNeige")) {
+    stop("'VarObs' cannot contain SCA or SWE if CemaNeige is not used")
+  }
+  if ("SCA" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"Gratio"   %in% RunOptions$Outputs_Sim) {
+    stop("'Gratio' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SCA with CemaNeige")
+  }
+  if ("SWE" %in% VarObs & inherits(RunOptions, "CemaNeige") & !"SnowPack" %in% RunOptions$Outputs_Sim) {
+    stop("'SnowPack' is missing in 'Outputs_Sim' of 'RunOptions', which is necessary to output SWE with CemaNeige")
+  }
+
+
+  ## check 'transfo'
+  if (missing(transfo)) {
+    listArgs$transfo <- as.list(rep("", times = length(listArgs$Obs)))
+    # if (warnings) {
+    #   warning("'transfo' automatically set to \"\"")
+    # }
+  }
+
+  ## check length of each args
+  if (length(unique(sapply(listArgs, FUN = length))) != 1) {
+    stopListArgs <- paste(sapply(names(listArgs), shQuote), collapse = ", ")
+    stop(sprintf("arguments %s must have the same length", stopListArgs))
+  }
+
+
+  ## check 'RunOptions'
+  if (!inherits(RunOptions , "RunOptions")) {
+    stop("'RunOptions' must be of class 'RunOptions'")
+  }
+
+
+  ## check 'Weights'
+  if (length(listArgs$Weights) > 1 & sum(unlist(listArgs$Weights)) == 0 & !any(sapply(listArgs$Weights, is.null))) {
+    stop("sum of 'Weights' cannot be equal to zero")
+  }
+
+
+  ## ---------- reformat
+
+  ## reformat list of arguments
+  listArgs2 <- lapply(seq_along(listArgs$FUN_CRIT), function(i) lapply(listArgs, "[[", i))
+
+  ## preparation of warning messages
+  inVarObs  <- c("Q", "SCA", "SWE", "ParamT")
+  msgVarObs <- "'VarObs' must be a (list of) character vector(s) and one of %s"
+  msgVarObs <- sprintf(msgVarObs, paste(sapply(inVarObs, shQuote), collapse = ", "))
+  inTransfo  <- c("", "sqrt", "log", "inv", "sort", "boxcox") # pow is not checked by inTransfo, but appears in the warning message and checkef after (see ## check 'transfo')
+  msgTransfo <- "'transfo' must be a (list of) character vector(s) and one of %s, or numeric value for power transformation"
+  msgTransfo <- sprintf(msgTransfo, paste(sapply(inTransfo, shQuote), collapse = ", "))
+
+
+  ## ---------- loop on the list of inputs
+
+  InputsCrit <- lapply(listArgs2, function(iListArgs2) {
+
+    ## define FUN_CRIT as a character string
+    iListArgs2$FUN_CRIT <- match.fun(iListArgs2$FUN_CRIT)
+
+    ## check 'FUN_CRIT'
+    if (!all(class(iListArgs2$FUN_CRIT) == c("FUN_CRIT", "function"))) {
+      stop("incorrect 'FUN_CRIT' for use in 'CreateInputsCrit'", call. = FALSE)
+    }
+    if (identical(iListArgs2$FUN_CRIT, ErrorCrit_RMSE) & length(listArgs$Weights) > 1 & all(!is.null(unlist(listArgs$Weights)))) {
+      stop("calculating a composite criterion with the RMSE is not allowed since RMSE is not a dimensionless metric", call. = FALSE)
+    }
+
+    ## check 'Obs'
+    if (iListArgs2$VarObs == "ParamT") {
+      # Parameter for regularisation
+      L2 <- RunOptions$FeatFUN_MOD$NbParam
+    } else {
+      # Observation time series
+      L2 <- LLL
+    }
+    if (!is.vector(iListArgs2$Obs) | length(iListArgs2$Obs) != L2 | !is.numeric(iListArgs2$Obs)) {
+      stop(sprintf("'Obs' must be a (list of) vector(s) of numeric values of length %i", L2), call. = FALSE)
+    }
+
+    ## check 'BoolCrit'
+    if (is.null(iListArgs2$BoolCrit)) {
+      iListArgs2$BoolCrit <- rep(TRUE, length(iListArgs2$Obs))
+    }
+    if (!is.logical(iListArgs2$BoolCrit)) {
+      stop("'BoolCrit' must be a (list of) vector(s) of boolean", call. = FALSE)
+    }
+    if (length(iListArgs2$BoolCrit) != L2) {
+      stop("'BoolCrit' and the period defined in 'RunOptions' must have the same length", call. = FALSE)
+    }
+
+    ## check 'VarObs'
+    if (!is.vector(iListArgs2$VarObs) | length(iListArgs2$VarObs) != 1 | !is.character(iListArgs2$VarObs) | !all(iListArgs2$VarObs %in% inVarObs)) {
+      stop(msgVarObs, call. = FALSE)
+    }
+
+    ## check 'VarObs' + 'Obs'
+    if (any(iListArgs2$VarObs %in% "SCA")) {
+      idSCA <- which(iListArgs2$VarObs == "SCA")
+      if (length(idSCA) == 1L) {
+        vecSCA <- iListArgs2$Obs
+      } else {
+        vecSCA <- unlist(iListArgs2$Obs[idSCA])
+      }
+      if (min(vecSCA, na.rm = TRUE) < 0 | max(vecSCA, na.rm = TRUE) > 1) {
+        stop("'Obs' outside [0,1] for \"SCA\"", call. = FALSE)
+      }
+    }
+    inPosVarObs <- c("Q", "SWE")
+    if (any(iListArgs2$VarObs %in% inPosVarObs)) {
+      idQSS <- which(iListArgs2$VarObs %in% inPosVarObs)
+      if (length(idQSS) == 1L) {
+        vecQSS <- iListArgs2$Obs
+      } else {
+        vecQSS <- unlist(iListArgs2$Obs[idQSS])
+      }
+      if (all(is.na(vecQSS))) {
+        stop("'Obs' contains only missing values", call. = FALSE)
+      }
+      if (min(vecQSS, na.rm = TRUE) < 0) {
+        stop(sprintf("'Obs' outside [0,Inf[ for \"%s\"", iListArgs2$VarObs), call. = FALSE)
+      }
+    }
+
+
+    ## check 'transfo'
+    if (is.null(iListArgs2$transfo) | !is.vector(iListArgs2$transfo) | length(iListArgs2$transfo) != 1 | !is.character(iListArgs2$transfo)) {
+      stop(msgTransfo, call. = FALSE)
+    }
+    isNotInTransfo <- !(iListArgs2$transfo %in% inTransfo)
+    if (any(isNotInTransfo)) {
+      powTransfo <- iListArgs2$transfo[isNotInTransfo]
+      powTransfo <- gsub("\\^|[[:alpha:]]", "", powTransfo)
+      numExpTransfo <- suppressWarnings(as.numeric(powTransfo))
+      if (any(is.na(numExpTransfo))) {
+        stop(msgTransfo, call. = FALSE)
+      }
+      iListArgs2$transfo <- paste0("^", iListArgs2$transfo)
+    }
+
+    ## check 'Weights'
+    if (!is.null(iListArgs2$Weights)) {
+      if (!is.vector(iListArgs2$Weights) | length(iListArgs2$Weights) != 1 | !is.numeric(iListArgs2$Weights) | any(iListArgs2$Weights < 0)) {
+        stop("'Weights' must be a single (list of) positive or equal to zero value(s)", call. = FALSE)
+      }
+    }
+
+    ## check 'epsilon'
+    if (!is.null(iListArgs2$epsilon)) {
+      if (!is.vector(iListArgs2$epsilon) | length(iListArgs2$epsilon) != 1 | !is.numeric(iListArgs2$epsilon) | any(iListArgs2$epsilon <= 0)) {
+        stop("'epsilon' must be a single (list of) positive value(s)", call. = FALSE)
+      }
+    } else if (iListArgs2$transfo %in% c("log", "inv") & any(iListArgs2$Obs %in% 0) & warnings) {
+      warning("zeroes detected in Obs: the corresponding time-steps will be excluded by the 'ErrorCrit*' functions as the epsilon argument was set to NULL", call. = FALSE)
+    }
+
+    ## check 'transfo' + 'FUN_CRIT'
+    if (iListArgs2$transfo == "log" & warnings) {
+      warn_log_kge <- "we do not advise using the %s with a log transformation on Obs (see the details section in the 'CreateInputsCrit' help)"
+      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE)) {
+        warning(sprintf(warn_log_kge, "KGE"), call. = FALSE)
+      }
+      if (identical(iListArgs2$FUN_CRIT, ErrorCrit_KGE2)) {
+        warning(sprintf(warn_log_kge, "KGE'"), call. = FALSE)
+      }
+    }
+
+
+    ## Create InputsCrit
+    iInputsCrit <- list(FUN_CRIT   = iListArgs2$FUN_CRIT,
+                        Obs        = iListArgs2$Obs,
+                        VarObs     = iListArgs2$VarObs,
+                        BoolCrit   = iListArgs2$BoolCrit,
+                        idLayer    = iListArgs2$idLayer,
+                        transfo    = iListArgs2$transfo,
+                        epsilon    = iListArgs2$epsilon,
+                        Weights    = iListArgs2$Weights)
+    class(iInputsCrit) <- c("Single", "InputsCrit", ObjectClass)
+    return(iInputsCrit)
+
+  })
+  names(InputsCrit) <- paste0("IC", seq_along(InputsCrit))
+
+  listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
+  inCnVarObs <- c("SCA", "SWE")
+  if (!"ZLayers" %in% names(InputsModel)) {
+    if (any(listVarObs %in% inCnVarObs)) {
+      stop(sprintf("'VarObs' can not be equal to %i if CemaNeige is not used",
+                   paste(sapply(inCnVarObs, shQuote), collapse = " or ")))
+    }
+  } else {
+    listGroupLayer0 <- sapply(InputsCrit, FUN = "[[", "idLayer")
+    listGroupLayer <- rep(listVarObs, times = listGroupLayer0)
+    tabGroupLayer  <- as.data.frame(table(listGroupLayer))
+    colnames(tabGroupLayer) <- c("VarObs", "freq")
+    nLayers <- length(InputsModel$ZLayers)
+    for (iInCnVarObs in inCnVarObs) {
+      if (any(listVarObs %in% iInCnVarObs)) {
+        if (tabGroupLayer[tabGroupLayer$VarObs %in% iInCnVarObs, "freq"] != nLayers) {
+          stop(sprintf("'Obs' must contain %i vector(s) about %s", nLayers, iInCnVarObs))
+        }
+      }
+    }
+  }
+
+  ## define idLayer as an index of the layer to use
+  for (iInCnVarObs in unique(listVarObs)) {
+    if (!iInCnVarObs %in% inCnVarObs) {
+      for (i in which(listVarObs == iInCnVarObs)) {
+        InputsCrit[[i]]$idLayer <- NA
+      }
+    } else {
+      aa <- listGroupLayer0[listVarObs == iInCnVarObs]
+      aa <- unname(aa)
+      bb <- cumsum(c(0, aa[-length(aa)]))
+      cc <- lapply(seq_along(aa), function(x) seq_len(aa[x]) + bb[x])
+      k <- 1
+      for (i in which(listVarObs == iInCnVarObs)) {
+        InputsCrit[[i]]$idLayer <- cc[[k]]
+        k <- k + 1
+      }
+    }
+  }
+
+
+  ## if only one criterion --> not a list of InputsCrit but directly an InputsCrit
+  if (length(InputsCrit) < 2) {
+    InputsCrit <- InputsCrit[[1L]]
+    InputsCrit["Weights"] <- list(Weights = NULL)
+  } else {
+    if (any(sapply(listArgs$Weights, is.null))) {
+      for (iListArgs in InputsCrit) {
+        iListArgs$Weights <- NULL
+      }
+      class(InputsCrit) <- c("Multi", "InputsCrit", ObjectClass)
+    } else {
+      class(InputsCrit) <- c("Compo", "InputsCrit", ObjectClass)
+    }
+    combInputsCrit <- combn(x = length(InputsCrit), m = 2)
+    apply(combInputsCrit, MARGIN = 2, function(i) {
+      equalInputsCrit <- identical(InputsCrit[[i[1]]], InputsCrit[[i[2]]])
+      if (equalInputsCrit) {
+        warning(sprintf("elements %i and %i of the criteria list are identical. This might not be necessary", i[1], i[2]), call. = FALSE)
+      }
+    })
+  }
+
+  return(InputsCrit)
+
+}

R/CreateInputsCrit_Lavenne.R

+CreateInputsCrit_Lavenne <- function(FUN_CRIT = ErrorCrit_KGE,
+                                      InputsModel,
+                                      RunOptions,
+                                      Obs,
+                                      VarObs = "Q",
+                                      AprParamR,
+                                      AprCrit = 1,
+                                      k = 0.15,
+                                      BoolCrit = NULL,
+                                      transfo = "sqrt",
+                                      epsilon = NULL) {
+
+  # Check parameters
+  if (!is.numeric(AprCrit) || length(AprCrit) != 1 || AprCrit > 1) {
+    stop("'AprCrit' must be a numeric of length 1 with a maximum value of 1")
+  }
+  if (!is.numeric(k) || length(k) != 1 || k < 0 || k > 1) {
+    stop("'k' must be a numeric of length 1 with a value between 0 and 1")
+  }
+  if (!is.null(BoolCrit) && !is.logical(BoolCrit)) {
+    stop("'BoolCrit must be logical")
+  }
+  if (!is.character(transfo)) {
+    stop("'transfo' must be character")
+  }
+  if (!is.null(epsilon) && !is.numeric(epsilon)) {
+    stop("'epsilon' must be numeric")
+  }
+  if (!is.numeric(AprParamR) || length(AprParamR) != RunOptions$FeatFUN_MOD$NbParam) {
+    stop("'AprParamR' must be a numeric vector of length ",
+         RunOptions$FeatFUN_MOD$NbParam)
+  }
+
+
+  FUN_TRANSFO <- .FunTransfo(RunOptions$FeatFUN_MOD)
+
+  AprParamT <- FUN_TRANSFO(AprParamR, "RT")
+
+  ErrorCrit_GAPX <- CreateErrorCrit_GAPX(FUN_TRANSFO)
+
+
+
+  CreateInputsCrit(FUN_CRIT = list(FUN_CRIT, ErrorCrit_GAPX),
+                   InputsModel = InputsModel,
+                   RunOptions = RunOptions,
+                   Obs = list(Obs, AprParamT),
+                   VarObs = c("Q", "ParamT"),
+                   Weights = c(1 - k, k * max(0, AprCrit)),
+                   BoolCrit = list(BoolCrit, NULL),
+                   transfo = list(transfo, ""),
+                   epsilon = list(epsilon, NULL))
+}