Skip to content
GitLab

Source

Target

Commits (1581)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 3402 additions and 1351 deletions
+3402 -1351
.Rbuildignore
View file @ 4bccfee2
^.*\.Rproj$ ^.*\.Rproj$
^\.Rproj\.user$ ^\.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 0 → 100644
View file @ 4bccfee2
{
"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
View file @ 4bccfee2
.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 .Rhistory
.Rapp.history
# Session Data files
.RData .RData
airGR.Rproj
# 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 0 → 100644
View file @ 4bccfee2
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 0 → 100644
View file @ 4bccfee2
# .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 0 → 100644
View file @ 4bccfee2
# 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
View file @ 4bccfee2
Package: airGR Package: airGR
Type: Package Type: Package
Title: Suite of GR hydrological models for precipitation-runoff modelling Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.0.1 Version: 1.7.6.9000
Date: 2016-04-21 Date: 2023-10-25
Authors@R: c( Authors@R: c(
person("Laurent", "Coron", role = c("aut", "cre", "trl")), person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Charles", "Perrin", role = c("aut", "ths")), 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("Claude", "Michel", role = c("aut", "ths")),
person("Vazken", "Andréassian", role = c("ctb", "ths")), person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
person("Pierre", "Brigode", role = c("ctb")), person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260")),
person("Olivier", "Delaigue", role = c("ctb")), person("Pierre", "Brigode", role = c("ctb"), comment = c(ORCID = "0000-0001-8257-0741")),
person("Nicolas", "Le Moine", role = c("ctb")), 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("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("Raji", "Pushpalatha", role = c("ctb")),
person("Guillaume", "Thirel", role = c("ctb")),
person("Audrey", "Valéry", role = c("ctb")) person("Audrey", "Valéry", role = c("ctb"))
) )
Author: Laurent Coron, Charles Perrin, with contributions Depends: R (>= 3.1.0)
from Vazken Andréassian, Pierre Brigode, Olivier Delaigue, Imports:
Nicolas Le Moine, Thibaut Mathevet, Safouane Mouelhi, graphics,
Ludovic Oudin, Raji Pushpalatha, Guillaume Thirel, Audrey Valéry. grDevices,
Based on earlier work by Claude Michel. stats,
Maintainer: Olivier Delaigue <airGR@irstea.fr> utils
Depends: R (>= 3.0.1) Suggests:
Description: This package brings into R the hydrological modelling tools developed knitr, markdown, rmarkdown,
at Irstea-Antony (HBAN Research Unit, France). The package includes several conceptual caRamel, coda, DEoptim, FME, ggmcmc, Rmalschains,
rainfall-runoff models (GR4H, GR4J, GR5J, GR6J, GR2M, GR1A), a snowmelt module (Cemaneige) GGally, ggplot2,
and the associated functions for their calibration and evaluation. Use help(airGR) for package description. 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 License: GPL-2
URL: http://webgr.irstea.fr/modeles/?lang=en URL: https://hydrogr.github.io/airGR/
BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues
NeedsCompilation: yes NeedsCompilation: yes
Encoding: UTF-8 Encoding: UTF-8
VignetteBuilder: knitr
RoxygenNote: 7.1.1
NAMESPACE
View file @ 4bccfee2
##################################### #####################################
## Load DLL ## ## Load DLL ##
##################################### #####################################
useDynLib(airGR) useDynLib(airGR, .registration = TRUE)
#####################################
## S3 methods ##
#####################################
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)
...@@ -11,7 +24,10 @@ useDynLib(airGR) ...@@ -11,7 +24,10 @@ useDynLib(airGR)
export(Calibration) export(Calibration)
export(Calibration_Michel) export(Calibration_Michel)
export(CreateCalibOptions) export(CreateCalibOptions)
export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit) export(CreateInputsCrit)
export(CreateInputsCrit_Lavenne)
export(CreateInputsModel) export(CreateInputsModel)
export(CreateRunOptions) export(CreateRunOptions)
export(DataAltiExtrapolation_Valery) export(DataAltiExtrapolation_Valery)
...@@ -20,26 +36,44 @@ export(ErrorCrit_KGE) ...@@ -20,26 +36,44 @@ export(ErrorCrit_KGE)
export(ErrorCrit_KGE2) export(ErrorCrit_KGE2)
export(ErrorCrit_NSE) export(ErrorCrit_NSE)
export(ErrorCrit_RMSE) export(ErrorCrit_RMSE)
export(PEdaily_Oudin) export(Imax)
export(PE_Oudin)
export(plot.OutputsModel) ### to remove from version 2.0
export(RunModel) export(RunModel)
export(RunModel_CemaNeige) export(RunModel_CemaNeige)
export(RunModel_CemaNeigeGR4H)
export(RunModel_CemaNeigeGR5H)
export(RunModel_CemaNeigeGR4J) export(RunModel_CemaNeigeGR4J)
export(RunModel_CemaNeigeGR5J) export(RunModel_CemaNeigeGR5J)
export(RunModel_CemaNeigeGR6J) export(RunModel_CemaNeigeGR6J)
export(RunModel_GR1A) export(RunModel_GR1A)
export(RunModel_GR2M) export(RunModel_GR2M)
export(RunModel_GR4H) export(RunModel_GR4H)
export(RunModel_GR5H)
export(RunModel_GR4J) export(RunModel_GR4J)
export(RunModel_GR5J) export(RunModel_GR5J)
export(RunModel_GR6J) export(RunModel_GR6J)
export(RunModel_Lag)
export(SeriesAggreg) export(SeriesAggreg)
export(TransfoParam) export(TransfoParam)
export(TransfoParam_CemaNeige) export(TransfoParam_CemaNeige)
export(TransfoParam_CemaNeigeHyst)
export(TransfoParam_GR1A) export(TransfoParam_GR1A)
export(TransfoParam_GR2M) export(TransfoParam_GR2M)
export(TransfoParam_GR4H) export(TransfoParam_GR4H)
export(TransfoParam_GR5H)
export(TransfoParam_GR4J) export(TransfoParam_GR4J)
export(TransfoParam_GR5J) export(TransfoParam_GR5J)
export(TransfoParam_GR6J) export(TransfoParam_GR6J)
export(plot_OutputsModel) export(TransfoParam_Lag)
#export(.ErrorCrit)
#export(.FeatModels)
#####################################
## Import ##
#####################################
import(stats)
import(graphics)
import(grDevices)
import(utils)
NEWS.md 0 → 100644
View file @ 4bccfee2
This diff is collapsed.
R/BasinData.R deleted 100644 → 0
View file @ 749e9fa9
#' @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
View file @ 4bccfee2
Calibration <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=NULL, verbose = TRUE){ Calibration <- function(InputsModel,
return( FUN_CALIB(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO, verbose = verbose) ) 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/Calibration_Michel.R
View file @ 4bccfee2
This diff is collapsed.
R/CreateCalibOptions.R
View file @ 4bccfee2
CreateCalibOptions <- function(FUN_MOD,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=NULL,OptimParam=NULL,FixedParam=NULL,SearchRanges=NULL, CreateCalibOptions <- function(FUN_MOD,
StartParam=NULL,StartParamList=NULL,StartParamDistrib=NULL){ FUN_CALIB = Calibration_Michel,
FUN_TRANSFO = NULL,
IsHyst = FALSE,
ObjectClass <- NULL; IsSD = FALSE,
FixedParam = NULL,
##check_FUN_MOD SearchRanges = NULL,
BOOL <- FALSE; StartParamList = NULL,
if(identical(FUN_MOD,RunModel_GR4H )){ ObjectClass <- c(ObjectClass,"GR4H" ); BOOL <- TRUE; } StartParamDistrib = NULL) {
if(identical(FUN_MOD,RunModel_GR4J )){ ObjectClass <- c(ObjectClass,"GR4J" ); BOOL <- TRUE; }
if(identical(FUN_MOD,RunModel_GR5J )){ ObjectClass <- c(ObjectClass,"GR5J" ); BOOL <- TRUE; } ObjectClass <- NULL
if(identical(FUN_MOD,RunModel_GR6J )){ ObjectClass <- c(ObjectClass,"GR6J" ); BOOL <- TRUE; }
if(identical(FUN_MOD,RunModel_GR2M )){ ObjectClass <- c(ObjectClass,"GR2M" ); BOOL <- TRUE; } FUN_MOD <- match.fun(FUN_MOD)
if(identical(FUN_MOD,RunModel_GR1A )){ ObjectClass <- c(ObjectClass,"GR1A" ); BOOL <- TRUE; } FUN_CALIB <- match.fun(FUN_CALIB)
if(identical(FUN_MOD,RunModel_CemaNeige )){ ObjectClass <- c(ObjectClass,"CemaNeige" ); BOOL <- TRUE; } if (!is.null(FUN_TRANSFO)) {
if(identical(FUN_MOD,RunModel_CemaNeigeGR4J)){ ObjectClass <- c(ObjectClass,"CemaNeigeGR4J"); BOOL <- TRUE; } FUN_TRANSFO <- match.fun(FUN_TRANSFO)
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 (!is.logical(IsHyst) | length(IsHyst) != 1L) {
if(!BOOL){ stop("incorrect FUN_MOD for use in CreateCalibOptions \n"); return(NULL); } stop("'IsHyst' must be a logical of length 1")
}
##check_FUN_CALIB if (!is.logical(IsSD) | length(IsSD) != 1L) {
BOOL <- FALSE; stop("'IsSD' must be a logical of length 1")
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_MOD
##check_FUN_TRANSFO FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD)
if(is.null(FUN_TRANSFO)){ FeatFUN_MOD$IsHyst <- IsHyst
##_set_FUN1 FeatFUN_MOD$IsSD <- IsSD
if(identical(FUN_MOD,RunModel_GR4H ) ){ FUN1 <- TransfoParam_GR4H ; } ObjectClass <- FeatFUN_MOD$Class
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_Lag) && IsSD) {
if(identical(FUN_MOD,RunModel_GR6J ) | identical(FUN_MOD,RunModel_CemaNeigeGR6J) ){ FUN1 <- TransfoParam_GR6J ; } stop("RunModel_Lag should not be used with 'IsSD=TRUE'")
if(identical(FUN_MOD,RunModel_GR2M ) ){ FUN1 <- TransfoParam_GR2M ; } }
if(identical(FUN_MOD,RunModel_GR1A ) ){ FUN1 <- TransfoParam_GR1A ; } if (IsHyst) {
if(identical(FUN_MOD,RunModel_CemaNeige) ){ FUN1 <- TransfoParam_CemaNeige; } ObjectClass <- c(ObjectClass, "hysteresis")
if(is.null(FUN1)){ stop("FUN1 was not found \n"); return(NULL); } }
##_set_FUN2 if (IsSD) {
FUN2 <- TransfoParam_CemaNeige; ObjectClass <- c(ObjectClass, "SD")
##_set_FUN_TRANSFO }
if(sum(ObjectClass %in% c("GR4H","GR4J","GR5J","GR6J","GR2M","GR1A","CemaNeige"))>0){
FUN_TRANSFO <- FUN1; ## check FUN_CALIB
} else { BOOL <- FALSE
FUN_TRANSFO <- function(ParamIn,Direction){
Bool <- is.matrix(ParamIn); if (identical(FUN_CALIB, Calibration_Michel)) {
if(Bool==FALSE){ ParamIn <- rbind(ParamIn); } ObjectClass <- c(ObjectClass, "HBAN")
ParamOut <- NA*ParamIn; BOOL <- TRUE
NParam <- ncol(ParamIn); }
if(NParam <= 3){ if (!BOOL) {
ParamOut[, 1:(NParam-2)] <- FUN1(cbind(ParamIn[,1:(NParam-2)]),Direction); stop("incorrect 'FUN_CALIB' for use in 'CreateCalibOptions'")
} else { return(NULL)
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); ## check FUN_TRANSFO
} if (is.null(FUN_TRANSFO)) {
} FUN_TRANSFO <- .FunTransfo(FeatFUN_MOD)
} }
if(is.null(FUN_TRANSFO)){ stop("FUN_TRANSFO was not found \n"); return(NULL); }
## NParam
##NParam NParam <- FeatFUN_MOD$NbParam
if("GR4H" %in% ObjectClass){ NParam <- 4; }
if("GR4J" %in% ObjectClass){ NParam <- 4; } if (IsHyst) {
if("GR5J" %in% ObjectClass){ NParam <- 5; } NParam <- NParam + 2
if("GR6J" %in% ObjectClass){ NParam <- 6; } }
if("GR2M" %in% ObjectClass){ NParam <- 2; } if (IsSD) {
if("GR1A" %in% ObjectClass){ NParam <- 1; } NParam <- NParam + 1
if("CemaNeige" %in% ObjectClass){ NParam <- 2; } }
if("CemaNeigeGR4J" %in% ObjectClass){ NParam <- 6; }
if("CemaNeigeGR5J" %in% ObjectClass){ NParam <- 7; } ## check FixedParam
if("CemaNeigeGR6J" %in% ObjectClass){ NParam <- 8; } if (is.null(FixedParam)) {
FixedParam <- rep(NA, NParam)
##check_OptimParam } else {
if(is.null(OptimParam)){ if (!is.vector(FixedParam)) {
OptimParam <- rep(TRUE,NParam); stop("FixedParam must be a vector")
} else {
if(!is.vector(OptimParam) ){ stop("OptimParam must be a vector of booleans \n"); return(NULL); }
if(length(OptimParam)!=NParam){ stop("Incompatibility between OptimParam length and FUN_MOD \n"); return(NULL); }
if(!is.logical(OptimParam) ){ stop("OptimParam must be a vector of booleans \n"); return(NULL); }
}
##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 OptimParam length and FUN_MOD \n"); return(NULL); }
if(sum(!OptimParam)>0){
if(!is.numeric(FixedParam[!OptimParam])){ stop("if OptimParam[i]==FALSE, FixedParam[i] must be a numeric value \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( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81,
+5.02, +0.61, +3.45, -8.68, +1.95, +2.27,
+5.58, +0.78, +4.18, -8.12, +3.59, +3.56),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( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81, -9.96, +6.63,
+5.02, +0.61, +3.45, -8.68, +1.95, +2.27, -9.14, +6.90,
+5.58, +0.78, +4.18, -8.12, +3.59, +3.56, +4.10, +7.21),ncol=NParam,byrow=TRUE); }
StartParamList <- NULL;
StartParamDistrib <- TransfoParam(ParamIn=ParamT,Direction="TR",FUN_TRANSFO=FUN_TRANSFO);
StartParam <- StartParamDistrib[2,];
}
##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); }
} }
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)
##Create_CalibOptions } else {
CalibOptions <- list(OptimParam=OptimParam,FixedParam=FixedParam,SearchRanges=SearchRanges); if (!is.matrix(SearchRanges)) {
if(!is.null(StartParam )){ CalibOptions <- c(CalibOptions,list(StartParam=StartParam)); } stop("'SearchRanges' must be a matrix")
if(!is.null(StartParamList )){ CalibOptions <- c(CalibOptions,list(StartParamList=StartParamList)); } }
if(!is.null(StartParamDistrib)){ CalibOptions <- c(CalibOptions,list(StartParamDistrib=StartParamDistrib)); } if (!is.numeric(SearchRanges)) {
class(CalibOptions) <- c("CalibOptions",ObjectClass); stop("'SearchRanges' must be a matrix of numeric values")
return(CalibOptions); }
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 0 → 100644
View file @ 4bccfee2
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 0 → 100644
View file @ 4bccfee2
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
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 (!(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'")
}
if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
stop("'InputsModel' must be of class 'GR'")
}
if ("CemaNeige" %in% ObjectClass &
!inherits(InputsModel, "CemaNeige")) {
stop("'InputsModel' must be of class 'CemaNeige'")
}
## 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")
}
} else if (!is.null(ExpStore)) {
if (verbose) {
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", FeatFUN_MOD$NameFunMod))
}
UH1 <- rep(Inf, UH1n)
}
if (!is.null(UH2)) {
if (verbose) {
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", 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", FeatFUN_MOD$NameFunMod))
}
}
ProdStore <- Inf
if (!is.null(RoutStore)) {
if (verbose) {
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", 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", FeatFUN_MOD$NameFunMod))
}
}
UH1 <- rep(Inf, UH1n)
if (!is.null(UH2)) {
if (verbose) {
warning(sprintf("'%s' does not require 'UH2'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
}
UH2 <- rep(Inf, UH2n)
}
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(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) | !is.null(GthrCemaNeigeLayers) | !is.null(GlocmaxCemaNeigeLayers))) {
if (verbose) {
warning(sprintf("'%s' does not require 'GCemaNeigeLayers' 'GCemaNeigeLayers', 'GthrCemaNeigeLayers' and 'GlocmaxCemaNeigeLayers'. Values set to NA", FeatFUN_MOD$NameFunMod))
}
GCemaNeigeLayers <- Inf
eTGCemaNeigeLayers <- Inf
GthrCemaNeigeLayers <- Inf
GlocmaxCemaNeigeLayers <- Inf
}
## set states
if ("CemaNeige" %in% ObjectClass) {
NLayers <- length(InputsModel$LayerPrecip)
} else {
NLayers <- 1
}
## manage NULL values
if (is.null(ExpStore)) {
ExpStore <- Inf
}
if (is.null(IntStore)) {
IntStore <- Inf
}
if (is.null(UH1)) {
if ("hourly" %in% ObjectClass) {
k <- 24
} else {
k <- 1
}
UH1 <- rep(Inf, UH1n * k)
}
if (is.null(UH2)) {
if ("hourly" %in% ObjectClass) {
k <- 24
} else {
k <- 1
}
UH2 <- rep(Inf, UH2n * k)
}
if (is.null(GCemaNeigeLayers)) {
GCemaNeigeLayers <- rep(Inf, NLayers)
}
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) | any(IntStore < 0) |
any(UH1 < 0) | any(UH2 < 0) |
any(GCemaNeigeLayers < 0)) {
stop("Negative values are not allowed for any of 'ProdStore', 'RoutStore', 'IntStore', 'UH1', 'UH2', 'GCemaNeigeLayers'")
}
## check length
if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
stop("'ProdStore' must be numeric of length one")
}
if (!is.numeric(RoutStore) || length(RoutStore) != 1L) {
stop("'RoutStore' must be numeric of length one")
}
if (!is.numeric(ExpStore) || length(ExpStore) != 1L) {
stop("'ExpStore' must be numeric of length one")
}
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) != UH1n)) {
stop(sprintf("'UH1' must be numeric of length %i", UH1n))
}
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) != UH2n)) {
stop(sprintf("'UH2' must be numeric of length %i (2 * %i)", UH2n, UH1n))
}
if (!is.numeric(GCemaNeigeLayers) || length(GCemaNeigeLayers) != NLayers) {
stop(sprintf("'GCemaNeigeLayers' must be numeric of length %i", NLayers))
}
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, Int = IntStore),
UH = list(UH1 = UH1, UH2 = UH2),
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
View file @ 4bccfee2
CreateInputsCrit <- function(FUN_CRIT,InputsModel,RunOptions,Qobs,BoolCrit=NULL,transfo="",Ind_zeroes=NULL,epsilon=NULL){ CreateInputsCrit <- function(FUN_CRIT,
InputsModel,
RunOptions,
Obs,
VarObs = "Q",
BoolCrit = NULL,
transfo = "",
Weights = NULL,
epsilon = NULL,
warnings = TRUE) {
ObjectClass <- NULL;
##check_FUN_CRIT ObjectClass <- NULL
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); } ## ---------- check arguments
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)); } ## check 'InputsModel'
if(!is.logical(BoolCrit)){ stop("BoolCrit must be a vector of boolean \n" ); return(NULL); } if (!inherits(InputsModel, "InputsModel")) {
if(length(BoolCrit)!=LLL){ stop("BoolCrit and InputsModel series must have the same length \n"); return(NULL); } stop("'InputsModel' must be of class 'InputsModel'")
}
if(is.null(transfo) ){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); }
if(!is.vector(transfo )){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); }
if(length(transfo)!=1 ){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); }
if(!is.character(transfo)){ stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); }
if(transfo %in% c("","sqrt","log","inv") == FALSE){
stop("transfo must be a chosen among the following: '', 'sqrt', 'log' or 'inv' \n"); return(NULL); }
if(!is.null(Ind_zeroes)){ ## length of index of period to be used for the model run
if(!is.vector( Ind_zeroes)){ stop("Ind_zeroes must be a vector of integers \n" ); return(NULL); } LLL <- length(InputsModel$DatesR[RunOptions$IndPeriod_Run])
if(!is.integer(Ind_zeroes)){ stop("Ind_zeroes must be a vector of integers \n" ); return(NULL); }
## 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(epsilon)){ }
if(!is.vector( epsilon) | length(epsilon)!=1 | !is.numeric(epsilon)){ if (!is.null(Obs2)) {
stop("epsilon must be single numeric value \n" ); return(NULL); } vecObs <- unlist(Obs2)
epsilon=as.double(epsilon); 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_InputsCrit
InputsCrit <- list(BoolCrit=BoolCrit,Qobs=Qobs,transfo=transfo,Ind_zeroes=Ind_zeroes,epsilon=epsilon);
class(InputsCrit) <- c("InputsCrit",ObjectClass);
return(InputsCrit);
## 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 0 → 100644
View file @ 4bccfee2
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))
}
R/CreateInputsModel.R
View file @ 4bccfee2
CreateInputsModel <- function(FUN_MOD,DatesR,Precip,PotEvap=NULL,TempMean=NULL,TempMin=NULL,TempMax=NULL,ZInputs=NULL,HypsoData=NULL,NLayers=5, verbose = TRUE){ CreateInputsModel <- function(FUN_MOD,
DatesR,
Precip, PrecipScale = TRUE,
PotEvap = NULL,
TempMean = NULL, TempMin = NULL, TempMax = NULL,
ZInputs = NULL, HypsoData = NULL, NLayers = 5,
Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
QupstrUnit = "mm",
verbose = TRUE) {
ObjectClass <- NULL;
##check_FUN_MOD ObjectClass <- NULL
BOOL <- FALSE;
if(identical(FUN_MOD,RunModel_GR4H)){ ## check DatesR
ObjectClass <- c(ObjectClass,"hourly","GR"); if (is.null(DatesR)) {
TimeStep <- as.integer(60*60); stop("'DatesR' is missing")
BOOL <- TRUE; }
if (!"POSIXlt" %in% class(DatesR) & !"POSIXct" %in% class(DatesR)) {
stop("'DatesR' must be defined as 'POSIXlt' or 'POSIXct'")
}
if (!"POSIXlt" %in% class(DatesR)) {
DatesR <- as.POSIXlt(DatesR)
}
if (any(duplicated(DatesR))) {
stop("'DatesR' must not include duplicated values")
}
LLL <- length(DatesR)
## check FUN_MOD
FUN_MOD <- match.fun(FUN_MOD)
FeatFUN_MOD <- .GetFeatModel(FUN_MOD = FUN_MOD, DatesR = DatesR)
ObjectClass <- FeatFUN_MOD$Class
TimeStep <- FeatFUN_MOD$TimeStep
##check_arguments
if ("GR" %in% ObjectClass) {
if (is.null(Precip)) {
stop("Precip is missing")
}
if (is.null(PotEvap)) {
stop("'PotEvap' is missing")
}
if (!is.vector(Precip) | !is.vector(PotEvap)) {
stop("'Precip' and 'PotEvap' must be vectors of numeric values")
} }
if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_GR6J)){ if (!is.numeric(Precip) | !is.numeric(PotEvap)) {
ObjectClass <- c(ObjectClass,"daily","GR"); stop("'Precip' and 'PotEvap' must be vectors of numeric values")
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
} }
if(identical(FUN_MOD,RunModel_GR2M)){ if (length(Precip) != LLL | length(PotEvap) != LLL) {
ObjectClass <- c(ObjectClass,"GR","monthly"); stop("'Precip', 'PotEvap' and 'DatesR' must have the same length")
TimeStep <- as.integer(c(28,29,30,31)*24*60*60);
BOOL <- TRUE;
} }
if(identical(FUN_MOD,RunModel_GR1A)){ }
ObjectClass <- c(ObjectClass,"GR","yearly"); if ("CemaNeige" %in% ObjectClass) {
TimeStep <- as.integer(c(365,366)*24*60*60); if (is.null(Precip)) {
BOOL <- TRUE; stop("'Precip' is missing")
} }
if(identical(FUN_MOD,RunModel_CemaNeige)){ if (is.null(TempMean)) {
ObjectClass <- c(ObjectClass,"daily","CemaNeige"); stop("'TempMean' is missing")
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
} }
if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){ if (!is.vector(Precip) | !is.vector(TempMean)) {
ObjectClass <- c(ObjectClass,"daily","GR","CemaNeige"); stop("'Precip' and 'TempMean' must be vectors of numeric values")
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
} }
if(!BOOL){ stop("incorrect FUN_MOD for use in CreateInputsModel \n"); return(NULL); } if (!is.numeric(Precip) | !is.numeric(TempMean)) {
stop("'Precip' and 'TempMean' must be vectors of numeric values")
}
if (length(Precip) != LLL | length(TempMean) != LLL) {
stop("'Precip', 'TempMean' and 'DatesR' must have the same length")
}
if (is.null(TempMin) != is.null(TempMax)) {
stop("'TempMin' and 'TempMax' must be both defined if not null")
}
if (!is.null(TempMin) & !is.null(TempMax)) {
if (!is.vector(TempMin) | !is.vector(TempMax)) {
stop("'TempMin' and 'TempMax' must be vectors of numeric values")
}
if (!is.numeric(TempMin) | !is.numeric(TempMax)) {
stop("'TempMin' and 'TempMax' must be vectors of numeric values")
}
if (length(TempMin) != LLL | length(TempMax) != LLL) {
stop("'TempMin', 'TempMax' and 'DatesR' must have the same length")
}
}
if (!is.null(HypsoData)) {
if (!is.vector(HypsoData)) {
stop("'HypsoData' must be a vector of numeric values if not null")
}
if (!is.numeric(HypsoData)) {
stop("'HypsoData' must be a vector of numeric values if not null")
}
if (length(HypsoData) != 101) {
stop("'HypsoData' must be of length 101 if not null")
}
if (sum(is.na(HypsoData)) != 0 & sum(is.na(HypsoData)) != 101) {
stop("'HypsoData' must not contain any NA if not null")
}
}
if (!is.null(ZInputs)) {
if (length(ZInputs) != 1) {
stop("'ZInputs' must be a single numeric value if not null")
}
if (is.na(ZInputs) | !is.numeric(ZInputs)) {
stop("'ZInputs' must be a single numeric value if not null")
}
}
if (is.null(HypsoData)) {
if (verbose) {
warning("'HypsoData' is missing: a single layer is used and no extrapolation is made", call. = FALSE)
}
HypsoData <- as.numeric(rep(NA, 101))
ZInputs <- as.numeric(NA)
NLayers <- as.integer(1)
##check_arguments }
if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){ if (is.null(ZInputs)) {
if(is.null(DatesR)){ stop("DatesR is missing \n"); return(NULL); } if (verbose & !identical(HypsoData, as.numeric(rep(NA, 101)))) {
if("POSIXlt" %in% class(DatesR) == FALSE & "POSIXct" %in% class(DatesR) == FALSE){ stop("DatesR must be defined as POSIXlt or POSIXct \n"); return(NULL); } warning("'ZInputs' is missing: HypsoData[51] is used")
if("POSIXlt" %in% class(DatesR) == FALSE){ DatesR <- as.POSIXlt(DatesR); }
if(difftime(tail(DatesR,1),tail(DatesR,2),units="secs")[[1]] %in% TimeStep==FALSE){ stop(paste("the time step of the model inputs must be ",TimeStep," seconds \n",sep="")); return(NULL); }
LLL <- length(DatesR);
}
if("GR" %in% ObjectClass){
if(is.null(Precip )){ stop("Precip is missing \n" ); return(NULL); }
if(is.null(PotEvap )){ stop("PotEvap is missing \n" ); return(NULL); }
if(!is.vector( Precip) | !is.vector( PotEvap)){ stop("Precip and PotEvap must be vectors of numeric values \n"); return(NULL); }
if(!is.numeric(Precip) | !is.numeric(PotEvap)){ stop("Precip and PotEvap must be vectors of numeric values \n"); return(NULL); }
if(length(Precip)!=LLL | length(PotEvap)!=LLL){ stop("Precip, PotEvap and DatesR must have the same length \n"); return(NULL); }
}
if("CemaNeige" %in% ObjectClass){
if(is.null(Precip )){ stop("Precip is missing \n" ); return(NULL); }
if(is.null(TempMean)){ stop("TempMean is missing \n"); return(NULL); }
if(!is.vector( Precip) | !is.vector( TempMean)){ stop("Precip and TempMean must be vectors of numeric values \n"); return(NULL); }
if(!is.numeric(Precip) | !is.numeric(TempMean)){ stop("Precip and TempMean must be vectors of numeric values \n"); return(NULL); }
if(length(Precip)!=LLL | length(TempMean)!=LLL){ stop("Precip, TempMean and DatesR must have the same length \n"); return(NULL); }
if(is.null(TempMin)!=is.null(TempMax)){ stop("TempMin and TempMax must be both defined if not null \n"); return(NULL); }
if(!is.null(TempMin) & !is.null(TempMax)){
if(!is.vector( TempMin) | !is.vector( TempMax)){ stop("TempMin and TempMax must be vectors of numeric values \n"); return(NULL); }
if(!is.numeric(TempMin) | !is.numeric(TempMax)){ stop("TempMin and TempMax must be vectors of numeric values \n"); return(NULL); }
if(length(TempMin)!=LLL | length(TempMax)!=LLL){ stop("TempMin, TempMax and DatesR must have the same length \n"); return(NULL); }
}
if(!is.null(HypsoData)){
if(!is.vector( HypsoData)){ stop("HypsoData must be a vector of numeric values if not null \n"); return(NULL); }
if(!is.numeric(HypsoData)){ stop("HypsoData must be a vector of numeric values if not null \n"); return(NULL); }
if(length(HypsoData)!=101){ stop("HypsoData must be of length 101 if not null \n"); return(NULL); }
if(sum(is.na(HypsoData))!=0 & sum(is.na(HypsoData))!=101){ stop("HypsoData must not contain any NA if not null \n"); return(NULL); }
}
if(!is.null(ZInputs)){
if(length(ZInputs)!=1 ){ stop("\t ZInputs must be a single numeric value if not null \n"); return(NULL); }
if(is.na(ZInputs) | !is.numeric(ZInputs)){ stop("\t ZInputs must be a single numeric value if not null \n"); return(NULL); }
}
if(is.null(HypsoData)){
if(verbose){ warning("\t HypsoData is missing => a single layer is used and no extrapolation is made \n"); }
HypsoData <- as.numeric(rep(NA,101)); ZInputs <- as.numeric(NA); NLayers <- as.integer(1);
}
if(is.null(ZInputs)){
if(verbose & !identical(HypsoData,as.numeric(rep(NA,101)))){ warning("\t ZInputs is missing => HypsoData[51] is used \n"); }
ZInputs <- HypsoData[51];
} }
ZInputs <- HypsoData[51L]
}
if (NLayers <= 0) {
stop("'NLayers' must be a positive integer value")
}
if (NLayers != as.integer(NLayers)) {
warning("Coerce 'NLayers' to be of integer type (", NLayers, ": ", as.integer(NLayers), ")")
NLayers <- as.integer(NLayers)
} }
}
## check semi-distributed mode
if (!is.null(Qupstream) & !is.null(LengthHydro) & !is.null(BasinAreas)) {
ObjectClass <- c(ObjectClass, "SD")
} else if (verbose & !all(c(is.null(Qupstream), is.null(LengthHydro), is.null(BasinAreas)))) {
warning("Missing argument: 'Qupstream', 'LengthHydro' and 'BasinAreas' must all be set to run in a semi-distributed mode. The lumped mode will be used")
}
if ("SD" %in% ObjectClass) {
if (!("daily" %in% ObjectClass) & !("hourly" %in% ObjectClass)) {
stop("Only daily and hourly time steps can be used in a semi-distributed mode")
}
if (!is.matrix(Qupstream) | !is.numeric(Qupstream)) {
stop("'Qupstream' must be a matrice of numeric values")
}
if (!is.vector(LengthHydro) | !is.vector(BasinAreas) | !is.numeric(LengthHydro) | !is.numeric(BasinAreas)) {
stop("'LengthHydro' and 'BasinAreas' must be vectors of numeric values")
}
if (ncol(Qupstream) != length(LengthHydro)) {
stop("'Qupstream' number of columns and 'LengthHydro' length must be equal")
}
if (length(LengthHydro) + 1 != length(BasinAreas)) {
stop("'BasinAreas' must have one more element than 'LengthHydro'")
}
if (nrow(Qupstream) != LLL) {
stop("'Qupstream' must have same number of rows as 'DatesR' length")
}
if (any(is.na(Qupstream))) {
warning("'Qupstream' contains NA values: model outputs will contain NAs")
}
if (any(LengthHydro > 1000)) {
warning("The unit of 'LengthHydro' has changed from m to km in airGR >= 1.6.12: values superior to 1000 km seem unrealistic")
}
QupstrUnit <- tolower(QupstrUnit)
QupstrUnit <- match.arg(arg = QupstrUnit, choices = c("mm", "m3", "m3/s", "l/s"))
}
##check_NA_values ##check_NA_values
BOOL_NA <- rep(FALSE,length(DatesR)); BOOL_NA <- rep(FALSE, length(DatesR))
if("GR" %in% ObjectClass){
BOOL_NA_TMP <- (Precip < 0) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in Precip series \n"); } } if ("GR" %in% ObjectClass) {
BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in PotEvap series \n"); } } BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
} if (sum(BOOL_NA_TMP) != 0) {
if("CemaNeige" %in% ObjectClass){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
BOOL_NA_TMP <- (Precip < 0 ) | is.na(Precip ); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < 0 or NA values detected in Precip series \n"); } } if (verbose) {
BOOL_NA_TMP <- (TempMean<(-150)) | is.na(TempMean); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMean series \n"); } } warning("Values < 0 or NA values detected in 'Precip' series")
if(!is.null(TempMin) & !is.null(TempMax)){ }
BOOL_NA_TMP <- (TempMin<(-150)) | is.na(TempMin); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMin series \n"); } } }
BOOL_NA_TMP <- (TempMax<(-150)) | is.na(TempMax); if(sum(BOOL_NA_TMP)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP; if(verbose){ warning("\t Values < -150) or NA values detected in TempMax series \n"); } } } BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap)
} if (sum(BOOL_NA_TMP) != 0) {
if(sum(BOOL_NA)!=0){ BOOL_NA <- BOOL_NA | BOOL_NA_TMP
WTxt <- NULL; if (verbose) {
WTxt <- paste(WTxt,"\t Missing values are not allowed in InputsModel \n",sep=""); warning("Values < 0 or NA values detected in 'PotEvap' series")
Select <- (max(which(BOOL_NA))+1):length(BOOL_NA); }
if(Select[1]>Select[2]){ stop(paste("time series could not be trunced since missing values were detected at the list time-step \n",sep="")); return(NULL); } }
if("GR" %in% ObjectClass){ }
Precip <- Precip[Select]; PotEvap <- PotEvap[Select]; } if ("CemaNeige" %in% ObjectClass) {
if("CemaNeige" %in% ObjectClass){ BOOL_NA_TMP <- (Precip < 0) | is.na(Precip)
Precip <- Precip[Select]; TempMean <- TempMean[Select]; if(!is.null(TempMin) & !is.null(TempMax)){ TempMin <- TempMin[Select]; TempMax <- TempMax[Select]; } } if (sum(BOOL_NA_TMP) != 0) {
DatesR <- DatesR[Select]; BOOL_NA <- BOOL_NA | BOOL_NA_TMP
WTxt <- paste(WTxt,"\t -> data were trunced to keep the most recent available time-steps \n",sep=""); if (verbose) {
WTxt <- paste(WTxt,"\t -> ",length(Select)," time-steps were kept \n",sep=""); warning("Values < 0 or NA values detected in 'Precip' series")
if(!is.null(WTxt) & verbose){ warning(WTxt); } }
}
BOOL_NA_TMP <- (TempMean < (-150)) | is.na(TempMean)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
warning("Values < -150 or NA values detected in 'TempMean' series")
}
}
if (!is.null(TempMin) & !is.null(TempMax)) {
BOOL_NA_TMP <- (TempMin < (-150)) | is.na(TempMin)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
warning("Values < -150 or NA values detected in 'TempMin' series")
}
}
BOOL_NA_TMP <- (TempMax < (-150)) | is.na(TempMax)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
warning("Values < -150 or NA values detected in 'TempMax' series")
}
}
} }
}
if (sum(BOOL_NA) != 0) {
WTxt <- NULL
WTxt <- paste(WTxt, "\t Missing values are not allowed in 'InputsModel'", sep = "")
Select <- (max(which(BOOL_NA)) + 1):length(BOOL_NA)
if (Select[1L] > Select[2L]) {
stop("time series could not be trunced since missing values were detected at the last time-step")
}
if ("GR" %in% ObjectClass) {
Precip <- Precip[Select]
PotEvap <- PotEvap[Select]
}
if ("CemaNeige" %in% ObjectClass) {
Precip <- Precip[Select]
TempMean <- TempMean[Select]
if (!is.null(TempMin) & !is.null(TempMax)) {
TempMin <- TempMin[Select]
TempMax <- TempMax[Select]
}
}
DatesR <- DatesR[Select]
WTxt <- paste0(WTxt, "\t -> data were trunced to keep the most recent available time-steps")
WTxt <- paste0(WTxt, "\t -> ", length(Select), " time-steps were kept")
if (!is.null(WTxt) & verbose) {
warning(WTxt)
}
}
##DataAltiExtrapolation_Valery ##DataAltiExtrapolation_Valery
if("CemaNeige" %in% ObjectClass){ if ("CemaNeige" %in% ObjectClass) {
RESULT <- DataAltiExtrapolation_Valery(DatesR=DatesR,Precip=Precip,TempMean=TempMean,TempMin=TempMin,TempMax=TempMax,ZInputs=ZInputs,HypsoData=HypsoData,NLayers=NLayers, verbose = verbose); RESULT <- DataAltiExtrapolation_Valery(DatesR = DatesR,
if(verbose){ if(NLayers==1){ cat(paste("\t Input series were successfully created on 1 elevation layer for use by CemaNeige \n",sep="")); Precip = Precip, PrecipScale = PrecipScale,
} else { cat(paste("\t Input series were successfully created on ",NLayers," elevation layers for use by CemaNeige \n",sep="")); } } TempMean = TempMean, TempMin = TempMin, TempMax = TempMax,
ZInputs = ZInputs, HypsoData = HypsoData, NLayers = NLayers,
verbose = verbose)
if (verbose) {
if (NLayers == 1) {
message("input series were successfully created on 1 elevation layer for use by CemaNeige")
} else {
message( "input series were successfully created on ", NLayers, " elevation layers for use by CemaNeige")
}
} }
}
##Create_InputsModel ##Create_InputsModel
InputsModel <- list(DatesR=DatesR); InputsModel <- list(DatesR = DatesR)
if("GR" %in% ObjectClass){ if ("GR" %in% ObjectClass) {
InputsModel <- c(InputsModel,list(Precip=as.double(Precip),PotEvap=as.double(PotEvap))); } InputsModel <- c(InputsModel, list(Precip = as.double(Precip), PotEvap = as.double(PotEvap)))
if("CemaNeige" %in% ObjectClass){ }
InputsModel <- c(InputsModel,list(LayerPrecip=RESULT$LayerPrecip,LayerTempMean=RESULT$LayerTempMean, if ("CemaNeige" %in% ObjectClass) {
LayerFracSolidPrecip=RESULT$LayerFracSolidPrecip,ZLayers=RESULT$ZLayers)); } InputsModel <- c(InputsModel, list(LayerPrecip = RESULT$LayerPrecip,
LayerTempMean = RESULT$LayerTempMean,
LayerFracSolidPrecip = RESULT$LayerFracSolidPrecip,
ZLayers = RESULT$ZLayers))
}
if ("SD" %in% ObjectClass) {
# Qupstream is internally stored in m3/time step
if (QupstrUnit == "mm") {
iConvBasins <- which(!is.na(BasinAreas[seq.int(length(LengthHydro))]))
Qupstream[, iConvBasins] <- Qupstream[, iConvBasins] * rep(BasinAreas[iConvBasins], each = LLL) * 1e3
} else if (QupstrUnit == "m3/s") {
Qupstream <- Qupstream * TimeStep
} else if (QupstrUnit == "l/s") {
Qupstream <- Qupstream * TimeStep / 1e3
}
InputsModel <- c(InputsModel, list(Qupstream = Qupstream,
LengthHydro = LengthHydro,
BasinAreas = BasinAreas))
}
class(InputsModel) <- c("InputsModel",ObjectClass); class(InputsModel) <- c("InputsModel", ObjectClass)
return(InputsModel);
return(InputsModel)
}
}
R/CreateRunOptions.R
View file @ 4bccfee2
This diff is collapsed.
R/DataAltiExtrapolation_Valery.R
View file @ 4bccfee2
This diff is collapsed.