Calibration_Michel.R 17.1 KB
Newer Older
1
Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions,
2
                               FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
3
4
5
6
  
  
  ##_____Arguments_check_____________________________________________________________________
  if (!inherits(InputsModel, "InputsModel")) {
7
    stop("InputsModel must be of class 'InputsModel'")
8
9
10
    return(NULL)
  }  
  if (!inherits(RunOptions, "RunOptions")) {
11
    stop("RunOptions must be of class 'RunOptions'")
12
13
14
    return(NULL)
  }  
  if (!inherits(InputsCrit, "InputsCrit")) {
15
    stop("InputsCrit must be of class 'InputsCrit'")
16
    return(NULL)
17
18
  }
  if (inherits(InputsCrit, "Multi")) {
19
    stop("InputsCrit must be of class 'Single' or 'Compo'")
20
21
    return(NULL)
  }
22
  if (!inherits(CalibOptions, "CalibOptions")) {
23
    stop("CalibOptions must be of class 'CalibOptions'")
24
25
26
    return(NULL)
  }  
  if (!inherits(CalibOptions, "HBAN")) {
27
    stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used")
28
    return(NULL)
29
30
31
  }
  if (!missing(FUN_CRIT)) {
    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object", call. = FALSE)
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
  
  
  ##_check_FUN_TRANSFO
  if (is.null(FUN_TRANSFO)) {
    if (identical(FUN_MOD, RunModel_GR4H         )) {
      FUN_TRANSFO <- TransfoParam_GR4H
    }
    if (identical(FUN_MOD, RunModel_GR4J         )) {
      FUN_TRANSFO <- TransfoParam_GR4J
    }
    if (identical(FUN_MOD, RunModel_GR5J         )) {
      FUN_TRANSFO <- TransfoParam_GR5J
    }
    if (identical(FUN_MOD, RunModel_GR6J         )) {
      FUN_TRANSFO <- TransfoParam_GR6J
    }
    if (identical(FUN_MOD, RunModel_GR2M         )) {
      FUN_TRANSFO <- TransfoParam_GR2M
    }
    if (identical(FUN_MOD, RunModel_GR1A         )) {
      FUN_TRANSFO <- TransfoParam_GR1A
    }
    if (identical(FUN_MOD, RunModel_CemaNeige    )) {
57
58
59
60
61
      if (inherits(FUN_MOD, "hysteresis")) {
        FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
      } else {
        FUN_TRANSFO <- TransfoParam_CemaNeige
      }
62
63
64
65
    }
    if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
      if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
        FUN1 <- TransfoParam_GR4J
66
      }
67
68
      if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
        FUN1 <- TransfoParam_GR5J
69
      }
70
71
      if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
        FUN1 <- TransfoParam_GR6J
72
73
74
75
      }
      if (inherits(FUN_MOD, "hysteresis")) {
        FUN2 <- TransfoParam_CemaNeigeHyst
      } else {
76
        FUN2 <- TransfoParam_CemaNeige
77
      }
78
79
80
81
      FUN_TRANSFO <- function(ParamIn, Direction) {
        Bool <- is.matrix(ParamIn)
        if (Bool == FALSE) {
          ParamIn <- rbind(ParamIn)
82
        }
83
84
85
86
87
88
        ParamOut <- NA * ParamIn
        NParam   <- ncol(ParamIn)
        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, ]
89
        }
90
        return(ParamOut)
91
92
      }
    }
