hollorol/RBBGCMuso
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

plot optimized muso

Browse Source
This commit is contained in:
Roland Hollós 2019-03-11 13:06:37 +01:00
parent 760bbaef9a
commit 4b1c5fe5a3
3 changed files with 72 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
RBBGCMuso/R/calibration.R
View File

@ -25,7 +25,7 @@
#' @importFrom ggplot2 ggplot aes_string geom_point ggsave #' @importFrom ggplot2 ggplot aes_string geom_point ggsave
#' @importFrom magrittr '%>%' #' @importFrom magrittr '%>%'
#' @importFrom gridExtra grid.arrange #' @importFrom gridExtra grid.arrange
#' @export #' @export
optiMuso <- function(measuredData, parameters = NULL, startDate, optiMuso <- function(measuredData, parameters = NULL, startDate,
endDate, formatString = "%Y-%m-%d", endDate, formatString = "%Y-%m-%d",
leapYearHandling = TRUE, leapYearHandling = TRUE,
@ -44,8 +44,9 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
modelVar = 3009, modelVar = 3009,
postProcString = NULL) postProcString = NULL)
{ {
mdata <- measuredData
dataCol <- grep(dataVar, colnames(measuredData)) dataCol <- grep(dataVar, colnames(measuredData))
if(is.null(parameters)){ if(is.null(parameters)){
parameters <- tryCatch(read.csv("parameters.csv", stringsAsFactor=FALSE), error = function (e) { parameters <- tryCatch(read.csv("parameters.csv", stringsAsFactor=FALSE), error = function (e) {
stop("You need to specify a path for the parameters.csv, or a matrix.") stop("You need to specify a path for the parameters.csv, or a matrix.")
@ -56,26 +57,26 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
stop("Cannot find neither parameters file neither the parameters matrix") stop("Cannot find neither parameters file neither the parameters matrix")
}) })
}} }}
outLoc <- normalizePath(outLoc) outLoc <- normalizePath(outLoc)
outLocPlain <- basename(outLoc) outLocPlain <- basename(outLoc)
currDir <- getwd() currDir <- getwd()
if(!dir.exists(outLoc)){ if(!dir.exists(outLoc)){
dir.create(outLoc) dir.create(outLoc)
warning(paste(outLoc," is not exists, so it was created")) warning(paste(outLoc," is not exists, so it was created"))
} }
outLoc <- normalizePath(outLoc) outLoc <- normalizePath(outLoc)
if(is.null(settings)){ if(is.null(settings)){
settings <- setupMuso() settings <- setupMuso()
} }
parameterNames <- parameters[,1] parameterNames <- parameters[,1]
pretag <- file.path(outLoc,preTag) pretag <- file.path(outLoc,preTag)
npar <- length(settings$calibrationPar) npar <- length(settings$calibrationPar)
##reading the original epc file at the specified ##reading the original epc file at the specified
## row numbers ## row numbers
print("optiMuso is randomizing the epc parameters now...",quote = FALSE) print("optiMuso is randomizing the epc parameters now...",quote = FALSE)
@ -85,17 +86,17 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
} else { } else {
randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = iterations) randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = iterations)
} }
origEpc <- readValuesFromFile(settings$epc[2],parameters[,2]) origEpc <- readValuesFromFile(settings$epc[2],parameters[,2])
## Prepare the preservedCalib matrix for the faster ## Prepare the preservedCalib matrix for the faster
## run. ## run.
pretag <- file.path(outLoc,preTag) pretag <- file.path(outLoc,preTag)
## Creating function for generating separate ## Creating function for generating separate
## csv files for each run ## csv files for each run
progBar <- txtProgressBar(1,iterations,style=3) progBar <- txtProgressBar(1,iterations,style=3)
colNumb <- which(settings$dailyVarCodes == modelVar) colNumb <- which(settings$dailyVarCodes == modelVar)
settings$iniInput[2] %>% settings$iniInput[2] %>%
@ -104,12 +105,13 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
randValues <- randVals[[2]] randValues <- randVals[[2]]
settings$calibrationPar <- randVals[[1]] settings$calibrationPar <- randVals[[1]]
list2env(alignData(measuredData,dataCol = dataCol,modellSettings = settings,startDate = startDate,endDate = endDate,leapYear = leapYearHandling, continious = continious),envir=environment()) list2env(alignData(measuredData,dataCol = dataCol,modellSettings = settings,startDate = startDate,endDate = endDate,leapYear = leapYearHandling, continious = continious),envir=environment())
## modIndex and measuredData are created.
modellOut <- numeric(iterations + 1) # single variable solution modellOut <- numeric(iterations + 1) # single variable solution
rmse <- numeric(iterations + 1) rmse <- numeric(iterations + 1)
origModellOut <- calibMuso(settings=settings,silent=TRUE, skipSpinup = skipSpinup) origModellOut <- calibMuso(settings=settings,silent=TRUE, skipSpinup = skipSpinup)
write.csv(x=origModellOut, file=paste0(pretag,1,".csv")) write.csv(x=origModellOut, file=paste0(pretag,1,".csv"))
modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb]) modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb])
print("Running the model with the random epc values...", quote = FALSE) print("Running the model with the random epc values...", quote = FALSE)
@ -117,13 +119,13 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
if(!is.null(postProcString)){ if(!is.null(postProcString)){
colNumb <- length(settings$dailyVarCodes) + 1 colNumb <- length(settings$dailyVarCodes) + 1
} }
for(i in 2:(iterations+1)){ for(i in 2:(iterations+1)){
tmp <- tryCatch(calibMuso(settings = settings, tmp <- tryCatch(calibMuso(settings = settings,
parameters = randValues[(i-1),], parameters = randValues[(i-1),],
silent= TRUE, silent= TRUE,
skipSpinup = skipSpinup, postProcString = postProcString)[modIndex,colNumb], error = function (e) NA) skipSpinup = skipSpinup, postProcString = postProcString)[modIndex,colNumb], error = function (e) NA)
modellOut[i]<- likelihood(measuredData,tmp) modellOut[i]<- likelihood(measuredData,tmp)
rmse[i] <- sqrt(mean((measuredData-tmp)^2)) rmse[i] <- sqrt(mean((measuredData-tmp)^2))
write.csv(x=tmp, file=paste0(pretag,(i+1),".csv")) write.csv(x=tmp, file=paste0(pretag,(i+1),".csv"))
@ -133,13 +135,13 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
paramLines <- order(paramLines) paramLines <- order(paramLines)
randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])] randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])]
randInd <- order(randInd) randInd <- order(randInd)
epcStrip <- rbind(origEpc[order(parameters[,2])], epcStrip <- rbind(origEpc[order(parameters[,2])],
randValues[,randVals[[1]] %in% parameters[,2]][,randInd]) randValues[,randVals[[1]] %in% parameters[,2]][,randInd])
preservedCalib <- cbind(epcStrip,rmse, preservedCalib <- cbind(epcStrip,rmse,
modellOut) modellOut)
columNames <- c(parameterNames[paramLines],"rmse", "likelihood") columNames <- c(parameterNames[paramLines],"rmse", "likelihood")
@ -149,12 +151,21 @@ optiMuso <- function(measuredData, parameters = NULL, startDate,
preservedCalib <- preservedCalib[-1,] preservedCalib <- preservedCalib[-1,]
dontInclude <-c((ncol(preservedCalib)-1),ncol(preservedCalib)) dontInclude <-c((ncol(preservedCalib)-1),ncol(preservedCalib))
for(i in seq_along(colnames(preservedCalib)[-dontInclude])){ for(i in seq_along(colnames(preservedCalib)[-dontInclude])){
p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood"))+geom_point(shape='.',size=1,alpha=0.8) p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood")) +
geom_point(shape='.',size=1,alpha=0.8)
} }
ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))),dpi = 3000) ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))),dpi = 300)
maxLikelihoodPlace <- which(preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"],na.rm = TRUE))
return(preservedCalib[preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"],na.rm = TRUE),]) resPlot <- plotMusoWithData(mdata = mdata, startDate = startDate, endDate = endDate,
dataVar = dataVar, modelVar = modelVar, settings = settings, continious = continious) +
plotMuso(settings = settings, parameters = randValues[maxLikelihoodPlace,],
postProcString = postProcString, skipSpinup = FALSE, variable = colNumb, layerPlot = TRUE, colour = "green")
print(resPlot)
tempEpc <- paste0(tools::file_path_sans_ext(basename(settings$epcInput[2])),"-tmp.",tools::file_ext(settings$epcInput[2]))
file.rename(tempEpc, "optimizedEpc.epc")
return(preservedCalib[maxLikelihoodPlace,])
} }

