adding convergence tracker, change graph output, filtering behavior elems GLUE intervals

RBBGCMuso/NAMESPACE

@@ -9,6 +9,7 @@ export(copyMusoExampleTo)
 export(corrigMuso)
 export(createSoilFile)
 export(getAnnualOutputList)
+export(getConstMatrix)
 export(getDailyOutputList)
 export(getMeteoData1BGC)
 export(getyearlycum)

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)
+                                  skipSpinup = skipSpinup, postProcString = postProcString)[modIndex,colNumb], error = function (e) NULL        )
-        # browser()
+        if(is.null(tmp)){
-
+           tmp <- rmse[i] <- modellOut[i] <- NA
+        } else {
             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")
 
-
+    preservedCalib<- read.csv("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")
     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")) +
+        plot(preservedCalib[,i],preservedCalib[,"likelihood"],pch=19,cex=.1, ylab="likelihood",
-            geom_point(shape='.',size=1,alpha=0.8)
+             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)
+    # ggsave(plotName,grid.arrange(grobs = p, ncol = floor(sqrt(ncol(preservedCalib)-1))), dpi = 1000)
-    maxLikelihoodPlace  <- which(preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"],na.rm = TRUE))
+    # maxLikelihoodPlace  <- which(preservedCalib[,"likelihood"]==max(preservedCalib[,"likelihood"],na.rm = TRUE))
-    resPlot <- plotMusoWithData(mdata = mdata, startDate = startDate, endDate = endDate,
+    # resPlot <- plotMusoWithData(mdata = measuredData, startDate = startDate, endDate = endDate,
-                                dataVar = dataVar, modelVar = modelVar, settings = settings, continious = continious) +
+    #                             dataVar = dataVar, modelVar = modelVar, settings = settings, continious = continious) +
-        plotMuso(settings = settings, parameters = randValues[maxLikelihoodPlace,],
+    #     plotMuso(settings = settings, parameters = randValues[maxLikelihoodPlace,],
-                 postProcString = postProcString, skipSpinup = FALSE, variable = colNumb, layerPlot = TRUE, colour = "green")
+    #              postProcString = postProcString, skipSpinup = FALSE, variable = colNumb, layerPlot = TRUE, colour = "green")
-
+    #
-    print(resPlot)
+    # print(resPlot)
-    tempEpc <- paste0(tools::file_path_sans_ext(basename(settings$epcInput[2])),"-tmp.",tools::file_ext(settings$epcInput[2]))
+    # tempEpc <- paste0(tools::file_path_sans_ext(basename(settings$epcInput[2])),"-tmp.",tools::file_ext(settings$epcInput[2]))
-    file.rename(tempEpc, "optimizedEpc.epc")
+    # file.rename(tempEpc, "optimizedEpc.epc")
-    return(preservedCalib[maxLikelihoodPlace,])
+    # return(preservedCalib[maxLikelihoodPlace,])
+    write.csv(optRanges,"optRanges.csv")
+    return(head(optRanges,n=-2))
 }
 
 
+
+

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]]
+}

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.
+}

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.}

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}