93
    if (is.null(FUN_TRANSFO)) {
94
      stop("FUN_TRANSFO was not found (in Calibration function)")
95
96
      return(NULL)
    }
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
  }
  
  ##_variables_initialisation 
  ParamFinalR <- NULL
  ParamFinalT <- NULL
  CritFinal   <- NULL
  NRuns <- 0
  NIter <- 0
  if ("StartParamDistrib" %in% names(CalibOptions)) {
    PrefilteringType <- 2
  } else {
    PrefilteringType <- 1
  }
  if (PrefilteringType == 1) {
    NParam <- ncol(CalibOptions$StartParamList)
  }
  if (PrefilteringType == 2) {
    NParam <- ncol(CalibOptions$StartParamDistrib)
  }
  if (NParam > 20) {
117
    stop("Calibration_Michel can handle a maximum of 20 parameters")
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
    return(NULL)
  }
  HistParamR    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
  HistParamT    <- matrix(NA, nrow = 500 * NParam, ncol = NParam)
  HistCrit      <- matrix(NA, nrow = 500 * NParam, ncol = 1)
  CritName      <- NULL
  CritBestValue <- NULL
  Multiplier    <- NULL
  CritOptim     <- +1e100
  ##_temporary_change_of_Outputs_Sim
  RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal  ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
  
  
  
  ##_____Parameter_Grid_Screening____________________________________________________________
  
  
  ##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
  ProposeCandidatesGrid <- function(DistribParam) {
    NewCandidates <- expand.grid(lapply(seq_len(ncol(DistribParamR)), function(x) DistribParam[, x]))
    NewCandidates <- unique(NewCandidates) # to avoid duplicates when a parameter is set
    Output <- list(NewCandidates = NewCandidates)
  }    
  
  
  ##Creation_of_new_candidates_______________________________________________
  OptimParam <- is.na(CalibOptions$FixedParam)
  if (PrefilteringType == 1) {
    CandidatesParamR <- CalibOptions$StartParamList
  }
  if (PrefilteringType == 2) {
    DistribParamR <- CalibOptions$StartParamDistrib
    DistribParamR[, !OptimParam] <- NA
    CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
  }
  ##Remplacement_of_non_optimised_values_____________________________________
  CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
    x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
    return(x)
  })
  if (NParam > 1) {
    CandidatesParamR <- t(CandidatesParamR)
  } else {
    CandidatesParamR <- cbind(CandidatesParamR)
  }
  
  ##Loop_to_test_the_various_candidates______________________________________
  iNewOptim <- 0
  Ncandidates <- nrow(CandidatesParamR)
  if (verbose & Ncandidates > 1) {
168
    if (PrefilteringType == 1) {
169
      message("List-Screening in progress (", appendLF = FALSE)
170
171
    }
    if (PrefilteringType == 2) {
172
      message("Grid-Screening in progress (", appendLF = FALSE)
173
    }
174
175
176
    message("0%", appendLF = FALSE)
  }
  for (iNew in 1:nrow(CandidatesParamR)) {
177
    if (verbose & Ncandidates > 1) {
178
179
180
      for (k in c(2, 4, 6, 8)) {
        if (iNew == round(k / 10 * Ncandidates)) {
          message(" ", 10 * k, "%", appendLF = FALSE)
181
        }
182
183
184
185
186
187
      } 
    }
    ##Model_run
    Param <- CandidatesParamR[iNew, ]
    OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
    ##Calibration_criterion_computation
188
    OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)      
189
190
191
192
    if (!is.na(OutputsCrit$CritValue)) {
      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
        iNewOptim <- iNew
Delaigue Olivier's avatar
Delaigue Olivier committed
193
194
      }
    }
195
196
197
198
199
    ##Storage_of_crit_info
    if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
      CritName      <- OutputsCrit$CritName
      CritBestValue <- OutputsCrit$CritBestValue
      Multiplier    <- OutputsCrit$Multiplier
200
    }
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
  }
  if (verbose & Ncandidates > 1) {
    message(" 100%)\n", appendLF = FALSE)
  }
  
  
  ##End_of_first_step_Parameter_Screening____________________________________
  ParamStartR <- CandidatesParamR[iNewOptim, ]
  if (!is.matrix(ParamStartR)) {
    ParamStartR <- matrix(ParamStartR, nrow = 1)
  }
  ParamStartT <- FUN_TRANSFO(ParamStartR, "RT")
  CritStart   <- CritOptim
  NRuns       <- NRuns+nrow(CandidatesParamR)
  if (verbose) {
    if (Ncandidates > 1) {
      message(sprintf("\t Screening completed (%s runs)", NRuns))
218
    }
219
220
    if (Ncandidates == 1) {
      message("\t Starting point for steepest-descent local search:")
Delaigue Olivier's avatar
Delaigue Olivier committed
221
    }
222
    message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
223
    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritStart * Multiplier), "\n")
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
  }
  ##Results_archiving________________________________________________________
  HistParamR[1, ] <- ParamStartR
  HistParamT[1, ] <- ParamStartT
  HistCrit[1, ]   <- CritStart
  
  
  
  
  ##_____Steepest_Descent_Local_Search_______________________________________________________
  
  
  ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
  ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
    ##Format_checking
    if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
240
      stop("each input set must be a matrix of one single line")
241
242
243
      return(NULL)
    }
    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
