Calibration_Michel.R 18.6 KB
Newer Older
1
2
3
4
5
6
7
8
Calibration_Michel <- function(InputsModel, 
                               RunOptions, 
                               InputsCrit, 
                               CalibOptions,
                               FUN_MOD, 
                               FUN_CRIT,           # deprecated
                               FUN_TRANSFO = NULL, 
                               verbose = TRUE) {
9
10
  
  
11
  FUN_MOD  <- match.fun(FUN_MOD)
12
13
14
  if (!missing(FUN_CRIT)) {
    FUN_CRIT <- match.fun(FUN_CRIT)
  }
15
16
17
18
19
  if (!is.null(FUN_TRANSFO)) {
    FUN_TRANSFO <- match.fun(FUN_TRANSFO)
  }
  
  
20
21
  ##_____Arguments_check_____________________________________________________________________
  if (!inherits(InputsModel, "InputsModel")) {
22
    stop("'InputsModel' must be of class 'InputsModel'")
23
24
  }  
  if (!inherits(RunOptions, "RunOptions")) {
25
    stop("'RunOptions' must be of class 'RunOptions'")
26
27
  }  
  if (!inherits(InputsCrit, "InputsCrit")) {
28
    stop("'InputsCrit' must be of class 'InputsCrit'")
29
30
  }
  if (inherits(InputsCrit, "Multi")) {
31
    stop("'InputsCrit' must be of class 'Single' or 'Compo'")
32
  }
33
  if (inherits(InputsCrit, "Single")) {
34
    listVarObs <- InputsCrit$VarObs
35
36
  }
  if (inherits(InputsCrit, "Compo")) {
37
    listVarObs <- sapply(InputsCrit, FUN = "[[", "VarObs")
38
39
  }
  if ("SCA" %in% listVarObs & !"Gratio" %in% RunOptions$Outputs_Cal) {
40
41
42
    warning("Missing 'Gratio' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SCA")
    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "Gratio")
  }
43
  if ("SWE" %in% listVarObs & !"SnowPack" %in% RunOptions$Outputs_Cal) {
44
45
46
    warning("Missing 'SnowPack' is automatically added to 'Output_Cal' in 'RunOptions' as it is necessary in the objective function for comparison with SWE")
    RunOptions$Outputs_Cal <- c(RunOptions$Outputs_Cal, "SnowPack")
  }
47
  if (!inherits(CalibOptions, "CalibOptions")) {
48
    stop("'CalibOptions' must be of class 'CalibOptions'")
49
50
  }  
  if (!inherits(CalibOptions, "HBAN")) {
51
    stop("'CalibOptions' must be of class 'HBAN' if 'Calibration_Michel' is used")
52
53
  }
  if (!missing(FUN_CRIT)) {
54
    warning("argument 'FUN_CRIT' is deprecated. The error criterion function is now automatically get from the 'InputsCrit' object")
55
  }  
56
  
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
  
  ##_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    )) {
79
      if (inherits(CalibOptions, "hysteresis")) {
80
81
82
83
        FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
      } else {
        FUN_TRANSFO <- TransfoParam_CemaNeige
      }
84
85
86
87
    }
    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
88
      }
89
90
      if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
        FUN1 <- TransfoParam_GR5J
91
      }
92
      if (identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
93
        FUN1 <- TransfoParam_GR6J
94
      }
95
      if (inherits(CalibOptions, "hysteresis")) {
96
97
        FUN2 <- TransfoParam_CemaNeigeHyst
      } else {
98
        FUN2 <- TransfoParam_CemaNeige
99
      }
