An error occurred while loading the file. Please try again.
-
unknown authored
v1.0.5.22 remove the using of dev.new() in the plot.OutputsModel() when RStudio is not used because it causes some disorders when airGRteachin uses it in shiny (not plot in the interface when a script is called from Rscript.exe or R.exe)
2d033329
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); } ##_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); } 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);7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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_______________________________________________
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){
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)){
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;
}
}
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);
if(verbose){
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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
if(verbose){
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)){
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
##Exit_loop_when_Pace_becomes_too_small___________________________________
if(Pace<0.01){ break; }
##Creation_of_new_candidates______________________________________________
CandidatesParamT <- ProposeCandidatesLoc(NewParamOptimT,OldParamOptimT,RangesT,OptimParam,Pace)$NewCandidatesT;
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##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;
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(OptimParam[iC]){
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(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]; }
} }
CandidatesParamR <- FUN_TRANSFO(CandidatesParamT,"TR");
##Model_run
Param <- CandidatesParamR[iNew,];
OutputsModel <- FUN_MOD(InputsModel,RunOptions,Param);
##Calibration_criterion_computation
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347
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;
### 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__________
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;
if(verbose){
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);
}