Adding sensitivity analysis

RBBGCMuso/R/calibMuso.R

@@ -19,7 +19,6 @@
 #' @import utils
 #' @export
 
-
 calibMuso <- function(settings,parameters=NULL, timee="d", debugging=FALSE, logfilename=NULL, keepEpc=FALSE, export=FALSE, silent=FALSE, aggressive=FALSE, leapYear=FALSE,keepBinary=FALSE, binayPlace="./", fileToChange="epc"){
 
 

RBBGCMuso/R/musoMonte.R

@@ -3,147 +3,150 @@
 #' This funcion is fundamental for the BiomBGC-MuSo modell related functions like spinupMuso, normalMuso, rungetMuso, because it sets the modells environment.
 #' inputDir = "./",
 #' @author Roland Hollos
-#' @param settings Do you want to run multiple modell paralelly, if yes, set this variable to TRUE
+#' @param settings A list of montecarlos environmental variables. It is generated by the setupMuso() function. In default the settings parameter is generated automatically. 
-#' @param parameters stores the place of the modell-executable file. In normal usage, you don't have to be set this, because a RBBgcmuso package contains allways the latest modell executable. In spite of this, if you would like to use this package for modell development or just want to use different models (for example for comparison), you will find it useful
+#' @param parameters This is a dataframe (heterogen data-matrix), which first column is the name of the parameters, the second is a numeric vector of the rownumbers of the given variable in the epc-fie, the last two column consist the endpont of the parameter-ranges, where the parameters will be randomized.
 #' @param calibrationPar You may want to change some parameters in your epc file, before you run the modell. You have to select the appropirate modell parameters. You can refence to these with the number of the line in the epc file where the variables are. It indexes from one. You should use a vector for this, like: c(1,5,8)
-#' @param outLoc Where should the modell puts its outputs. You should give a location for it via this variable, for example: outputLoc="/place/of/the/outputs/"
+#' @param inputDir The location of the input directory, this directory must content a viable pack of all inputfiles and the executable file. 
-#' @param iterations Usually it is the root directory, where you put the iniFiles for the modell
+#' @param iterations number of the monteCarlo run.
-#' @param preTag Via metInput parameter, you can tell the modell where are the meteorological files. As default it reads this from the iniFiles.
+#' @param preTag It will be the name of the output files. For example preTag-1.csv, pretag-2csv...
-#' @param inputName Via CO2 parameter, you can tell the modell where are the CO2 data files. As default it reads this from the iniFiles.
 #' @param outputType Via planting parameter, you can tell the modell where are the data files, which contains the planting informations. As default it reads this from the iniFiles.
 #' @param fun Via thining parameter, you can tell the modell where are the data files, which contains the thining informations. As default it reads this from the iniFiles.
 #' @param varIndex Via mowing parameter, you can tell the modell where are the data files, which contains the mowing informations. As default it reads this from the iniFiles.
 #' @param doSensitivity Via grazing parameter, you can tell the modell where are the data files, which contains the grazing informations. As default it reads this from the iniFiles.
 #' @param onDisk Via harvesting parameter, you can tell the modell where are the data files, which contains the harvesting informations. As default it reads this from the iniFiles.
 #' @export
+
 musoMonte <- function(settings=NULL,
                      parameters,
-                      inputDir = "./",
+                     inputDir = "./",
-                      outLoc = "./calib",
+                     outLoc = "./calib",
-                      iterations = 10,
+                     iterations = 10,
-                      preTag = "mont-",
+                     preTag = "mont-",
-                      inputName = paste0(pretag,"epcs.csv"),
+                     outputType = "moreCsv",
-                      outputType = "moreCsv",
+                     fun=mean,
-                      fun=mean,
+                     varIndex = 1,
-                      varIndex=8,
+                     ...){
-                      doSensitivity=FALSE,
+    
-                      onDisk=FALSE,
-                     ...)
-{
     outLocPlain <- basename(outLoc)
     currDir <- getwd()
     inputDir <- normalizePath(inputDir)
-  tmp <- file.path(outLoc,"tmp/")
+    tmp <- file.path(outLoc,"tmp/")
 
     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"))
-  }
-
-  if(dir.exists(tmp)){
-    stop("There is a tmp directory inside the output location, please replace it. tmp is an important temporary directory for the function")
     }
-dir.create(tmp)
-  outLoc <- normalizePath(outLoc)
-  tmp <- normalizePath(tmp)
 
-  inputFiles <- file.path(inputDir,grep(basename(outLoc),list.files(inputDir),invert = TRUE,value = TRUE)) 
+    if(dir.exists(tmp)){
-
+        stop("There is a tmp directory inside the output location, please replace it. tmp is an important temporary directory for the function")
-  
-  for(i in inputFiles){
-    file.copy(i,tmp)
-  }
-  setwd(tmp)
-  
-  if(is.null(settings)){
-    settings <- setupMuso()
-  }
-  
-  parameterNames <- parameters[,1]
-  parReal <- parameters[,-1]
-  Otable <- OtableMaker(parReal)
-  A <- as.matrix(Otable[[1]][,c(2,4,5,6)])
-  B <- as.matrix(Otable[[2]])
-  settings$calibrationPar <- A[,1] 
-  pretag <- file.path(outLoc,preTag)
-  
-  ##reading the original epc file at the specified
-  ## row numbers
-  
-  origEpcFile <- readLines(settings$epcInput[2])
-  
-  origEpc <- unlist(lapply(settings$calibrationPar, function (x) {
-    as.numeric(unlist(strsplit(origEpcFile[x],split="[\t ]"))[1])
-  }))
-  
-  ## Prepare the preservedEpc matrix for the faster
-  ##  run.
-  preservedEpc <- matrix(nrow = (iterations +1 ), ncol = length(settings$calibrationPar))
-  preservedEpc[1,] <- origEpc
-  names(origEpc)<-colnames(preservedEpc)
-  colnames(preservedEpc) <- Otable[[1]][,1]
-  write.table(t(origEpc),row.names = FALSE,"preservedEpc.csv",sep=",")
-  ## Save the backupEpc, while change the settings
-  ## variable and set the output.
-  file.copy(settings$epc[2],"savedEpc",overwrite = TRUE) # do I need this? 
-  pretag <- file.path(outLoc,preTag)
-  
-  ## Creating function for generating separate
-  ## csv files for each run
-  moreCsv <- function(){
-      a <- list()
-      for(i in 1:iterations){
-          parVar <- musoRandomizer(A,B)[,2]
-                                        #preservedEpc[(i+1),] <- parVar
-          write.table(x=t(parVar),file="preservedEpc.csv",row.names=FALSE,col.names=FALSE, append=TRUE,sep=",")
-          exportName <- paste0(preTag,i,".csv")
-          tempData <- calibMuso(settings,debugging = "stamplog",
-                               parameters = parVar,
-                               keepEpc = TRUE)
-          write.csv(x=tempData,file=exportName)
-          a[[i]]<-fun(tempData[,varIndex])
-      
-      }
-      return(a)
-  }
-  
-  ## Creating function for generating one
-  ## csv files for each run
-  
-  oneCsv <- function () {
-    numDays <- settings$numdata[1]
-    if(!onDisk){
-      for(i in 1:iterations){
-      
-        parVar <- apply(parameters,1,function (x) {
-            runif(1, as.numeric(x[3]), as.numeric(x[4]))})
-          
-        preservedEpc[(i+1),] <- parVar
-        exportName <- paste0(preTag,".csv")
-        write.csv(parvar,"preservedEpc.csv",append=TRUE)
-        calibMuso(settings,debugging = "stamplog",
-                  parameters = parVar,keepEpc = TRUE) %>%
-                  {mutate(.,iD = i)} %>%
-                  {write.csv(.,file=exportName,append=TRUE)}
-      }
-      
-    return(preservedEpc)
-      } else {
-      
     }
-  }
+    dir.create(tmp)
-  
+    outLoc <- normalizePath(outLoc)
-  netCDF <- function () {
+    tmp <- normalizePath(tmp)
-    stop("This function is not inplemented yet")
+
-  }
+    inputFiles <- file.path(inputDir,grep(basename(outLoc),list.files(inputDir),invert = TRUE,value = TRUE)) 
-  
+
-  ## Call one function according to the outputType
+    
-  switch(outputType,
+    for(i in inputFiles){
-         "oneCsv" = (a <- oneCsv()),
+        file.copy(i,tmp)
-         "moreCsv" = (a <- moreCsv()),
+    }
-         "netCDF" = (a <- netCDF()))
+    setwd(tmp)
-  
+    
-  ## Change back the epc file to the original
+    if(is.null(settings)){
+        settings <- setupMuso()
+    }
+    
+    parameterNames <- parameters[,1]
+    parReal <- parameters[,-1]
+    Otable <- OtableMaker(parReal)
+    A <- as.matrix(Otable[[1]][,c(2,4,5,6)])
+    B <- as.matrix(Otable[[2]])
+    settings$calibrationPar <- A[,1] 
+    pretag <- file.path(outLoc,preTag)
+    npar <- length(settings$calibrationPar)
+    
+    ##reading the original epc file at the specified
+    ## row numbers
+    
+    origEpcFile <- readLines(settings$epcInput[2])
+    
+    origEpc <- unlist(lapply(settings$calibrationPar, function (x) {
+        as.numeric(unlist(strsplit(origEpcFile[x],split="[\t ]"))[1])
+    }))
+    
+    ## Prepare the preservedEpc matrix for the faster
+    ##  run.
+    preservedEpc <- matrix(nrow = (iterations +1 ), ncol = npar)
+    preservedEpc[1,] <- origEpc
+    colnames(preservedEpc) <- Otable[[1]][,1]
+    preservedEpc <- cbind(preservedEpc,rep(NA,(iterations+1)))
+    colnames(preservedEpc)[(npar+1)] <- "y"
+    ## Save the backupEpc, while change the settings
+    ## variable and set the output.
+    file.copy(settings$epc[2],"savedEpc",overwrite = TRUE) # do I need this? 
+    pretag <- file.path(outLoc,preTag)
+    
+    ## Creating function for generating separate
+    ## csv files for each run
+    moreCsv <- function(){
+        a <- numeric(iterations+1)
+        tempData <- calibMuso(settings, debugging = "stamplog", parameters = origEpc,keepEpc = TRUE)
+        a[1] <- tryCatch(fun(tempData[,varIndex]),error=function(e){return(NA)})
+        preservedEpc[1,(npar+1)] <- a[1]
+        write.table(t(preservedEpc[1,]),row.names = FALSE,"preservedEpc.csv",sep=",")
+        write.csv(x=tempData, file=paste0(preTag,1,".csv"))
+        for(i in 1:iterations){
+            parVar <- musoRandomizer(A,B)[,2]
+            preservedEpc[(i+1),] <- c(parVar,NA)
+            exportName <- paste0(preTag,(i+1),".csv")
+            tempData <- calibMuso(settings,debugging = "stamplog",
+                                 parameters = parVar,
+                                 keepEpc = TRUE)
+            write.csv(x=tempData,file=exportName)
+            
+            preservedEpc[(i+1),(npar+1)] <- a[i+1]<- tryCatch(fun(tempData[,varIndex]),error=function(e){return(NA)})
+            write.table(t(preservedEpc[(i+1),]),file="preservedEpc.csv",row.names=FALSE,col.names=FALSE, append=TRUE,sep=",")
+        }
+        return(preservedEpc)
+    }
+    
+    ## Creating function for generating one
+    ## csv files for each run
+    
+    oneCsv <- function () {
+        numDays <- settings$numdata[1]
+        if(!onDisk){
+            for(i in 1:iterations){
+                
+                parVar <- apply(parameters,1,function (x) {
+                    runif(1, as.numeric(x[3]), as.numeric(x[4]))})
+                
+                preservedEpc[(i+1),] <- parVar
+                exportName <- paste0(preTag,".csv")
+                write.csv(parvar,"preservedEpc.csv",append=TRUE)
+                calibMuso(settings,debugging = "stamplog",
+                          parameters = parVar,keepEpc = TRUE) %>%
+                    {mutate(.,iD = i)} %>%
+                    {write.csv(.,file=exportName,append=TRUE)}
+            }
+            
+            return(preservedEpc)
+        } else {
+            
+        }
+    }
+    
+    netCDF <- function () {
+        stop("This function is not inplemented yet")
+    }
+    
+    ## Call one function according to the outputType
+    switch(outputType,
+           "oneCsv" = (a <- oneCsv()),
+           "moreCsv" = (a <- moreCsv()),
+           "netCDF" = (a <- netCDF()))
+    
+    ## Change back the epc file to the original
     for(i in file.path("./",grep(outLocPlain, list.files(inputDir), invert = TRUE, value = TRUE))){
         file.remove(i,recursive=TRUE)
     }
@@ -154,6 +157,6 @@ dir.create(tmp)
     unlink(tmp,recursive = TRUE)
     setwd(currDir)
     file.copy("savedEpc",settings$epc[2],overwrite = TRUE)
-    return(do.call("c",a))
+    return(a)
-
 }
+

RBBGCMuso/R/musoSensi.R

@@ -0,0 +1,72 @@
+musoSensi <- function(monteCarloFile = NULL,
+                     parameters,
+                     settings = NULL,
+                     inputDir = "./",
+                     outLoc = "./calib",
+                     iterations = 30,
+                     preTag = "mount-",
+                     outputType = "moreCsv",
+                     fun = mean,
+                     varIndex = 1,
+                     outputFile = "sensitivity.csv",
+                     plotName = "sensitivity.jpg"){
+
+    if(is.null(monteCarloFile)){
+        M <- musoMonte(parameters = parameters,
+                      settings = settings,
+                      inputDir = inputDir,
+                      outLoc = outLoc,
+                      iterations = iterations,
+                      preTag = preTag,
+                      outputType = outputType,
+                      fun = fun,
+                      varIndex = varIndex
+                      )
+        npar <- ncol(M)-1
+        M %<>%
+            tbl_df() %>%
+            filter(.,!is.na(y)) %>%
+            as.data.frame()
+        y <- M[,(npar+1)]
+        M <- apply(M[,1:npar],2,function(x){x-mean(x)})
+        w <- lm(y~M)$coefficients[-1]
+        Sv <- apply(M,2,var)
+        overalVar <- sum(Sv^2*w^2)
+        S=numeric(npar)
+        for(i in 1:npar){
+            S[i] <- ((w[i]^2*Sv[i]^2)/overalVar)*100
+        }
+
+        jpg(plotName)
+        barplot(S)
+        dev.off()
+        
+        write.csv(file = outputFile, x = S)
+        barplot(S)
+        return(S)
+    } else {
+        M <- read.csv(monteCarloFile)
+        npar <- ncol(M)-1
+        M %<>%
+            tbl_df() %>%
+            filter(.,!is.na(y)) %>%
+            as.data.frame()
+        y <- M[,(npar+1)]
+        M <- apply(M[,1:npar],2,function(x){x-mean(x)})
+        w <- lm(y~M)$coefficients[-1]
+        Sv <- apply(M,2,var)
+        overalVar <- sum(Sv^2*w^2)
+        S=numeric(npar)
+        for(i in 1:npar){
+            S[i] <- ((w[i]^2*Sv[i]^2)/overalVar)*100
+        }
+
+        jpg(plotName)
+        barplot(S)
+        dev.off()
+        
+        write.csv(file = outputFile, x = S)
+        barplot(S)
+        return(S)        
+    }
+}

RBBGCMuso/man/musoMonte.Rd

@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/musoMonte.R
+\name{musoMonte}
+\alias{musoMonte}
+\title{musoMonte}
+\usage{
+musoMonte(settings = NULL, parameters, inputDir = "./",
+  outLoc = "./calib", iterations = 10, preTag = "mont-",
+  inputName = paste0(pretag, "epcs.csv"), outputType = "moreCsv",
+  fun = mean, varIndex = 8, doSensitivity = FALSE, onDisk = FALSE, ...)
+}
+\arguments{
+\item{settings}{Do you want to run multiple modell paralelly, if yes, set this variable to TRUE}
+
+\item{parameters}{stores the place of the modell-executable file. In normal usage, you don't have to be set this, because a RBBgcmuso package contains allways the latest modell executable. In spite of this, if you would like to use this package for modell development or just want to use different models (for example for comparison), you will find it useful}
+
+\item{outLoc}{Where should the modell puts its outputs. You should give a location for it via this variable, for example: outputLoc="/place/of/the/outputs/"}
+
+\item{iterations}{Usually it is the root directory, where you put the iniFiles for the modell}
+
+\item{preTag}{Via metInput parameter, you can tell the modell where are the meteorological files. As default it reads this from the iniFiles.}
+
+\item{inputName}{Via CO2 parameter, you can tell the modell where are the CO2 data files. As default it reads this from the iniFiles.}
+
+\item{outputType}{Via planting parameter, you can tell the modell where are the data files, which contains the planting informations. As default it reads this from the iniFiles.}
+
+\item{fun}{Via thining parameter, you can tell the modell where are the data files, which contains the thining informations. As default it reads this from the iniFiles.}
+
+\item{varIndex}{Via mowing parameter, you can tell the modell where are the data files, which contains the mowing informations. As default it reads this from the iniFiles.}
+
+\item{doSensitivity}{Via grazing parameter, you can tell the modell where are the data files, which contains the grazing informations. As default it reads this from the iniFiles.}
+
+\item{onDisk}{Via harvesting parameter, you can tell the modell where are the data files, which contains the harvesting informations. As default it reads this from the iniFiles.}
+
+\item{calibrationPar}{You may want to change some parameters in your epc file, before you run the modell. You have to select the appropirate modell parameters. You can refence to these with the number of the line in the epc file where the variables are. It indexes from one. You should use a vector for this, like: c(1,5,8)}
+}
+\description{
+This funcion is fundamental for the BiomBGC-MuSo modell related functions like spinupMuso, normalMuso, rungetMuso, because it sets the modells environment.
+inputDir = "./",
+}
+\author{
+Roland Hollos
+}

RBBGCMuso/tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(RBBGCMuso)
+
+test_check("RBBGCMuso")