Newer
Older
Delaigue Olivier
committed
Calibration_Michel <- function(InputsModel,RunOptions,InputsCrit,CalibOptions,FUN_MOD,FUN_CRIT,FUN_TRANSFO=NULL, verbose = TRUE){
##_____Arguments_check_____________________________________________________________________
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if(inherits(RunOptions,"RunOptions")==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); }
if(inherits(InputsCrit,"InputsCrit")==FALSE){ stop("InputsCrit must be of class 'InputsCrit' \n"); return(NULL); }
if(inherits(CalibOptions,"CalibOptions")==FALSE){ stop("CalibOptions must be of class 'CalibOptions' \n"); return(NULL); }
if(inherits(CalibOptions,"HBAN")==FALSE){ stop("CalibOptions must be of class 'HBAN' if Calibration_Michel is used \n"); return(NULL); }
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
##_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; }
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); }
}
##_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); }
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
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){
##Managing_matrix_sizes
Nvalmax <- nrow(DistribParam);
NParam <- ncol(DistribParam);
##we_add_columns_to_MatDistrib_until_it_has_20_columns
DistribParam2 <- matrix(NA,nrow=Nvalmax,ncol=20);
DistribParam2[1:Nvalmax,1:NParam] <- DistribParam;
##we_check_the_number_of_values_to_test_for_each_param
NbDistrib <- rep(1,20);
for(iC in 1:20){ NbDistrib[iC] <- max( 1 , Nvalmax-sum(is.na(DistribParam2[,iC])) ); }
##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
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]){
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],
DistribParam2[iL16,16],DistribParam2[iL17,17],DistribParam2[iL18,18],DistribParam2[iL19,19],DistribParam2[iL20,20]);
} } } } }
} } } } }
} } } } }
} } } } }
MAT <- matrix(VECT,ncol=20,byrow=TRUE)[,1:NParam];
if(is.matrix(MAT)==FALSE){ MAT <- cbind(MAT); }
Output <- NULL;
Output$NewCandidates <- MAT;
return(Output);
}
##Creation_of_new_candidates_______________________________________________
if(PrefilteringType==1){ CandidatesParamR <- CalibOptions$StartParamList; }
if(PrefilteringType==2){ DistribParamR <- CalibOptions$StartParamDistrib; DistribParamR[,!CalibOptions$OptimParam] <- NA; CandidatesParamR <- ProposeCandidatesGrid(DistribParamR)$NewCandidates; }
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); });
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
##Loop_to_test_the_various_candidates______________________________________
iNewOptim <- 0;
Ncandidates <- nrow(CandidatesParamR);
Delaigue Olivier
committed
if(verbose & Ncandidates>1){
if(PrefilteringType==1){ cat(paste("\t List-Screening in progress (",sep="")); }
if(PrefilteringType==2){ cat(paste("\t Grid-Screening in progress (",sep="")); }
cat("0%");
}
for(iNew in 1:nrow(CandidatesParamR)){
Delaigue Olivier
committed
if(verbose & Ncandidates>1){
for(k in c(2,4,6,8)){ if(iNew==round(k/10*Ncandidates)){ cat(paste(" ",10*k,"%",sep="")); } }
}
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
iNewOptim <- iNew;
} }
##Storage_of_crit_info
if(is.null(CritName) | is.null(CritBestValue) | is.null(Multiplier)){
CritName <- OutputsCrit$CritName;
CritBestValue <- OutputsCrit$CritBestValue;
Multiplier <- OutputsCrit$Multiplier;
}
}
Delaigue Olivier
committed
if(verbose & Ncandidates>1){ cat(" 100%) \n"); }
##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);
Delaigue Olivier
committed
if(verbose){
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
if(Ncandidates> 1){ cat(paste("\t Screening completed (",NRuns," runs): \n",sep="")); }
if(Ncandidates==1){ cat(paste("\t Starting point for steepest-descent local search: \n",sep="")); }
cat(paste("\t Param = ",paste(formatC(ParamStartR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritStart*Multiplier,format="f",digits=4),"\n",sep=""));
}
##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==TRUE){ VECT <- c(VECT,PotentialCandidateT); }
}
}
}
Output <- NULL;
Output$NewCandidatesT <- matrix(VECT,ncol=NParam,byrow=TRUE);
return(Output);
}
##Initialisation_of_variables
Delaigue Olivier
committed
if(verbose){
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
243
244
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
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
cat("\t Steepest-descent local search in progress \n");
}
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)){
##Exit_loop_when_Pace_becomes_too_small___________________________________
if(Pace<0.01){ break; }
##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT,OldParamOptimT,RangesT,CalibOptions$OptimParam,Pace)$NewCandidatesT;
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##Remplacement_of_non_optimised_values_____________________________________
CandidatesParamR <- apply(CandidatesParamR,1,function(x){ x[!CalibOptions$OptimParam] <- CalibOptions$FixedParam[!CalibOptions$OptimParam]; return(x); });
if(NParam>1){ CandidatesParamR <- t(CandidatesParamR); } else { CandidatesParamR <- cbind(CandidatesParamR); }
##Loop_to_test_the_various_candidates_____________________________________
iNewOptim <- 0;
for(iNew in 1:nrow(CandidatesParamR)){
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
if(!is.na(OutputsCrit$CritValue)){ if(OutputsCrit$CritValue*OutputsCrit$Multiplier < CritOptim){
CritOptim <- OutputsCrit$CritValue*OutputsCrit$Multiplier;
iNewOptim <- iNew;
} }
}
NRuns <- NRuns+nrow(CandidatesParamR);
##When_a_progress_has_been_achieved_______________________________________
if(iNewOptim!=0){
##We_store_the_optimal_set
OldParamOptimT <- NewParamOptimT;
NewParamOptimT <- matrix(CandidatesParamT[iNewOptim,1:NParam],nrow=1);
Compt <- Compt+1;
##When_necessary_we_increase_the_pace ### if_successive_progress_occur_in_a_row
if(Compt>2*NParam){
Pace <- Pace*2;
Compt <- 0;
}
##We_update_PaceDiag
VectPace <- NewParamOptimT-OldParamOptimT;
for(iC in 1:NParam){ if(CalibOptions$OptimParam[iC]==TRUE){
if(VectPace[iC]!=0){ PaceDiag[iC] <- CLG*PaceDiag[iC]+(1-CLG)*VectPace[iC]; }
if(VectPace[iC]==0){ PaceDiag[iC] <- CLG*PaceDiag[iC]; }
} }
} else {
##When_no_progress_has_been_achieved_we_decrease_the_pace_________________
Pace <- Pace/2;
Compt <- 0;
}
##Test_of_an_additional_candidate_using_diagonal_progress_________________
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(CalibOptions$OptimParam[iC]==TRUE){
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]; }
} }
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
OutputsCrit <- FUN_CRIT(InputsCrit,OutputsModel);
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);
}
}
##Results_archiving_______________________________________________________
NewParamOptimR <- FUN_TRANSFO(NewParamOptimT,"TR");
HistParamR[ITER+1,] <- NewParamOptimR;
HistParamT[ITER+1,] <- NewParamOptimT;
HistCrit[ITER+1,] <- CritOptim;
Delaigue Olivier
committed
### 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="")); }
} ##END_LOOP_ITER_________________________________________________________
ITER <- ITER-1;
##Case_when_the_starting_parameter_set_remains_the_best_solution__________
Delaigue Olivier
committed
if(CritOptim==CritStart & verbose){
cat("\t No progress achieved \n");
}
##End_of_Steepest_Descent_Local_Search____________________________________
ParamFinalR <- NewParamOptimR;
ParamFinalT <- NewParamOptimT;
CritFinal <- CritOptim;
NIter <- 1+ITER;
Delaigue Olivier
committed
if(verbose){
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
cat(paste("\t Calibration completed (",NIter," iterations, ",NRuns," runs): \n",sep=""));
cat(paste("\t Param = ",paste(formatC(ParamFinalR,format="f",width=8,digits=3),collapse=" , "),"\n",sep=""));
cat(paste("\t Crit ",format(CritName,width=12,justify="left")," = ",formatC(CritFinal*Multiplier,format="f",digits=4),"\n",sep=""));
}
##Results_archiving_______________________________________________________
HistParamR <- cbind(HistParamR[1:NIter,]); colnames(HistParamR) <- paste("Param",1:NParam,sep="");
HistParamT <- cbind(HistParamT[1:NIter,]); colnames(HistParamT) <- paste("Param",1:NParam,sep="");
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(as.double(ParamFinalR),CritFinal*Multiplier,NIter,NRuns,HistParamR,HistCrit*Multiplier,MatBoolCrit,CritName,CritBestValue);
names(OutputsCalib) <- c("ParamFinalR","CritFinal","NIter","NRuns","HistParamR","HistCrit","MatBoolCrit","CritName","CritBestValue");
class(OutputsCalib) <- c("OutputsCalib","HBAN");
return(OutputsCalib);
}