244
      stop("each input set must have the same number of values")
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
      return(NULL)
    }
    ##Proposal_of_new_parameter_sets ###(local search providing 2 * NParam-1 new sets)
    NParam <- ncol(NewParamOptimT)
    VECT <- NULL
    for (I in 1:NParam) {
      ##We_check_that_the_current_parameter_should_indeed_be_optimised
      if (OptimParam[I] == TRUE) {
        for (J in 1:2) {
          Sign <- 2 * J - 3   #Sign can be equal to -1 or +1
          ##We_define_the_new_potential_candidate
          Add <- TRUE
          PotentialCandidateT <- NewParamOptimT
          PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
          ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
          if (PotentialCandidateT[1, I] < RangesT[1, I] ) {
            PotentialCandidateT[1,I] <- RangesT[1, I]
          }
          if (PotentialCandidateT[1, I] > RangesT[2, I]) {
            PotentialCandidateT[1,I] <- RangesT[2,I]
          }
          ##We_check_the_set_is_not_outside_the_range_of_possible_values
          if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
            Add <- FALSE
          }
          if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
            Add <- FALSE
          }
          ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
          if (identical(PotentialCandidateT, OldParamOptimT)) {
            Add <- FALSE
          }
          ##We_add_the_candidate_to_our_list
          if (Add) {
            VECT <- c(VECT, PotentialCandidateT)
Delaigue Olivier's avatar
Delaigue Olivier committed
280
281
282
283
          }
        }
      }
    }
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
    Output <- NULL
    Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
    return(Output)
  }
  
  
  ##Initialisation_of_variables
  if (verbose) {
    message("Steepest-descent local search in progress") 
  }
  Pace <- 0.64
  PaceDiag <- rep(0, NParam)
  CLG <- 0.7^(1 / NParam)
  Compt <- 0
  CritOptim <- CritStart
  ##Conversion_of_real_parameter_values
  RangesR <- CalibOptions$SearchRanges
  RangesT <- FUN_TRANSFO(RangesR, "RT")
  NewParamOptimT <- ParamStartT
  OldParamOptimT <- ParamStartT
  
  
  ##START_LOOP_ITER_________________________________________________________
  for (ITER in 1:(100 * NParam)) {
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
310
    ##Exit_loop_when_Pace_becomes_too_small___________________________________
311
312
313
    if (Pace < 0.01) {
      break
    }
314
315
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
316
    ##Creation_of_new_candidates______________________________________________
317
318
    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
Delaigue Olivier's avatar
Delaigue Olivier committed
319
    ##Remplacement_of_non_optimised_values_____________________________________
320
321
    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
322
323
      return(x)
    })
324
325
326
327
328
    if (NParam > 1) {
      CandidatesParamR <- t(CandidatesParamR)
    } else {
      CandidatesParamR <- cbind(CandidatesParamR)
    }
329
330
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
331
    ##Loop_to_test_the_various_candidates_____________________________________
332
333
    iNewOptim <- 0
    for (iNew in 1:nrow(CandidatesParamR)) {
Delaigue Olivier's avatar
Delaigue Olivier committed
334
      ##Model_run
335
336
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Delaigue Olivier's avatar
Delaigue Olivier committed
337
      ##Calibration_criterion_computation
338
      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)      
339
      if (!is.na(OutputsCrit$CritValue)) {
340
341
        if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
          CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
342
343
344
          iNewOptim <- iNew
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
345
    }
346
    NRuns <- NRuns + nrow(CandidatesParamR)
347
348
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
349
    ##When_a_progress_has_been_achieved_______________________________________
