Calibration_Michel.R 16 KB
Newer Older
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
Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions,
                               FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
  
  
  ##_____Arguments_check_____________________________________________________________________
  if (!inherits(InputsModel, "InputsModel")) {
    stop("InputsModel must be of class 'InputsModel' \n")
    return(NULL)
  }  
  if (!inherits(RunOptions, "RunOptions")) {
    stop("RunOptions must be of class 'RunOptions' \n")
    return(NULL)
  }  
  if (!inherits(InputsCrit, "InputsCrit")) {
    stop("InputsCrit must be of class 'InputsCrit' \n")
    return(NULL)
  }  
  if (!inherits(CalibOptions, "CalibOptions")) {
    stop("CalibOptions must be of class 'CalibOptions' \n")
    return(NULL)
  }  
  if (!inherits(CalibOptions, "HBAN")) {
    stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \n")
    return(NULL)
  }  
  
  
  ##_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    )) {
      FUN_TRANSFO <- TransfoParam_CemaNeige
    }
    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
        FUN2 <- TransfoParam_CemaNeige
55
      }
56
57
58
      if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
        FUN1 <- TransfoParam_GR5J
        FUN2 <- TransfoParam_CemaNeige
59
      }
60
61
62
      if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
        FUN1 <- TransfoParam_GR6J
        FUN2 <- TransfoParam_CemaNeige
63
      }
64
65
66
67
      FUN_TRANSFO <- function(ParamIn, Direction) {
        Bool <- is.matrix(ParamIn)
        if (Bool == FALSE) {
          ParamIn <- rbind(ParamIn)
68
        }
69
70
71
72
73
74
        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, ]
75
        }
76
        return(ParamOut)
77
78
      }
    }
79
80
    if (is.null(FUN_TRANSFO)) {
      stop("FUN_TRANSFO was not found (in Calibration function) \n")
81
82
      return(NULL)
    }
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
  }
  
  ##_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) {
    stop("Calibration_Michel can handle a maximum of 20 parameters \n")
    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) {
154
    if (PrefilteringType == 1) {
155
      message("List-Screening in progress (", appendLF = FALSE)
156
157
    }
    if (PrefilteringType == 2) {
158
      message("Grid-Screening in progress (", appendLF = FALSE)
159
    }
160
161
162
    message("0%", appendLF = FALSE)
  }
  for (iNew in 1:nrow(CandidatesParamR)) {
163
    if (verbose & Ncandidates > 1) {
164
165
166
      for (k in c(2, 4, 6, 8)) {
        if (iNew == round(k / 10 * Ncandidates)) {
          message(" ", 10 * k, "%", appendLF = FALSE)
167
        }
168
169
170
171
172
173
174
175
176
177
178
      } 
    }
    ##Model_run
    Param <- CandidatesParamR[iNew, ]
    OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
    ##Calibration_criterion_computation
    OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)      
    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
179
180
      }
    }
181
182
183
184
185
    ##Storage_of_crit_info
    if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
      CritName      <- OutputsCrit$CritName
      CritBestValue <- OutputsCrit$CritBestValue
      Multiplier    <- OutputsCrit$Multiplier
186
    }
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  }
  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))
204
    }
205
206
    if (Ncandidates == 1) {
      message("\t Starting point for steepest-descent local search:")
Delaigue Olivier's avatar
Delaigue Olivier committed
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
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
    message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
    message(sprintf("\t     Crit %-12s = %.4f", CritName, CritStart * Multiplier))
  }
  ##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) {
      stop("each input set must be a matrix of one single line \n")
      return(NULL)
    }
    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
      stop("each input set must have the same number of values \n")
      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
266
267
268
269
          }
        }
      }
    }
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
    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
296
    ##Exit_loop_when_Pace_becomes_too_small___________________________________
297
298
299
    if (Pace < 0.01) {
      break
    }
300
301
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
302
    ##Creation_of_new_candidates______________________________________________
303
304
    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
Delaigue Olivier's avatar
Delaigue Olivier committed
305
    ##Remplacement_of_non_optimised_values_____________________________________
306
307
    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
308
309
      return(x)
    })
310
311
312
313
314
    if (NParam > 1) {
      CandidatesParamR <- t(CandidatesParamR)
    } else {
      CandidatesParamR <- cbind(CandidatesParamR)
    }
315
316
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
317
    ##Loop_to_test_the_various_candidates_____________________________________
318
319
    iNewOptim <- 0
    for (iNew in 1:nrow(CandidatesParamR)) {
Delaigue Olivier's avatar
Delaigue Olivier committed
320
      ##Model_run
321
322
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Delaigue Olivier's avatar
Delaigue Olivier committed
323
      ##Calibration_criterion_computation
324
325
      OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)      
      if (!is.na(OutputsCrit$CritValue)) {
326
327
        if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
          CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
328
329
330
          iNewOptim <- iNew
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
331
    }
332
    NRuns <- NRuns + nrow(CandidatesParamR)
333
334
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
335
    ##When_a_progress_has_been_achieved_______________________________________
336
    if (iNewOptim != 0) {
Delaigue Olivier's avatar
Delaigue Olivier committed
337
      ##We_store_the_optimal_set
338
339
      OldParamOptimT <- NewParamOptimT
      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
340
      Compt <- Compt + 1
Delaigue Olivier's avatar
Delaigue Olivier committed
341
      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
342
      if (Compt > 2 * NParam) {
343
344
        Pace <- Pace * 2
        Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
345
346
      }
      ##We_update_PaceDiag
347
348
349
      VectPace <- NewParamOptimT-OldParamOptimT
      for (iC in 1:NParam) {
        if (OptimParam[iC]) { 
350
          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
351
352
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
353
    } else {
354
      ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
355
356
      Pace <- Pace / 2
      Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
357
    }
358
359
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
360
    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
361
    if (ITER > 4 * NParam) {
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
      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
378
        }
379
380
381
382
383
384
      }
      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
      ##Model_run
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
      ##Calibration_criterion_computation
385
      OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)
386
387
      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
388
389
390
391
392
393
394
        iNewOptim <- iNew
      }
      ##When_a_progress_has_been_achieved
      if (iNewOptim != 0) {
        OldParamOptimT <- NewParamOptimT
        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
      }
395
      
Delaigue Olivier's avatar
Delaigue Olivier committed
396
397
    }
    
398
    
Delaigue Olivier's avatar
Delaigue Olivier committed
399
    ##Results_archiving_______________________________________________________
400
401
402
403
404
    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
405
406
    
    
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
    
  } ##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 = " , "))
    message(sprintf("\t     Crit %-12s = %.4f", CritName, CritFinal * Multiplier))
  }
  ##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
452
453
454
455
}



456

Delaigue Olivier's avatar
Delaigue Olivier committed
457