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\.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
.Rapp.history
# Session Data files
.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
Type: Package
Title: Suite of GR hydrological models for precipitation-runoff modelling
Version: 1.0.1
Date: 2016-04-21
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.7.6.9000
Date: 2023-10-25
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "cre", "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("ctb")),
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"))
)
Author: Laurent Coron, Charles Perrin, with contributions
from Vazken Andréassian, Pierre Brigode, Olivier Delaigue,
Nicolas Le Moine, Thibaut Mathevet, Safouane Mouelhi,
Ludovic Oudin, Raji Pushpalatha, Guillaume Thirel, Audrey Valéry.
Based on earlier work by Claude Michel.
Maintainer: Olivier Delaigue <airGR@irstea.fr>
Depends: R (>= 3.0.1)
Description: This package brings into R the 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 snowmelt module (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: http://webgr.irstea.fr/modeles/?lang=en
URL: https://hydrogr.github.io/airGR/
BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues
NeedsCompilation: yes
Encoding: UTF-8
VignetteBuilder: knitr
RoxygenNote: 7.1.1
NAMESPACE
View file @ 4bccfee2
#####################################
## 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)
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
export(CreateErrorCrit_GAPX)
export(CreateIniStates)
export(CreateInputsCrit)
export(CreateInputsCrit_Lavenne)
export(CreateInputsModel)
export(CreateRunOptions)
export(DataAltiExtrapolation_Valery)
......@@ -20,26 +36,44 @@ 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(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){
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/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,
StartParam=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_OptimParam
if(is.null(OptimParam)){
OptimParam <- rep(TRUE,NParam);
} 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); }
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)
##Create_CalibOptions
CalibOptions <- list(OptimParam=OptimParam,FixedParam=FixedParam,SearchRanges=SearchRanges);
if(!is.null(StartParam )){ CalibOptions <- c(CalibOptions,list(StartParam=StartParam)); }
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);
} 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 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
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); }
ObjectClass <- 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); }
## ---------- check arguments
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); }
## check 'InputsModel'
if (!inherits(InputsModel, "InputsModel")) {
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)){
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); }
## 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(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);
}
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_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
BOOL <- FALSE;
if(identical(FUN_MOD,RunModel_GR4H)){
ObjectClass <- c(ObjectClass,"hourly","GR");
TimeStep <- as.integer(60*60);
BOOL <- TRUE;
ObjectClass <- NULL
## check DatesR
if (is.null(DatesR)) {
stop("'DatesR' is missing")
}
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)){
ObjectClass <- c(ObjectClass,"daily","GR");
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
if (!is.numeric(Precip) | !is.numeric(PotEvap)) {
stop("'Precip' and 'PotEvap' must be vectors of numeric values")
}
if(identical(FUN_MOD,RunModel_GR2M)){
ObjectClass <- c(ObjectClass,"GR","monthly");
TimeStep <- as.integer(c(28,29,30,31)*24*60*60);
BOOL <- TRUE;
if (length(Precip) != LLL | length(PotEvap) != LLL) {
stop("'Precip', 'PotEvap' and 'DatesR' must have the same length")
}
if(identical(FUN_MOD,RunModel_GR1A)){
ObjectClass <- c(ObjectClass,"GR","yearly");
TimeStep <- as.integer(c(365,366)*24*60*60);
BOOL <- TRUE;
}
if ("CemaNeige" %in% ObjectClass) {
if (is.null(Precip)) {
stop("'Precip' is missing")
}
if(identical(FUN_MOD,RunModel_CemaNeige)){
ObjectClass <- c(ObjectClass,"daily","CemaNeige");
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
if (is.null(TempMean)) {
stop("'TempMean' is missing")
}
if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
ObjectClass <- c(ObjectClass,"daily","GR","CemaNeige");
TimeStep <- as.integer(24*60*60);
BOOL <- TRUE;
if (!is.vector(Precip) | !is.vector(TempMean)) {
stop("'Precip' and 'TempMean' must be vectors of numeric values")
}
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(DatesR)){ stop("DatesR is missing \n"); return(NULL); }
if("POSIXlt" %in% class(DatesR) == FALSE & "POSIXct" %in% class(DatesR) == FALSE){ stop("DatesR must be defined as POSIXlt or POSIXct \n"); return(NULL); }
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];
}
if (is.null(ZInputs)) {
if (verbose & !identical(HypsoData, as.numeric(rep(NA, 101)))) {
warning("'ZInputs' is missing: HypsoData[51] is used")
}
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
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"); } }
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"); } }
}
if("CemaNeige" %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"); } }
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"); } }
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"); } } }
}
if(sum(BOOL_NA)!=0){
WTxt <- NULL;
WTxt <- paste(WTxt,"\t Missing values are not allowed in InputsModel \n",sep="");
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){
Precip <- Precip[Select]; TempMean <- TempMean[Select]; if(!is.null(TempMin) & !is.null(TempMax)){ TempMin <- TempMin[Select]; TempMax <- TempMax[Select]; } }
DatesR <- DatesR[Select];
WTxt <- paste(WTxt,"\t -> data were trunced to keep the most recent available time-steps \n",sep="");
WTxt <- paste(WTxt,"\t -> ",length(Select)," time-steps were kept \n",sep="");
if(!is.null(WTxt) & verbose){ warning(WTxt); }
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("Values < 0 or NA values detected in 'Precip' series")
}
}
BOOL_NA_TMP <- (PotEvap < 0) | is.na(PotEvap)
if (sum(BOOL_NA_TMP) != 0) {
BOOL_NA <- BOOL_NA | BOOL_NA_TMP
if (verbose) {
warning("Values < 0 or NA values detected in 'PotEvap' series")
}
}
}
if ("CemaNeige" %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("Values < 0 or NA values detected in 'Precip' series")
}
}
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
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);
if(verbose){ if(NLayers==1){ cat(paste("\t Input series were successfully created on 1 elevation layer for use by CemaNeige \n",sep=""));
} else { cat(paste("\t Input series were successfully created on ",NLayers," elevation layers for use by CemaNeige \n",sep="")); } }
if ("CemaNeige" %in% ObjectClass) {
RESULT <- DataAltiExtrapolation_Valery(DatesR = DatesR,
Precip = Precip, PrecipScale = PrecipScale,
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
InputsModel <- list(DatesR=DatesR);
if("GR" %in% ObjectClass){
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,
LayerFracSolidPrecip=RESULT$LayerFracSolidPrecip,ZLayers=RESULT$ZLayers)); }
InputsModel <- list(DatesR = DatesR)
if ("GR" %in% ObjectClass) {
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,
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);
return(InputsModel);
class(InputsModel) <- c("InputsModel", ObjectClass)
return(InputsModel)
}
}
R/CreateRunOptions.R
View file @ 4bccfee2
This diff is collapsed.
R/DataAltiExtrapolation_Valery.R
View file @ 4bccfee2
This diff is collapsed.