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}