CreateCalibOptions.R 16 KB
Newer Older
Delaigue Olivier's avatar
Delaigue Olivier committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
#*************************************************************************************************
#' Creation of the CalibOptions object required to the Calibration functions.
#'
#' Users wanting to use FUN_MOD, FUN_CALIB or FUN_TRANSFO functions that are not included in 
#' the package must create their own CalibOptions object accordingly.
#*************************************************************************************************
#' @title  Creation of the CalibOptions object required to the Calibration functions
#' @author Laurent Coron (June 2014)
#' @seealso \code{\link{Calibration}}, \code{\link{RunModel}}
#' @example tests/example_Calibration.R
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________
#' @param  FUN_MOD             [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)
#' @param  FUN_CALIB           (optional) [function] calibration algorithm function (e.g. Calibration_HBAN, Calibration_optim), default=Calibration_HBAN
#' @param  FUN_TRANSFO         (optional) [function] model parameters transformation function, if the FUN_MOD used is native in the package FUN_TRANSFO is automatically defined
#' @param  OptimParam          (optional) [boolean] vector of booleans indicating which parameters must be optimised (NParam columns, 1 line) 
#' @param  FixedParam          (optional) [numeric] vector giving the values to allocate to non-optimised parameter values (NParam columns, 1 line)
#' @param  SearchRanges        (optional) [numeric] matrix giving the ranges of real parameters (NParam columns, 2 lines)
#'                             \tabular{llllll}{
#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
#'                                  [1,]     \tab    0   \tab   -1   \tab    0   \tab  ...    \tab  0.0 \cr
#'                                  [2,]     \tab 3000   \tab   +1   \tab  100   \tab  ...    \tab  3.0 \cr
#'                             }
#' @param  StartParam          (optional) [numeric] vector of parameter values used to start global search calibration procedure (e.g. Calibration_optim)
#'                             \tabular{llllll}{
#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
#'                                           \tab 1000   \tab -0.5   \tab   22   \tab  ...    \tab  1.1 \cr
#'                             }
#' @param  StartParamList      (optional) [numeric] matrix of parameter sets used for grid-screening calibration procedure (values in columns, sets in line)
#'                             \tabular{llllll}{
#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
#'                                  [set1]   \tab  800   \tab -0.7   \tab   25   \tab  ...    \tab  1.0 \cr
#'                                  [set2]   \tab 1000   \tab -0.5   \tab   22   \tab  ...    \tab  1.1 \cr
#'                                  [...]    \tab  ...   \tab  ...   \tab  ...   \tab  ...    \tab  ... \cr
#'                                  [set n]  \tab  200   \tab -0.3   \tab   17   \tab  ...    \tab  1.0 \cr
#'                             }
#' @param  StartParamDistrib   (optional) [numeric] matrix of parameter values used for grid-screening calibration procedure (values in columns, percentiles in line) \cr
#'                             \tabular{llllll}{
#'                                           \tab [X1]   \tab [X2]   \tab [X3]   \tab [...]   \tab [Xi] \cr
#'                                  [value1] \tab  800   \tab -0.7   \tab   25   \tab  ...    \tab  1.0 \cr
#'                                  [value2] \tab 1000   \tab   NA   \tab   50   \tab  ...    \tab  1.2 \cr
#'                                  [value3] \tab 1200   \tab   NA   \tab   NA   \tab  ...    \tab  1.6 \cr
#'                             }
#_FunctionOutputs_________________________________________________________________________________
#' @return  [list] object of class \emph{CalibOptions} containing the data required to evaluate the model outputs; it can include the following:
#'          \tabular{ll}{
#'          \emph{$OptimParam       }  \tab   [boolean] vector of booleans indicating which parameters must be optimised \cr
#'          \emph{$FixedParam       }  \tab   [numeric] vector giving the values to allocate to non-optimised parameter values \cr
#'          \emph{$SearchRanges     }  \tab   [numeric] matrix giving the ranges of real parameters \cr
#'          \emph{$StartParam       }  \tab   [numeric] vector of parameter values used to start global search calibration procedure \cr
#'          \emph{$StartParamList   }  \tab   [numeric] matrix of parameter sets used for grid-screening calibration procedure \cr
#'          \emph{$StartParamDistrib}  \tab   [numeric] matrix of parameter values used for grid-screening calibration procedure \cr
#'          }
#**************************************************************************************************
CreateCalibOptions <- function(FUN_MOD,FUN_CALIB=Calibration_HBAN,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_HBAN  )){ ObjectClass <- c(ObjectClass,"HBAN"  ); BOOL <- TRUE; }
    if(identical(FUN_CALIB,Calibration_optim )){ ObjectClass <- c(ObjectClass,"optim" ); 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)) |
        ("optim" %in% ObjectClass & is.null(StartParam)) ){

      if("GR4H"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -2.00, +3.40, -9.10,
                             +3.90, -0.90, +4.10, -8.70,
                             +4.50, -0.10, +5.00, -8.10),ncol=NParam,byrow=TRUE); }
      if("GR4J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -2.00, +3.40, -9.10,
                             +3.90, -0.90, +4.10, -8.70,
                             +4.50, -0.10, +5.00, -8.10),ncol=NParam,byrow=TRUE); }
      if("GR5J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -1.70, +3.30, -9.10, -0.70,
                             +3.90, -0.60, +4.10, -8.70, +0.30,
                             +4.50, -0.10, +5.00, -8.10, +0.50),ncol=NParam,byrow=TRUE); }
      if("GR6J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
                             +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
                             +4.50, +0.50, +5.00, -8.10, +1.10, +5.00),ncol=NParam,byrow=TRUE); }
      if("GR2M"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -5.00,
                             +3.90, +0.00,
                             +4.50, +5.00),ncol=NParam,byrow=TRUE); }
      if("GR1A"%in% ObjectClass){ 
        ParamT <- matrix( c( -5.00,
                             +0.00,
                             +5.00),ncol=NParam,byrow=TRUE); }
      if("CemaNeige"%in% ObjectClass){ 
        ParamT <- matrix( c( -6.26, +0.55,
                             -2.13, +0.92,
                             +4.86, +1.40),ncol=NParam,byrow=TRUE); }
      if("CemaNeigeGR4J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -2.00, +3.40, -9.10, -6.26, +0.55,
                             +3.90, -0.90, +4.10, -8.70, -2.13, +0.92,
                             +4.50, -0.10, +5.00, -8.10, +4.86, +1.40),ncol=NParam,byrow=TRUE); }
      if("CemaNeigeGR5J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -1.70, +3.30, -9.10, -0.70, -6.26, +0.55,
                             +3.90, -0.60, +4.10, -8.70, +0.30, -2.13, +0.92,
                             +4.50, -0.10, +5.00, -8.10, +0.50, +4.86, +1.40),ncol=NParam,byrow=TRUE); }
      if("CemaNeigeGR6J"%in% ObjectClass){ 
        ParamT <- matrix( c( +3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -6.26, +0.55,
                             +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -2.13, +0.92,
                             +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.86, +1.40),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("optim" %in% ObjectClass & !is.null(StartParam)){
      if(!is.vector( StartParam)    ){ stop("StartParam must be a vector \n"); return(NULL);  }
      if(!is.numeric(StartParam)    ){ stop("StartParam must be a vector of numeric values \n"); return(NULL);  }
      if(sum(is.na(StartParam))!=0  ){ stop("StartParam must not include NA values \n"); return(NULL);  }
      if(length(StartParam)!=NParam ){ stop("Incompatibility between StartParam length and FUN_MOD \n"); return(NULL);  }
    }


  ##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);


}