Calibration_Michel.R 18.5 KB
Newer Older
1
Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE) {
Delaigue Olivier's avatar
Delaigue Olivier committed
2
3
4


##_____Arguments_check_____________________________________________________________________
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
    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)
      }  
Delaigue Olivier's avatar
Delaigue Olivier committed
25
26
27


   ##_check_FUN_TRANSFO
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
    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
Delaigue Olivier's avatar
Delaigue Olivier committed
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
        if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
          FUN1 <- TransfoParam_GR5J
          FUN2 <- TransfoParam_CemaNeige
        }
        if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
          FUN1 <- TransfoParam_GR6J
          FUN2 <- TransfoParam_CemaNeige
        }
        FUN_TRANSFO <- function(ParamIn, Direction) {
          Bool <- is.matrix(ParamIn)
          if (Bool == FALSE) {
            ParamIn <- rbind(ParamIn)
          }
          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, ]
          }
          return(ParamOut)
        }
      }
      if (is.null(FUN_TRANSFO)) {
        stop("FUN_TRANSFO was not found (in Calibration function) \n")
        return(NULL)
Delaigue Olivier's avatar
Delaigue Olivier committed
81
82
83
84
      }
    }

    ##_variables_initialisation 
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
    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
Delaigue Olivier's avatar
Delaigue Olivier committed
112
    ##_temporary_change_of_Outputs_Sim
113
    RunOptions$Outputs_Sim <- RunOptions$Outputs_Cal  ### this reduces the size of the matrix exchange with fortran and therefore speeds the calibration
Delaigue Olivier's avatar
Delaigue Olivier committed
114
115
116
117
118
119
120



##_____Parameter_Grid_Screening____________________________________________________________


    ##Definition_of_the_function_creating_all_possible_parameter_sets_from_different_values_for_each_parameter
121
    ProposeCandidatesGrid <- function(DistribParam) {
Delaigue Olivier's avatar
Delaigue Olivier committed
122
      ##Managing_matrix_sizes
123
124
        Nvalmax <- nrow(DistribParam)
        NParam <- ncol(DistribParam)
Delaigue Olivier's avatar
Delaigue Olivier committed
125
        ##we_add_columns_to_MatDistrib_until_it_has_20_columns
126
127
        DistribParam2 <- matrix(NA, nrow = Nvalmax, ncol = 20)
        DistribParam2[1:Nvalmax, 1:NParam] <- DistribParam
Delaigue Olivier's avatar
Delaigue Olivier committed
128
        ##we_check_the_number_of_values_to_test_for_each_param
129
130
131
132
        NbDistrib <- rep(1, 20)
        for (iC in 1:20) {
          NbDistrib[iC] <- max(1, Nvalmax-sum(is.na(DistribParam2[, iC])))
        }
Delaigue Olivier's avatar
Delaigue Olivier committed
133
134
      ##Loop_on_the_various_values_to_test ###(if 4 param and 3 values for each => 3^4 sets)
      ##NB_we_always_do_20_loops ###which_is_here_the_max_number_of_param_that_can_be_optimised
135
136
137
138
139
        VECT <- NULL
        for (iL01 in 1:NbDistrib[01]) { for (iL02 in 1:NbDistrib[02]) { for (iL03 in 1:NbDistrib[03]) { for (iL04 in 1:NbDistrib[04]) { for (iL05 in 1:NbDistrib[05]) { 
        for (iL06 in 1:NbDistrib[06]) { for (iL07 in 1:NbDistrib[07]) { for (iL08 in 1:NbDistrib[08]) { for (iL09 in 1:NbDistrib[09]) { for (iL10 in 1:NbDistrib[10]) {
        for (iL11 in 1:NbDistrib[11]) { for (iL12 in 1:NbDistrib[12]) { for (iL13 in 1:NbDistrib[13]) { for (iL14 in 1:NbDistrib[14]) { for (iL15 in 1:NbDistrib[15]) { 
        for (iL16 in 1:NbDistrib[16]) { for (iL17 in 1:NbDistrib[17]) { for (iL18 in 1:NbDistrib[18]) { for (iL19 in 1:NbDistrib[19]) { for (iL20 in 1:NbDistrib[20]) {
Delaigue Olivier's avatar
Delaigue Olivier committed
140
141
142
143
          VECT <- c(VECT,                         
            DistribParam2[iL01,01],DistribParam2[iL02,02],DistribParam2[iL03,03],DistribParam2[iL04,04],DistribParam2[iL05,05],
            DistribParam2[iL06,06],DistribParam2[iL07,07],DistribParam2[iL08,08],DistribParam2[iL09,09],DistribParam2[iL10,10],
            DistribParam2[iL11,11],DistribParam2[iL12,12],DistribParam2[iL13,13],DistribParam2[iL14,14],DistribParam2[iL15,15],
144
            DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20])
Delaigue Olivier's avatar
Delaigue Olivier committed
145
146
147
148
        } } } } }
        } } } } }
        } } } } }
        } } } } }
