From 035bfe12c6492c6a7b3a048762ac6cea703336ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roland=20Holl=C3=B3s?= Date: Sun, 17 May 2020 21:06:31 +0200 Subject: [PATCH] adding convergence tracker, change graph output, filtering behavior elems GLUE intervals --- RBBGCMuso/NAMESPACE | 1 + RBBGCMuso/R/calibration.R | 107 +++++++++++++++++++--------- RBBGCMuso/R/otherUsefullFunctions.R | 11 +++ RBBGCMuso/man/getConstMatrix.Rd | 17 +++++ RBBGCMuso/man/musoRand.Rd | 3 +- RBBGCMuso/man/optiMuso.Rd | 2 +- 6 files changed, 107 insertions(+), 34 deletions(-) create mode 100644 RBBGCMuso/man/getConstMatrix.Rd diff --git a/RBBGCMuso/NAMESPACE b/RBBGCMuso/NAMESPACE index 72b15fd..ede2399 100644 --- a/RBBGCMuso/NAMESPACE +++ b/RBBGCMuso/NAMESPACE @@ -9,6 +9,7 @@ export(copyMusoExampleTo) export(corrigMuso) export(createSoilFile) export(getAnnualOutputList) +export(getConstMatrix) export(getDailyOutputList) export(getMeteoData1BGC) export(getyearlycum) diff --git a/RBBGCMuso/R/calibration.R b/RBBGCMuso/R/calibration.R index fcf4135..105cc58 100644 --- a/RBBGCMuso/R/calibration.R +++ b/RBBGCMuso/R/calibration.R @@ -42,6 +42,7 @@ optiMuso <- function(measuredData, parameters = NULL, startDate = NULL, }, continious, modelVar = 3009, + naVal = NULL, postProcString = NULL) { mdata <- measuredData @@ -104,9 +105,14 @@ optiMuso <- function(measuredData, parameters = NULL, startDate = NULL, unlink randValues <- randVals[[2]] settings$calibrationPar <- randVals[[1]] + if(!is.null(naVal)){ + measuredData <- as.data.frame(measuredData) + measuredData[measuredData == naVal] <- NA + } 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 rmse <- numeric(iterations + 1) origModellOut <- calibMuso(settings=settings,silent=TRUE, skipSpinup = skipSpinup,postProcString=postProcString) @@ -114,59 +120,96 @@ optiMuso <- function(measuredData, parameters = NULL, startDate = NULL, write.csv(x=origModellOut, file=paste0(pretag,1,".csv")) modellOut[1] <- likelihood(measuredData,origModellOut[modIndex,colNumb]) + rmse[1] <- sqrt(mean((measuredData-origModellOut[modIndex,colNumb])^2)) print("Running the model with the random epc values...", quote = FALSE) if(!is.null(postProcString)){ colNumb <- length(settings$dailyVarCodes) + 1 } + + partialTemplate <- matrix(ncol=length(randVals[[1]])+2) + colN <- randVals[[1]] + colN[match(parameters[,2],randVals[[1]])] <- parameters[,1] + colnames(partialTemplate) <- c(colN, "rmse","likelihood") + partialTemplate[1:((length(randVals[[1]]))+2)] <- c(readValuesFromFile(settings$epc[2],randVals[[1]]),modellOut[1], rmse[1]) + write.csv(x=partialTemplate, file="preservedCalib.csv",row.names=FALSE) for(i in 2:(iterations+1)){ tmp <- tryCatch(calibMuso(settings = settings, parameters = randValues[(i-1),], silent= TRUE, - skipSpinup = skipSpinup, postProcString = postProcString)[modIndex,colNumb], error = function (e) NA) - # browser() + skipSpinup = skipSpinup, postProcString = postProcString)[modIndex,colNumb], error = function (e) NULL ) + if(is.null(tmp)){ + tmp <- rmse[i] <- modellOut[i] <- NA + } else { + modellOut[i]<- likelihood(measuredData,tmp) + rmse[i] <- sqrt(mean((measuredData-tmp)^2)) + } - modellOut[i]<- likelihood(measuredData,tmp) - rmse[i] <- sqrt(mean((measuredData-tmp)^2)) - write.csv(x=tmp, file=paste0(pretag,(i+1),".csv")) + partialTemplate[1:(length(randVals[[1]])+2)] <- c(randValues[(i-1),], rmse[i], modellOut[i]) + write.table(x=partialTemplate, file="preservedCalib.csv",append=TRUE,row.names=FALSE,sep=",",col.names=FALSE) + write.csv(x=tmp, file=paste0(pretag, (i+1),".csv")) setTxtProgressBar(progBar,i) } - paramLines <- parameters[,2] - paramLines <- order(paramLines) - randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])] - randInd <- order(randInd) + # paramLines <- parameters[,2] + # paramLines <- order(paramLines) + # randInd <- randVals[[1]][(randVals[[1]] %in% parameters[,2])] + # randInd <- order(randInd) + # + # + # + # epcStrip <- rbind(origEpc[order(parameters[,2])], + # randValues[,randVals[[1]] %in% parameters[,2]][,randInd]) + # + # + # preservedCalib <- cbind(epcStrip,rmse, + # modellOut) + # columNames <- c(parameterNames[paramLines],"rmse", "likelihood") + # colnames(preservedCalib) <- columNames + # write.csv(preservedCalib,"preservedCalib.csv") - - epcStrip <- rbind(origEpc[order(parameters[,2])], - randValues[,randVals[[1]] %in% parameters[,2]][,randInd]) - - - preservedCalib <- cbind(epcStrip,rmse, - modellOut) - columNames <- c(parameterNames[paramLines],"rmse", "likelihood") - colnames(preservedCalib) <- columNames - write.csv(preservedCalib,"preservedCalib.csv") + preservedCalib<- read.csv("preservedCalib.csv") p<-list() preservedCalib <- preservedCalib[-1,] dontInclude <-c((ncol(preservedCalib)-1),ncol(preservedCalib)) + + # 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) + # } + + unfilteredLikelihood <- preservedCalib$likelihood + preservedCalib <- preservedCalib[preservedCalib$likelihood>quantile(preservedCalib$likelihood,0.95),] + optRanges <-t(apply(preservedCalib,2,function(x) quantile(x,c(0.05,0.5,0.95)))) + +pdf("dotplot.pdf") + plot(Reduce(min, -log(unfilteredLikelihood), accumulate=TRUE),type="l", ylab="-log(likelihood)",xlab="iterations") 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) + plot(preservedCalib[,i],preservedCalib[,"likelihood"],pch=19,cex=.1, ylab="likelihood", + main = colnames(preservedCalib)[i], xlab="") + abline(v=optRanges[i,1],col="blue") + abline(v=optRanges[i,2],col="green") + abline(v=optRanges[i,3],col="red") + } +dev.off() - ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))), dpi = 1000) - maxLikelihoodPlace <- which(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,]) + # ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))), dpi = 1000) + # maxLikelihoodPlace <- which(preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"],na.rm = TRUE)) + # resPlot <- plotMusoWithData(mdata = measuredData, 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,]) + write.csv(optRanges,"optRanges.csv") + return(head(optRanges,n=-2)) } + + diff --git a/RBBGCMuso/R/otherUsefullFunctions.R b/RBBGCMuso/R/otherUsefullFunctions.R index 2385911..f70bdff 100644 --- a/RBBGCMuso/R/otherUsefullFunctions.R +++ b/RBBGCMuso/R/otherUsefullFunctions.R @@ -172,3 +172,14 @@ while(n<=N){ } return(randomNorm) } + +#' getConstMatrix +#' +#' getConstMatrix is a function whith wich you can get the default constrain matrix for your choosen type and version. +#' @param filetype It can be "epc" or "soil". +#' @param version The version of the MuSo environment +#' @export + +getConstMatrix <- function (filetype="epc", version = as.character(getOption("RMuso_version"))) { + getOption("RMuso_constMatrix")[[filetype]][[version]] +} diff --git a/RBBGCMuso/man/getConstMatrix.Rd b/RBBGCMuso/man/getConstMatrix.Rd new file mode 100644 index 0000000..234fb65 --- /dev/null +++ b/RBBGCMuso/man/getConstMatrix.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/otherUsefullFunctions.R +\name{getConstMatrix} +\alias{getConstMatrix} +\title{getConstMatrix} +\usage{ +getConstMatrix(filetype = "epc", + version = as.character(getOption("RMuso_version"))) +} +\arguments{ +\item{filetype}{It can be "epc" or "soil".} + +\item{version}{The version of the MuSo environment} +} +\description{ +getConstMatrix is a function whith wich you can get the default constrain matrix for your choosen type and version. +} diff --git a/RBBGCMuso/man/musoRand.Rd b/RBBGCMuso/man/musoRand.Rd index 208748b..5fff6d5 100644 --- a/RBBGCMuso/man/musoRand.Rd +++ b/RBBGCMuso/man/musoRand.Rd @@ -5,7 +5,8 @@ \title{musoRand} \usage{ musoRand(parameters, iterations = 3000, fileType = "epc", - constrains = NULL) + version = as.character(getOption("RMuso_version")), + constrains = getConstMatrix(fileType = fileType, version = version)) } \arguments{ \item{parameters}{This is a dataframe (heterogeneous data-matrix), where the first column is the name of the parameter, the second is a numeric vector of the rownumbers of the given variable in the input EPC file, and the last two columns describe the minimum and the maximum of the parameter (i.e. the parameter ranges), defining the interval for the randomization.} diff --git a/RBBGCMuso/man/optiMuso.Rd b/RBBGCMuso/man/optiMuso.Rd index 4e4a64d..191b22e 100644 --- a/RBBGCMuso/man/optiMuso.Rd +++ b/RBBGCMuso/man/optiMuso.Rd @@ -11,7 +11,7 @@ optiMuso(measuredData, parameters = NULL, startDate = NULL, iterations = 30, skipSpinup = TRUE, constrains = NULL, plotName = "calib.jpg", likelihood = function(x, y) { exp(-sqrt(mean((x - y)^2))) }, continious, modelVar = 3009, - postProcString = NULL) + naVal = NULL, postProcString = NULL) } \arguments{ \item{parameters}{b}