17
RBBGCMuso/R/musoTime.R
View File

@ -59,17 +59,22 @@ musoDate <- function(startYear, endYears = NULL, numYears, combined = TRUE, leap
#' This function align the data to the model and the model to the data #' This function align the data to the model and the model to the data
#' @importFrom lubridate leap_year #' @importFrom lubridate leap_year
#' @keywords internal #' @keywords internal
alignData <- function(mdata, dataCol, modellSettings = NULL, startDate=NULL, endDate=NULL, formatString = "%Y-%m-%d", leapYear = TRUE, continious = TRUE){ alignData <- function(mdata, dataCol, modellSettings = NULL, startDate=NULL, endDate=NULL, formatString = "%Y-%m-%d", leapYear = TRUE, continious = FALSE){
if(continious){
if((is.null(startDate) | is.null(endDate))){
stop("If your date is continuous, you have to provide both startDate and endDate. ")
}
startDate <- as.Date(startDate, format = formatString)
endDate <- as.Date(endDate, format = formatString)
}
startDate <- as.Date(startDate, format = formatString)
endDate <- as.Date(endDate, format = formatString)
mdata <- as.data.frame(mdata)
if(is.null(modellSettings)){ if(is.null(modellSettings)){
modellSettings <- setupMuso() modellSettings <- setupMuso()
} }
mdata <- as.data.frame(mdata)
if(continious){ if(continious){
dates <- seq(startDate, to = endDate, by= "day") dates <- seq(startDate, to = endDate, by= "day")
} else{ } else{

32
RBBGCMuso/R/plotMuso.R
View File

@ -251,8 +251,12 @@ plotMuso <- function(settings = NULL, variable = 1,
#' @importFrom ggplot2 ggplot geom_line geom_point aes aes_string labs theme element_blank #' @importFrom ggplot2 ggplot geom_line geom_point aes aes_string labs theme element_blank
#' @export #' @export
plotMusoWithData <- function(mdata, plotName=NULL, plotMusoWithData <- function(mdata, plotName=NULL,
startDate, endDate, startDate = NULL, endDate = NULL,
colour=c("black","blue"),dataVar, modelVar, settings = setupMuso(), silent = TRUE, continious = TRUE){ colour=c("black","blue"), dataVar, modelVar, settings = setupMuso(), silent = TRUE, continious = FALSE){
if(continious & (is.null(startDate) | is.null(endDate))){
stop("If your date is continuous, you have to provide both startDate and endDate. ")
}
dataCol<- grep(paste0("^",dataVar,"$"), colnames(mdata)) dataCol<- grep(paste0("^",dataVar,"$"), colnames(mdata))
selVar <- grep(modelVar,(settings$dailyVarCodes))+4 selVar <- grep(modelVar,(settings$dailyVarCodes))+4
@ -261,24 +265,36 @@ plotMusoWithData <- function(mdata, plotName=NULL,
modellSettings = settings, modellSettings = settings,
startDate = startDate, startDate = startDate,
endDate = endDate, leapYear = FALSE, continious = continious),envir=environment()) endDate = endDate, leapYear = FALSE, continious = continious),envir=environment())
mesData <- numeric(settings$numYears*365)
k <- 1
for(i in seq(mesData)){
if(i %in% modIndex){
mesData[i] <- measuredData[k]
k <- k + 1
} else {
mesData[i] <- NA
}
}
rm(k)
# modIndex and measuredData are created.
## measuredData is created ## measuredData is created
baseData <- calibMuso(settings = settings, silent = silent, prettyOut = TRUE)[modIndex,] ## baseData <- calibMuso(settings = settings, silent = silent, prettyOut = TRUE)[modIndex,]
baseData <- calibMuso(settings = settings, silent = silent, prettyOut = TRUE)
baseData[,1] <- as.Date(baseData[,1],format = "%d.%m.%Y") baseData[,1] <- as.Date(baseData[,1],format = "%d.%m.%Y")
selVarName <- colnames(baseData)[selVar] selVarName <- colnames(baseData)[selVar]
if(!all.equal(colnames(baseData),unique(colnames(baseData)))){ if(!all.equal(colnames(baseData),unique(colnames(baseData)))){
notUnique <- setdiff((unlist(settings$dailyVarCodes)),unique(unlist(settings$dailyVarCodes))) notUnique <- setdiff((unlist(settings$dailyVarCodes)),unique(unlist(settings$dailyVarCodes)))
stop(paste0("Error: daily output variable list in the ini file must contain unique numbers. Check your ini files! Not unique codes: ",notUnique)) stop(paste0("Error: daily output variable list in the ini file must contain unique numbers. Check your ini files! Not unique codes: ",notUnique))
} }
mesData<-cbind.data.frame(baseData[,1],mesData)
colnames(mesData) <- c("date", "measured")
p <- baseData %>% p <- baseData %>%
ggplot(aes_string("date",selVarName)) + ggplot(aes_string("date",selVarName)) +
geom_line(colour=colour[1]) + geom_line(colour=colour[1]) +
geom_point(colour=colour[2], aes(date,measuredData)) + geom_point(data = mesData, colour=colour[2], aes(date,measured)) +
labs(y = paste0(selVarName,"_measured"))+ labs(y = paste0(selVarName,"_measured"))+
theme(axis.title.x = element_blank()) theme(axis.title.x = element_blank())
if(!is.null(plotName)){ if(!is.null(plotName)){
ggsave(plotName,p) ggsave(plotName,p)
return(p) return(p)
} else { } else {