149
150
151
152
153
154
155
        MAT <- matrix(VECT, ncol = 20, byrow = TRUE)[, 1:NParam]
        if (!is.matrix(MAT)) {
          MAT <- cbind(MAT)
        }
        Output <- NULL
        Output$NewCandidates <- MAT
        return(Output)
Delaigue Olivier's avatar
Delaigue Olivier committed
156
157
158
159
    }


    ##Creation_of_new_candidates_______________________________________________
160
    OptimParam <- is.na(CalibOptions$FixedParam)
161
162
163
164
165
166
167
168
    if (PrefilteringType == 1) {
      CandidatesParamR <- CalibOptions$StartParamList
    }
    if (PrefilteringType == 2) {
      DistribParamR <- CalibOptions$StartParamDistrib
      DistribParamR[,!OptimParam] <- NA
      CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates
    }
Delaigue Olivier's avatar
Delaigue Olivier committed
169
    ##Remplacement_of_non_optimised_values_____________________________________
170
171
172
173
174
175
176
    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
      return(x)})
    if (NParam>1) {
      CandidatesParamR <- t(CandidatesParamR)
    } else { CandidatesParamR <- cbind(CandidatesParamR)
    }
Delaigue Olivier's avatar
Delaigue Olivier committed
177
178

    ##Loop_to_test_the_various_candidates______________________________________
179
180
181
182
183
184
185
186
187
    iNewOptim <- 0
    Ncandidates <- nrow(CandidatesParamR)
    if (verbose & Ncandidates > 1) {
      if (PrefilteringType == 1) {
        message("List-Screening in progress (", appendLF = FALSE)
      }
      if (PrefilteringType == 2) {
        message("Grid-Screening in progress (", appendLF = FALSE)
      }
188
      message("0%", appendLF = FALSE)
Delaigue Olivier's avatar
Delaigue Olivier committed
189
    }
190
191
192
193
194
195
196
    for (iNew in 1:nrow(CandidatesParamR)) {
      if (verbose & Ncandidates > 1) {
        for (k in c(2, 4, 6, 8)) {
          if (iNew == round(k/10*Ncandidates)) {
            message(" ", 10*k, "%", appendLF = FALSE)
          }
        } 
Delaigue Olivier's avatar
Delaigue Olivier committed
197
198
      }
      ##Model_run
199
200
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Delaigue Olivier's avatar
Delaigue Olivier committed
201
      ##Calibration_criterion_computation
202
203
204
205
206
207
208
      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
209
      ##Storage_of_crit_info
210
211
212
213
      if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
        CritName      <- OutputsCrit$CritName
        CritBestValue <- OutputsCrit$CritBestValue
        Multiplier    <- OutputsCrit$Multiplier
Delaigue Olivier's avatar
Delaigue Olivier committed
214
215
      }
    }
216
217
218
    if (verbose & Ncandidates > 1) {
      message(" 100%)\n", appendLF = FALSE)
    }
Delaigue Olivier's avatar
Delaigue Olivier committed
219
220
221
	  

    ##End_of_first_step_Parameter_Screening____________________________________
222
223
224
225
226
227
228
229
230
    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) {
231
        message(sprintf("\t Screening completed (%s runs)", NRuns))
232
      }
233
      if (Ncandidates == 1) {
234
235
        message("\t Starting point for steepest-descent local search:")
      }
236
237
      message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = " , "))
      message(sprintf("\t     Crit %-12s = %.4f", CritName, CritStart*Multiplier))
Delaigue Olivier's avatar
Delaigue Olivier committed
238
239
    }
    ##Results_archiving________________________________________________________