350
    if (iNewOptim != 0) {
Delaigue Olivier's avatar
Delaigue Olivier committed
351
      ##We_store_the_optimal_set
352
353
      OldParamOptimT <- NewParamOptimT
      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
354
      Compt <- Compt + 1
Delaigue Olivier's avatar
Delaigue Olivier committed
355
      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
356
      if (Compt > 2 * NParam) {
357
358
        Pace <- Pace * 2
        Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
359
360
      }
      ##We_update_PaceDiag
361
362
363
      VectPace <- NewParamOptimT-OldParamOptimT
      for (iC in 1:NParam) {
        if (OptimParam[iC]) { 
364
          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
365
366
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
367
    } else {
368
      ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
369
370
      Pace <- Pace / 2
      Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
371
    }
372
373
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
374
    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
375
    if (ITER > 4 * NParam) {
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
      NRuns <- NRuns + 1
      iNewOptim <- 0
      iNew <- 1
      CandidatesParamT <- NewParamOptimT+PaceDiag
      if (!is.matrix(CandidatesParamT)) {
        CandidatesParamT <- matrix(CandidatesParamT, nrow = 1)
      }
      ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
      for (iC in 1:NParam) {
        if (OptimParam[iC]) {
          if (CandidatesParamT[iNew, iC] < RangesT[1, iC]) {
            CandidatesParamT[iNew, iC] <- RangesT[1, iC]
          }
          if (CandidatesParamT[iNew, iC] > RangesT[2, iC]) {
            CandidatesParamT[iNew, iC] <- RangesT[2, iC]
          }
Delaigue Olivier's avatar
Delaigue Olivier committed
392
        }
393
394
395
396
397
398
      }
      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
      ##Model_run
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
      ##Calibration_criterion_computation
399
      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
400
401
      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
402
403
404
405
406
407
408
        iNewOptim <- iNew
      }
      ##When_a_progress_has_been_achieved
      if (iNewOptim != 0) {
        OldParamOptimT <- NewParamOptimT
        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
      }
409
      
Delaigue Olivier's avatar
Delaigue Olivier committed
410
411
    }
    
412
    
Delaigue Olivier's avatar
Delaigue Olivier committed
413
    ##Results_archiving_______________________________________________________
414
415
416
417
418
    NewParamOptimR       <- FUN_TRANSFO(NewParamOptimT, "TR")
    HistParamR[ITER+1, ] <- NewParamOptimR
    HistParamT[ITER+1, ] <- NewParamOptimT
    HistCrit[ITER+1, ]   <- CritOptim
    ### if (verbose) { cat(paste("\t     Iter ",formatC(ITER,format="d",width=3), "    Crit ",formatC(CritOptim,format="f",digits=4), "    Pace ",formatC(Pace,format="f",digits=4), "\n",sep=""))}
Delaigue Olivier's avatar
Delaigue Olivier committed
419
420
    
    
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
    
  } ##END_LOOP_ITER_________________________________________________________
  ITER <- ITER - 1
  
  
  ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
  if (CritOptim == CritStart & verbose) { 
    message("\t No progress achieved")
  }
  
  ##End_of_Steepest_Descent_Local_Search____________________________________
  ParamFinalR <- NewParamOptimR
  ParamFinalT <- NewParamOptimT
  CritFinal   <- CritOptim
  NIter       <- 1 + ITER
  if (verbose) { 
    message(sprintf("\t Calibration completed (%s iterations, %s runs)", NIter, NRuns))
    message("\t     Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = " , "))
439
440
441
442
443
444
445
446
447
448
449
    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritFinal * Multiplier), "\n")
    if (inherits(InputsCrit, "Compo")) {
      listweights  <- OutputsCrit$CritCompo$MultiCritWeights
      listNameCrit <- OutputsCrit$CritCompo$MultiCritNames
      msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
      msgForm <- unlist(strsplit(msgForm, split = ","))
      msgFormSep <- rep(c(",", ",", ",\n\t\t     "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
      msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
      msgForm <- gsub("\\,\\\n\\\t\\\t     $|\\,$", "", msgForm)
      message("\tFormula: mean(", msgForm, ")\n")
    }
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
  }
  ##Results_archiving_______________________________________________________
  HistParamR <- cbind(HistParamR[1:NIter, ])
  colnames(HistParamR) <- paste0("Param", 1:NParam)
  HistParamT <- cbind(HistParamT[1:NIter, ])
  colnames(HistParamT) <- paste0("Param", 1:NParam)
  HistCrit   <- cbind(HistCrit[1:NIter, ])
  ###colnames(HistCrit) <- paste("HistCrit")
  
  BoolCrit_Actual <- InputsCrit$BoolCrit
  BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
  MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
  colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
  
  
  ##_____Output______________________________________________________________________________
  OutputsCalib <- list(ParamFinalR = as.double(ParamFinalR), CritFinal = CritFinal * Multiplier,
                       NIter = NIter, NRuns = NRuns,
                       HistParamR = HistParamR, HistCrit = HistCrit * Multiplier,
                       MatBoolCrit = MatBoolCrit,
                       CritName = CritName, CritBestValue = CritBestValue)
  class(OutputsCalib) <- c("OutputsCalib", "HBAN")
  return(OutputsCalib)
  
  
  
Delaigue Olivier's avatar
Delaigue Olivier committed
476
477
478
479
}



480

Delaigue Olivier's avatar
Delaigue Olivier committed
481