100
101
102
      if (inherits(CalibOptions, "hysteresis")) {
        FUN_TRANSFO <- function(ParamIn, Direction) {
          Bool <- is.matrix(ParamIn)
103
          if (!Bool) {
104
105
106
107
108
109
            ParamIn <- rbind(ParamIn)
          }
          ParamOut <- NA * ParamIn
          NParam   <- ncol(ParamIn)
          ParamOut[,          1:(NParam-4)] <- FUN1(ParamIn[,          1:(NParam-4)], Direction)
          ParamOut[, (NParam-3):NParam    ] <- FUN2(ParamIn[, (NParam-3):NParam    ], Direction)
110
          if (!Bool) {
111
112
113
            ParamOut <- ParamOut[1, ]
          }
          return(ParamOut)
114
        }
115
116
117
      } else {
        FUN_TRANSFO <- function(ParamIn, Direction) {
          Bool <- is.matrix(ParamIn)
118
          if (!Bool) {
119
120
121
122
123
124
            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)
125
          if (!Bool) {
126
127
128
            ParamOut <- ParamOut[1, ]
          }
          return(ParamOut)
129
130
        }
      }
131
      
132
    }
133
    if (is.null(FUN_TRANSFO)) {
134
      stop("'FUN_TRANSFO' was not found (in 'Calibration' function)")
135
    }
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
  }
  
  ##_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) {
156
    stop("Calibration_Michel can handle a maximum of 20 parameters")
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
  }
  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) {
206
    if (PrefilteringType == 1) {
207
      message("List-Screening in progress (", appendLF = FALSE)
208
209
    }
    if (PrefilteringType == 2) {
210
      message("Grid-Screening in progress (", appendLF = FALSE)
211
    }
212
213
214
    message("0%", appendLF = FALSE)
  }
  for (iNew in 1:nrow(CandidatesParamR)) {
215
    if (verbose & Ncandidates > 1) {
216
217
218
      for (k in c(2, 4, 6, 8)) {
        if (iNew == round(k / 10 * Ncandidates)) {
          message(" ", 10 * k, "%", appendLF = FALSE)
219
        }
220
221
222
223
      } 
    }
    ##Model_run
    Param <- CandidatesParamR[iNew, ]
224
225
    OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
    
226
    ##Calibration_criterion_computation
227
    OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
228
229
230
231
    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
232
233
      }
    }
234
235
236
237
238
    ##Storage_of_crit_info
    if (is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)) {
      CritName      <- OutputsCrit$CritName
      CritBestValue <- OutputsCrit$CritBestValue
      Multiplier    <- OutputsCrit$Multiplier
239
    }
240
241
242
243
244
245
246
247
248
249
250
251
252
  }
  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
253
  NRuns       <- NRuns + nrow(CandidatesParamR)
254
255
256
  if (verbose) {
    if (Ncandidates > 1) {
      message(sprintf("\t Screening completed (%s runs)", NRuns))
257
    }
258
259
    if (Ncandidates == 1) {
      message("\t Starting point for steepest-descent local search:")
Delaigue Olivier's avatar
Delaigue Olivier committed
260
    }
261
    message("\t     Param = ", paste(sprintf("%8.3f", ParamStartR), collapse = ", "))
262
    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritStart * Multiplier))
263
264
265
266
267
268
269
270
271
272
273
274
275
  }
  ##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
276
  ProposeCandidatesLoc <- function(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace) {
277
278
    ##Format_checking
    if (nrow(NewParamOptimT) != 1 | nrow(OldParamOptimT) != 1) {
279
      stop("each input set must be a matrix of one single line")
280
281
    }
    if (ncol(NewParamOptimT)!=ncol(OldParamOptimT) | ncol(NewParamOptimT) != length(OptimParam)) {
282
      stop("each input set must have the same number of values")
283
284
285
286
287
288
    }
    ##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
289
      if (OptimParam[I]) {
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
        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
317
318
319
320
          }
        }
      }
    }
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
    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
347
    ##Exit_loop_when_Pace_becomes_too_small___________________________________
348
349
350
    if (Pace < 0.01) {
      break
    }
351
352
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
353
    ##Creation_of_new_candidates______________________________________________
354
355
    CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT, OldParamOptimT, RangesT, OptimParam, Pace)$NewCandidatesT
    CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
Delaigue Olivier's avatar
Delaigue Olivier committed
356
    ##Remplacement_of_non_optimised_values_____________________________________