240
241
242
    HistParamR[1, ] <- ParamStartR
    HistParamT[1, ] <- ParamStartT
    HistCrit[1, ]   <- CritStart
Delaigue Olivier's avatar
Delaigue Olivier committed
243
244
245
246
247
248
249
250




##_____Steepest_Descent_Local_Search_______________________________________________________


    ##Definition_of_the_function_creating_new_parameter_sets_through_a_step_by_step_progression_procedure
251
    ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam,Pace) {
Delaigue Olivier's avatar
Delaigue Olivier committed
252
      ##Format_checking
253
254
255
256
257
258
259
260
      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)
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
261
      ##Proposal_of_new_parameter_sets ###(local search providing 2*NParam-1 new sets)
262
263
264
      NParam <- ncol(NewParamOptimT)
      VECT <- NULL
      for (I in 1:NParam) {
Delaigue Olivier's avatar
Delaigue Olivier committed
265
        ##We_check_that_the_current_parameter_should_indeed_be_optimised
266
267
268
        if (OptimParam[I] == TRUE) {
          for (J in 1:2) {
            Sign <- 2 * J - 3   #Sign can be equal to -1 or +1
Delaigue Olivier's avatar
Delaigue Olivier committed
269
            ##We_define_the_new_potential_candidate
270
271
272
            Add <- TRUE
            PotentialCandidateT <- NewParamOptimT
            PotentialCandidateT[1, I] <- NewParamOptimT[I] + Sign * Pace
Delaigue Olivier's avatar
Delaigue Olivier committed
273
            ##If_we_exit_the_range_of_possible_values_we_go_back_on_the_boundary
274
275
276
277
278
279
            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]
            }
Delaigue Olivier's avatar
Delaigue Olivier committed
280
            ##We_check_the_set_is_not_outside_the_range_of_possible_values
281
282
283
284
285
286
             if (NewParamOptimT[I] == RangesT[1, I] & Sign < 0) {
               Add <- FALSE
             }
             if (NewParamOptimT[I] == RangesT[2, I] & Sign > 0) {
               Add <- FALSE
             }
Delaigue Olivier's avatar
Delaigue Olivier committed
287
            ##We_check_that_this_set_has_not_been_tested_during_the_last_iteration
288
289
290
            if (identical(PotentialCandidateT, OldParamOptimT)) {
              Add <- FALSE
            }
Delaigue Olivier's avatar
Delaigue Olivier committed
291
            ##We_add_the_candidate_to_our_list
292
293
294
            if (Add == TRUE) {
              VECT <- c(VECT, PotentialCandidateT)
            }
Delaigue Olivier's avatar
Delaigue Olivier committed
295
296
297
          }
        }
      }
298
299
300
      Output <- NULL
      Output$NewCandidatesT <- matrix(VECT, ncol = NParam, byrow = TRUE)
      return(Output)
Delaigue Olivier's avatar
Delaigue Olivier committed
301
302
303
304
    }
      

    ##Initialisation_of_variables
305
    if (verbose) {
306
      message("Steepest-descent local search in progress") 
Delaigue Olivier's avatar
Delaigue Olivier committed
307
    }
308
309
310
311
312
    Pace <- 0.64
    PaceDiag <- rep(0, NParam)
    CLG <- 0.7^(1/NParam)
    Compt <- 0
    CritOptim <- CritStart
Delaigue Olivier's avatar
Delaigue Olivier committed
313
    ##Conversion_of_real_parameter_values
314
315
316
317
    RangesR <- CalibOptions$SearchRanges
    RangesT <- FUN_TRANSFO(RangesR, "RT")
    NewParamOptimT <- ParamStartT
    OldParamOptimT <- ParamStartT
Delaigue Olivier's avatar
Delaigue Olivier committed
318
319
320


    ##START_LOOP_ITER_________________________________________________________
321
    for (ITER in 1:(100*NParam)) {
Delaigue Olivier's avatar
Delaigue Olivier committed
322
323
324


    ##Exit_loop_when_Pace_becomes_too_small___________________________________
325
326
327
    if (Pace < 0.01) {
      break
    }
Delaigue Olivier's avatar
Delaigue Olivier committed
328
329
330
  

    ##Creation_of_new_candidates______________________________________________
331
332
    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
Delaigue Olivier's avatar
Delaigue Olivier committed
333
    ##Remplacement_of_non_optimised_values_____________________________________
334
335
336
337
338
339
340
341
    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
      return(x)})
    if (NParam > 1) {
      CandidatesParamR <- t(CandidatesParamR)
    } else {
      CandidatesParamR <- cbind(CandidatesParamR)
    }
Delaigue Olivier's avatar
Delaigue Olivier committed
342
343
344


    ##Loop_to_test_the_various_candidates_____________________________________
345
346
    iNewOptim <- 0
    for (iNew in 1:nrow(CandidatesParamR)) {
Delaigue Olivier's avatar
Delaigue Olivier committed
347
      ##Model_run
348
349
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Delaigue Olivier's avatar
Delaigue Olivier committed
350
      ##Calibration_criterion_computation
351
352
353
354
355
356
357
      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
358
    }
359
    NRuns <- NRuns + nrow(CandidatesParamR)
Delaigue Olivier's avatar
Delaigue Olivier committed
360
361
362


    ##When_a_progress_has_been_achieved_______________________________________
363
    if (iNewOptim != 0) {
Delaigue Olivier's avatar
Delaigue Olivier committed
364
      ##We_store_the_optimal_set
365
366
367
      OldParamOptimT <- NewParamOptimT
      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
      Compt <- Compt+1
Delaigue Olivier's avatar
Delaigue Olivier committed
368
      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
369
370
371
      if (Compt > 2*NParam) {
        Pace <- Pace * 2
        Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
372
373
      }
      ##We_update_PaceDiag
374
375
376
377
378
379
380
381
382
383
384
      VectPace <- NewParamOptimT-OldParamOptimT
      for (iC in 1:NParam) {
        if (OptimParam[iC]) { 
          if (VectPace[iC] != 0) {
            PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
          }
          if (VectPace[iC] == 0) {
            PaceDiag[iC] <- CLG*PaceDiag[iC]
          }
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
385
386
    } else {
    ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
387
388
      Pace <- Pace / 2
      Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
389
390
391
392
    }


    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
    if (ITER > 4*NParam) {
      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
410
        }
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
      }
      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
      ##Model_run
      Param <- CandidatesParamR[iNew, ]
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
      ##Calibration_criterion_computation
      OutputsCrit <- FUN_CRIT(InputsCrit, OutputsModel, verbose = FALSE)      
      if (OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim) {
        CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier
        iNewOptim <- iNew
      }
      ##When_a_progress_has_been_achieved
      if (iNewOptim != 0) {
        OldParamOptimT <- NewParamOptimT
        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
427
428
429
430
431

    }
    

    ##Results_archiving_______________________________________________________
432
433
434
435
436
    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
437
438
439
440



    } ##END_LOOP_ITER_________________________________________________________
441
    ITER <- ITER-1
Delaigue Olivier's avatar
Delaigue Olivier committed
442
443
444
    

    ##Case_when_the_starting_parameter_set_remains_the_best_solution__________
445
446
    if (CritOptim == CritStart & verbose) { 
      message("\t No progress achieved")
Delaigue Olivier's avatar
Delaigue Olivier committed
447
448
449
    }
    
    ##End_of_Steepest_Descent_Local_Search____________________________________
450
451
452
453
454
    ParamFinalR <- NewParamOptimR
    ParamFinalT <- NewParamOptimT
    CritFinal   <- CritOptim
    NIter       <- 1 + ITER
    if (verbose) { 
455
456
457
      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))
Delaigue Olivier's avatar
Delaigue Olivier committed
458
459
    }
    ##Results_archiving_______________________________________________________
460
461
462
463
464
465
  	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")
Delaigue Olivier's avatar
Delaigue Olivier committed
466

467
468
469
470
    BoolCrit_Actual <- InputsCrit$BoolCrit
	  BoolCrit_Actual[OutputsCrit$Ind_notcomputed] <- FALSE
    MatBoolCrit <- cbind(InputsCrit$BoolCrit, BoolCrit_Actual)
    colnames(MatBoolCrit) <- c("BoolCrit_Requested", "BoolCrit_Actual")
Delaigue Olivier's avatar
Delaigue Olivier committed
471
472
473


##_____Output______________________________________________________________________________
474
475
476
477
478
479
    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
480
481
482
483
484
485
486
487
488



}