357
358
    CandidatesParamR <- apply(CandidatesParamR, 1, function(x) {
      x[!OptimParam] <- CalibOptions$FixedParam[!OptimParam]
359
360
      return(x)
    })
361
362
363
364
365
    if (NParam > 1) {
      CandidatesParamR <- t(CandidatesParamR)
    } else {
      CandidatesParamR <- cbind(CandidatesParamR)
    }
366
367
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
368
    ##Loop_to_test_the_various_candidates_____________________________________
369
370
    iNewOptim <- 0
    for (iNew in 1:nrow(CandidatesParamR)) {
Delaigue Olivier's avatar
Delaigue Olivier committed
371
      ##Model_run
372
      Param <- CandidatesParamR[iNew, ]
373
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
Delaigue Olivier's avatar
Delaigue Olivier committed
374
      ##Calibration_criterion_computation
375
      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)      
376
      if (!is.na(OutputsCrit$CritValue)) {
377
378
        if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
          CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
379
380
381
          iNewOptim <- iNew
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
382
    }
383
    NRuns <- NRuns + nrow(CandidatesParamR)
384
385
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
386
    ##When_a_progress_has_been_achieved_______________________________________
387
    if (iNewOptim != 0) {
Delaigue Olivier's avatar
Delaigue Olivier committed
388
      ##We_store_the_optimal_set
389
390
      OldParamOptimT <- NewParamOptimT
      NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
391
      Compt <- Compt + 1
Delaigue Olivier's avatar
Delaigue Olivier committed
392
      ##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
393
      if (Compt > 2 * NParam) {
394
395
        Pace <- Pace * 2
        Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
396
397
      }
      ##We_update_PaceDiag
398
399
400
      VectPace <- NewParamOptimT-OldParamOptimT
      for (iC in 1:NParam) {
        if (OptimParam[iC]) { 
401
          PaceDiag[iC] <- CLG * PaceDiag[iC] + (1-CLG) * VectPace[iC]
402
403
        }
      }
Delaigue Olivier's avatar
Delaigue Olivier committed
404
    } else {
405
      ##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
406
407
      Pace <- Pace / 2
      Compt <- 0
Delaigue Olivier's avatar
Delaigue Olivier committed
408
    }
409
410
    
    
Delaigue Olivier's avatar
Delaigue Olivier committed
411
    ##Test_of_an_additional_candidate_using_diagonal_progress_________________
412
    if (ITER > 4 * NParam) {
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
      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
429
        }
430
431
432
433
      }
      CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
      ##Model_run
      Param <- CandidatesParamR[iNew, ]
434
      OutputsModel <- FUN_MOD(InputsModel, RunOptions, Param)
435
      ##Calibration_criterion_computation
436
      OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
437
438
      if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
        CritOptim <- OutputsCrit$CritValue * OutputsCrit$Multiplier
439
440
441
442
443
444
445
        iNewOptim <- iNew
      }
      ##When_a_progress_has_been_achieved
      if (iNewOptim != 0) {
        OldParamOptimT <- NewParamOptimT
        NewParamOptimT <- matrix(CandidatesParamT[iNewOptim, 1:NParam], nrow = 1)
      }
446
      
Delaigue Olivier's avatar
Delaigue Olivier committed
447
448
    }
    
449
    
Delaigue Olivier's avatar
Delaigue Olivier committed
450
    ##Results_archiving_______________________________________________________
451
452
453
454
455
    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
456
457
    
    
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
    
  } ##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))
475
476
    message("\t     Param = ", paste(sprintf("%8.3f", ParamFinalR), collapse = ", "))
    message(sprintf("\t     Crit. %-12s = %.4f", CritName, CritFinal * Multiplier))
477
478
479
480
481
    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 = ","))
482
      msgFormSep <- rep(c(",", ",", ",\n\t\t    "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
483
      msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
484
      msgForm <- gsub("\\,\\\n\\\t\\\t    $|\\,$", "", msgForm)
485
      message("\tFormula: sum(", msgForm, ")")
486
    }
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
  }
  ##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
513
514
}