Merge branch 'master' into Devel
This commit is contained in:
Hollos Roland
@@ -32,7 +32,10 @@ Imports:
|
||||
ncdf4,
|
||||
future,
|
||||
httr,
|
||||
tcltk
|
||||
tcltk,
|
||||
Boruta,
|
||||
rpart,
|
||||
rpart.plot
|
||||
Maintainer: Roland Hollo's <hollorol@gmail.com>
|
||||
RoxygenNote: 7.1.0
|
||||
Suggests: knitr,
|
||||
|
||||
@@ -3,17 +3,22 @@
|
||||
export(calibMuso)
|
||||
export(calibrateMuso)
|
||||
export(changemulline)
|
||||
export(checkFileSystem)
|
||||
export(checkMeteoBGC)
|
||||
export(cleanupMuso)
|
||||
export(compareMuso)
|
||||
export(copyMusoExampleTo)
|
||||
export(corrigMuso)
|
||||
export(createSoilFile)
|
||||
export(flatMuso)
|
||||
export(getAnnualOutputList)
|
||||
export(getConstMatrix)
|
||||
export(getDailyOutputList)
|
||||
export(getFilePath)
|
||||
export(getFilesFromIni)
|
||||
export(getyearlycum)
|
||||
export(getyearlymax)
|
||||
export(multiSiteCalib)
|
||||
export(musoDate)
|
||||
export(musoGlue)
|
||||
export(musoMapping)
|
||||
@@ -81,6 +86,9 @@ importFrom(magrittr,'%>%')
|
||||
importFrom(openxlsx,read.xlsx)
|
||||
importFrom(rmarkdown,pandoc_version)
|
||||
importFrom(rmarkdown,render)
|
||||
importFrom(rpart,rpart)
|
||||
importFrom(rpart,rpart.control)
|
||||
importFrom(rpart.plot,rpart.plot)
|
||||
importFrom(scales,percent)
|
||||
importFrom(stats,approx)
|
||||
importFrom(tcltk,tk_choose.files)
|
||||
|
||||
@@ -23,9 +23,13 @@
|
||||
RMuso_varTable[[version]] <<- varTable
|
||||
})
|
||||
|
||||
RMuso_depTree<- read.csv(file.path(system.file("data",package="RBBGCMuso"),"depTree.csv"), stringsAsFactors=FALSE)
|
||||
|
||||
|
||||
options(RMuso_version=RMuso_version,
|
||||
RMuso_constMatrix=RMuso_constMatrix,
|
||||
RMuso_varTable=RMuso_varTable)
|
||||
RMuso_varTable=RMuso_varTable,
|
||||
RMuso_depTree=RMuso_depTree
|
||||
)
|
||||
# getOption("RMuso_constMatrix")$soil[[as.character(getOption("RMuso_version"))]]
|
||||
}
|
||||
|
||||
@@ -56,13 +56,13 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s
|
||||
})
|
||||
|
||||
# musoSingleThread(measuredData, parameters, startDate,
|
||||
# endDate, formatString,
|
||||
# dataVar, outLoc,
|
||||
# preTag, settings,
|
||||
# outVars, iterations = threadCount[i],
|
||||
# skipSpinup, plotName,
|
||||
# modifyOriginal, likelihood, uncertainity,
|
||||
# naVal, postProcString, i)
|
||||
# endDate, formatString,
|
||||
# dataVar, outLoc,
|
||||
# preTag, settings,
|
||||
# outVars, iterations = threadCount[i],
|
||||
# skipSpinup, plotName,
|
||||
# modifyOriginal, likelihood, uncertainity,
|
||||
# naVal, postProcString, i)
|
||||
})
|
||||
})
|
||||
|
||||
|
||||
@@ -0,0 +1,267 @@
|
||||
getQueue <- function(depTree=options("RMuso_depTree")[[1]], startPoint){
|
||||
|
||||
if(length(startPoint) == 0){
|
||||
return(c())
|
||||
}
|
||||
parent <- depTree[depTree[,"name"] == startPoint,"parent"]
|
||||
c(getQueue(depTree, depTree[depTree[,"child"] == depTree[depTree[,"name"] == startPoint,"parent"],"name"]),parent)
|
||||
}
|
||||
|
||||
isRelative <- function(path){
|
||||
substr(path,1,1) != '/'
|
||||
}
|
||||
|
||||
#' getFilePath
|
||||
#'
|
||||
#' This function reads the ini file and for a chosen fileType it gives you the filePath
|
||||
#' @param iniName The name of the ini file
|
||||
#' @param filetype The type of the choosen file. For options see options("RMuso_depTree")[[1]]$name
|
||||
#' @param depTree The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]
|
||||
#' @export
|
||||
|
||||
getFilePath <- function(iniName, fileType, execPath = "./", depTree=options("RMuso_depTree")[[1]]){
|
||||
if(!file.exists(iniName) || dir.exists(iniName)){
|
||||
stop(sprintf("Cannot find iniFile: %s", iniName))
|
||||
}
|
||||
|
||||
startPoint <- fileType
|
||||
startRow <- depTree[depTree[,"name"] == startPoint,]
|
||||
startExt <- startRow$child
|
||||
|
||||
parentFile <- Reduce(function(x,y){
|
||||
tryCatch(file.path(execPath,gsub(sprintf("\\.%s.*",y),
|
||||
sprintf("\\.%s",y),
|
||||
grep(sprintf("\\.%s",y),readLines(x),value=TRUE,perl=TRUE))), error = function(e){
|
||||
stop(sprintf("Cannot find %s",x))
|
||||
})
|
||||
},
|
||||
getQueue(depTree,startPoint)[-1],
|
||||
init=iniName)
|
||||
if(startRow$mod > 0){
|
||||
tryCatch(
|
||||
gsub(sprintf("\\.%s.*", startExt),
|
||||
sprintf("\\.%s", startExt),
|
||||
grep(sprintf("\\.%s",startExt),readLines(parentFile),value=TRUE,perl=TRUE))[startRow$mod]
|
||||
,error = function(e){stop(sprintf("Cannot read %s",parentFile))})
|
||||
} else {
|
||||
res <- tryCatch(
|
||||
gsub(sprintf("\\.%s.*", startExt),
|
||||
sprintf("\\.%s",startExt),
|
||||
grep(sprintf("\\.%s",startExt),readLines(parentFile),value=TRUE, perl=TRUE))
|
||||
,error = function(e){stop(sprintf("Cannot read %s", parentFile))})
|
||||
unique(gsub(".*\\t","",res))
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
#' getFilesFromIni
|
||||
#'
|
||||
#' This function reads the ini file and gives yout back the path of all file involved in model run
|
||||
#' @param iniName The name of the ini file
|
||||
#' @param depTree The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]
|
||||
#' @export
|
||||
|
||||
getFilesFromIni <- function(iniName, execPath = "./", depTree=options("RMuso_depTree")[[1]]){
|
||||
res <- lapply(depTree$name,function(x){
|
||||
tryCatch(getFilePath(iniName,x,execPath,depTree), error = function(e){
|
||||
return(NA);
|
||||
})
|
||||
})
|
||||
names(res) <- depTree$name
|
||||
res
|
||||
}
|
||||
|
||||
#' flatMuso
|
||||
#'
|
||||
#' This function reads the ini file and creates a directory (named after the directory argument) with all the files the modell uses with this file. the directory will be flat.
|
||||
#' @param iniName The name of the ini file
|
||||
#' @param depTree The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]
|
||||
#' @param directory The destination directory for flattening. At default it will be flatdir
|
||||
#' @export
|
||||
|
||||
flatMuso <- function(iniName, execPath="./", depTree=options("RMuso_depTree")[[1]], directory="flatdir", d=TRUE,outE=TRUE){
|
||||
dir.create(directory, showWarnings=FALSE, recursive = TRUE)
|
||||
files <- getFilesFromIni(iniName,execPath,depTree)
|
||||
files <- sapply(unlist(files)[!is.na(files)], function(x){ifelse(isRelative(x),file.path(execPath,x),x)})
|
||||
file.copy(unlist(files), directory, overwrite=TRUE)
|
||||
file.copy(iniName, directory, overwrite=TRUE)
|
||||
|
||||
filesByName <- getFilesFromIni(iniName, execPath, depTree)
|
||||
for(i in seq_along(filesByName)){
|
||||
fileLines <- readLines(file.path(directory,list.files(directory, pattern = sprintf("*\\.%s", depTree$parent[i])))[1])
|
||||
|
||||
sapply(filesByName[[i]],function(origname){
|
||||
if(!is.na(origname)){
|
||||
fileLines <<- gsub(origname, basename(origname), fileLines, fixed=TRUE)
|
||||
}
|
||||
})
|
||||
|
||||
if(!is.na(filesByName[[i]][1])){
|
||||
writeLines(fileLines, file.path(directory,list.files(directory, pattern = sprintf("*\\.%s", depTree$parent[i])))[1])
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
iniLines <- readLines(file.path(directory, basename(iniName)))
|
||||
outPlace <- grep("OUTPUT_CONTROL", iniLines, perl=TRUE)+1
|
||||
if(outE){
|
||||
iniLines[outPlace] <- tools::file_path_sans_ext(basename(iniName))
|
||||
} else {
|
||||
iniLines[outPlace] <- basename(strsplit(iniLines[outPlace], split = "\\s+")[[1]][1])
|
||||
}
|
||||
if(d){
|
||||
iniLines[outPlace + 1] <- 1
|
||||
}
|
||||
writeLines(iniLines, file.path(directory, basename(iniName)))
|
||||
}
|
||||
|
||||
#' checkFileSystem
|
||||
#'
|
||||
#' This function checks the MuSo file system, if it is correct
|
||||
#' @param iniName The name of the ini file
|
||||
#' @param depTree The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]
|
||||
#' @export
|
||||
|
||||
checkFileSystem <- function(iniName,root = ".", depTree = options("RMuso_depTree")[[1]]){
|
||||
recoverAfterEval({
|
||||
setwd(root)
|
||||
fileNames <- getFilesFromIni(iniName, depTree)
|
||||
if(is.na(fileNames$management)){
|
||||
fileNames[getLeafs("management")] <- NA
|
||||
}
|
||||
fileNames <- fileNames[!is.na(fileNames)]
|
||||
errorFiles <- fileNames[!file.exists(unlist(fileNames))]
|
||||
})
|
||||
return(errorFiles)
|
||||
}
|
||||
|
||||
recoverAfterEval <- function(expr){
|
||||
wd <- getwd()
|
||||
tryCatch({
|
||||
eval(expr)
|
||||
setwd(wd)
|
||||
}, error=function(e){
|
||||
setwd(wd)
|
||||
stop(e)
|
||||
})
|
||||
}
|
||||
|
||||
getLeafs <- function(name, depTree=options("RMuso_depTree")[[1]]){
|
||||
|
||||
if(length(name) == 0){
|
||||
return(NULL)
|
||||
}
|
||||
|
||||
if(name[1] == "ini"){
|
||||
return(getLeafs(depTree$name))
|
||||
}
|
||||
|
||||
pname <- depTree[ depTree[,"name"] == name[1] , "child"]
|
||||
children <- depTree[depTree[,"parent"] == pname,"child"]
|
||||
|
||||
if(length(children)==0){
|
||||
if(length(name) == 1){
|
||||
return(NULL)
|
||||
} else{
|
||||
apname <- depTree[ depTree[,"name"] == name[2] , "child"]
|
||||
achildren <- depTree[depTree[,"parent"] == apname,"child"]
|
||||
if(length(achildren)!=0){
|
||||
return(c(name[1],name[2],getLeafs(name[-1])))
|
||||
} else{
|
||||
return(c(name[1], getLeafs(name[-1])))
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
childrenLogic <-depTree[,"child"] %in% children
|
||||
parentLogic <- depTree[,"parent"] ==pname
|
||||
res <- depTree[childrenLogic & parentLogic, "name"]
|
||||
getChildelem <- depTree[depTree[,"child"] == intersect(depTree[,"parent"], children), "name"]
|
||||
unique(c(res,getLeafs(getChildelem)))
|
||||
}
|
||||
|
||||
getParent <- function (name, depTree=options("RMuso_depTree")[[1]]) {
|
||||
parentExt <- depTree[depTree$name == name,"parent"]
|
||||
# if(length(parentExt) == 0){
|
||||
# browser()
|
||||
# }
|
||||
if(parentExt == "ini"){
|
||||
return("iniFile")
|
||||
}
|
||||
|
||||
depTree[depTree[,"child"] == parentExt,"name"]
|
||||
}
|
||||
|
||||
|
||||
|
||||
getFilePath2 <- function(iniName, fileType, depTree=options("RMuso_depTree")[[1]]){
|
||||
if(!file.exists(iniName) || dir.exists(iniName)){
|
||||
stop(sprintf("Cannot find iniFile: %s", iniName))
|
||||
}
|
||||
|
||||
startPoint <- fileType
|
||||
startRow <- depTree[depTree[,"name"] == startPoint,]
|
||||
startExt <- startRow$child
|
||||
|
||||
parentFile <- Reduce(function(x,y){
|
||||
tryCatch(gsub(sprintf("\\.%s.*",y),
|
||||
sprintf("\\.%s",y),
|
||||
grep(sprintf("\\.%s",y),readLines(x),value=TRUE,perl=TRUE)), error = function(e){
|
||||
stop(sprintf("Cannot find %s",x))
|
||||
})
|
||||
},
|
||||
getQueue(depTree,startPoint)[-1],
|
||||
init=iniName)
|
||||
res <- list()
|
||||
res["parent"] <- parentFile
|
||||
if(startRow$mod > 0){
|
||||
res["children"] <- tryCatch(
|
||||
gsub(sprintf("\\.%s.*", startExt),
|
||||
sprintf("\\.%s", startExt),
|
||||
grep(sprintf("\\.%s",startExt),readLines(parentFile),value=TRUE,perl=TRUE))[startRow$mod]
|
||||
,error = function(e){stop(sprintf("Cannot read %s",parentFile))})
|
||||
|
||||
} else {
|
||||
rows <- tryCatch(
|
||||
gsub(sprintf("\\.%s.*", startExt),
|
||||
sprintf("\\.%s",startExt),
|
||||
grep(sprintf("\\.%s",startExt),readLines(parentFile),value=TRUE, perl=TRUE))
|
||||
|
||||
,error = function(e){stop(sprintf("Cannot read %s", parentFile))})
|
||||
unique(gsub(".*\\t","",res))
|
||||
res["children"] <- unique(gsub(".*\\s+(.*\\.epc)","\\1",rows))
|
||||
}
|
||||
res
|
||||
}
|
||||
|
||||
getFilesFromIni2 <- function(iniName, depTree=options("RMuso_depTree")[[1]]){
|
||||
res <- lapply(depTree$name,function(x){
|
||||
tryCatch(getFilePath2(iniName,x,depTree), error = function(e){
|
||||
return(NA);
|
||||
})
|
||||
})
|
||||
names(res) <- depTree$name
|
||||
res
|
||||
}
|
||||
|
||||
checkFileSystemForNotif <- function(iniName,root = ".", depTree = options("RMuso_depTree")[[1]]){
|
||||
recoverAfterEval({
|
||||
setwd(root)
|
||||
fileNames <- suppressWarnings(getFilesFromIni2(iniName, depTree))
|
||||
if(is.atomic(fileNames$management)){
|
||||
fileNames[getLeafs("management")] <- NA
|
||||
}
|
||||
|
||||
hasparent <- sapply(fileNames, function(x){
|
||||
!is.atomic(x)
|
||||
})
|
||||
notNA <- ! sapply(fileNames[hasparent], function(x) {is.na(x$children)})
|
||||
errorIndex <- ! sapply(fileNames[hasparent & notNA], function(x) file.exists(x$children))
|
||||
|
||||
})
|
||||
return(fileNames[hasparent & notNA][errorIndex])
|
||||
}
|
||||
|
||||
|
||||
@@ -0,0 +1,625 @@
|
||||
`%between%` <- function(x, y){
|
||||
(x <= y[2]) & (x >= y[1])
|
||||
}
|
||||
|
||||
annualAggregate <- function(x, aggFun){
|
||||
tapply(x, rep(1:(length(x)/365), each=365), aggFun)
|
||||
}
|
||||
|
||||
SELECT <- function(x, selectPart){
|
||||
if(!is.function(selectPart)){
|
||||
index <- as.numeric(selectPart)
|
||||
tapply(x,rep(1:(length(x)/365),each=365), function(y){
|
||||
y[index]
|
||||
})
|
||||
} else {
|
||||
tapply(x,rep(1:(length(x)/365),each=365), selectPart)
|
||||
}
|
||||
}
|
||||
|
||||
bVectToInt<- function(bin_vector){
|
||||
bin_vector <- rev(as.integer(bin_vector))
|
||||
packBits(as.raw(c(bin_vector,numeric(32-length(bin_vector)))),"integer")
|
||||
}
|
||||
|
||||
constMatToDec <- function(constRes){
|
||||
tab <- table(apply(constRes,2,function(x){paste(x,collapse=" ")}))
|
||||
bitvect <- strsplit(names(tab[which.max(tab)]),split=" ")[[1]]
|
||||
bVectToInt(bitvect)
|
||||
}
|
||||
|
||||
compose <- function(expr){
|
||||
splt <- strsplit(expr,split="\\|")[[1]]
|
||||
lhs <- splt[1]
|
||||
rhs <- splt[2]
|
||||
penv <- parent.frame()
|
||||
lhsv <- eval(parse(text=lhs),envir=penv)
|
||||
penv[["lhsv"]] <- lhsv
|
||||
place <- regexpr("\\.[^0-9a-zA-Z]",rhs)
|
||||
|
||||
if(place != -1){
|
||||
finalExpression <- paste0(substr(rhs, 1, place -1),"lhsv",
|
||||
substr(rhs, place + 1, nchar(rhs)))
|
||||
} else {
|
||||
finalExpression <- paste0(rhs,"(lhsv)")
|
||||
}
|
||||
eval(parse(text=finalExpression),envir=penv)
|
||||
}
|
||||
|
||||
compoVect <- function(mod, constrTable, fileToWrite = "const_results.data"){
|
||||
with(as.data.frame(mod), {
|
||||
nexpr <- nrow(constrTable)
|
||||
filtered <- numeric(nexpr)
|
||||
vali <- numeric(nexpr)
|
||||
for(i in 1:nexpr){
|
||||
val <- compose(constrTable[i,1])
|
||||
filtered[i] <- (val <= constrTable[i,3]) &&
|
||||
(val >= constrTable[i,2])
|
||||
vali[i] <- val
|
||||
}
|
||||
|
||||
write(paste(vali,collapse=","), fileToWrite, append=TRUE)
|
||||
filtered
|
||||
})
|
||||
}
|
||||
|
||||
modCont <- function(expr, datf, interval, dumping_factor){
|
||||
tryCatch({
|
||||
if((with(datf,eval(parse(text=expr))) %between% interval)){
|
||||
return(NA)
|
||||
} else{
|
||||
return(dumping_factor)
|
||||
}
|
||||
},
|
||||
error = function(e){
|
||||
stop(sprintf("Cannot find the variable names in the dataframe, detail:\n%s",
|
||||
e))
|
||||
})
|
||||
}
|
||||
|
||||
copyToThreadDirs2 <- function(iniSource, thread_prefix = "thread", numCores, execPath="./",
|
||||
|
||||
executable = ifelse(Sys.info()[1]=="Linux", file.path(execPath, "muso"),
|
||||
file.path(execPath,"muso.exe"))){
|
||||
sapply(iniSource, function(x){
|
||||
flatMuso(x, execPath,
|
||||
directory=file.path("tmp", paste0(thread_prefix,"_1"),tools::file_path_sans_ext(basename(x)),""), d =TRUE)
|
||||
file.copy(executable,
|
||||
file.path("tmp", paste0(thread_prefix,"_1"),tools::file_path_sans_ext(basename(x))))
|
||||
tryCatch(file.copy(file.path(execPath,"cygwin1.dll"),
|
||||
file.path("tmp", paste0(thread_prefix,"_1"),tools::file_path_sans_ext(basename(x)))),
|
||||
error = function(e){"If you are in Windows..."})
|
||||
})
|
||||
sapply(2:numCores,function(thread){
|
||||
dir.create(sprintf("tmp/%s_%s",thread_prefix,thread), showWarnings=FALSE)
|
||||
file.copy(list.files(sprintf("tmp/%s_1",thread_prefix),full.names = TRUE),sprintf("tmp/%s_%s/",thread_prefix,thread),
|
||||
recursive=TRUE, overwrite = TRUE)
|
||||
})
|
||||
|
||||
}
|
||||
|
||||
|
||||
#' multiSiteCalib
|
||||
#'
|
||||
#' This funtion uses the Monte Carlo technique to uniformly sample the parameter space from user defined parameters of the Biome-BGCMuSo model. The sampling algorithm ensures that the parameters are constrained by the model logic which means that parameter dependencies are fully taken into account (parameter dependency means that e.g leaf C:N ratio must be smaller than C:N ratio of litter; more complicated rules apply to the allocation parameters where the allocation fractions to different plant compartments must sum up 1). This function implements a mathematically correct solution to provide uniform distriution of the random parameters on convex polytopes.
|
||||
#' @author Roland HOLLOS
|
||||
#' @importFrom future future
|
||||
#' @importFrom rpart rpart rpart.control
|
||||
#' @importFrom rpart.plot rpart.plot
|
||||
#' @param measuremets The table which contains the measurements
|
||||
#' @param calTable A dataframe which contantains the ini file locations and the domains they belongs to
|
||||
#' @param parameters A dataframe with the name, the minimum, and the maximum value for the parameters used in MonteCarlo experiment
|
||||
#' @param dataVar A named vector where the elements are the MuSo variable codes and the names are the same as provided in measurements and likelihood
|
||||
#' @param iterations The number of MonteCarlo experiments to be executed
|
||||
#' @param burnin Currently not used, altought it is the length of burnin period of the MCMC sampling used to generate random parameters
|
||||
#' @param likelihood A list of likelihood functions which names are linked to dataVar
|
||||
#' @param execPath If you are running the calibration from different location than the MuSo executable, you have to provide the path
|
||||
#' @param thread_prefix The prefix of thread directory names in the tmp directory created during the calibrational process
|
||||
#' @param numCores The number of processes used during the calibration. At default it uses one less than the number of threads available
|
||||
#' @param pb The progress bar function. If you use (web-)GUI you can provide a different function
|
||||
#' @param pbUpdate The update function for pb (progress bar)
|
||||
#' @param copyThread A boolean, recreate tmp directory for calibration or not (case of repeating the calibration)
|
||||
#' @param contsraints A dataframe containing the constraints logic the minimum and a maximum value for the calibration.
|
||||
#' @param th A trashold value for multisite calibration. What percentage of the site should satisfy the constraints.
|
||||
#' @param treeControl A list which controls (maximal complexity, maximal depth) the details of the decession tree making.
|
||||
#' @export
|
||||
multiSiteCalib <- function(measurements,
|
||||
calTable,
|
||||
parameters,
|
||||
dataVar,
|
||||
iterations = 100,
|
||||
burnin =ifelse(iterations < 3000, 3000, NULL),
|
||||
likelihood,
|
||||
execPath,
|
||||
thread_prefix="thread",
|
||||
numCores = (parallel::detectCores()-1),
|
||||
pb = txtProgressBar(min=0, max=iterations, style=3),
|
||||
pbUpdate = setTxtProgressBar,
|
||||
copyThread = TRUE,
|
||||
constraints=NULL, th = 10, treeControl=rpart.control()
|
||||
){
|
||||
future::plan(future::multisession)
|
||||
# file.remove(list.files(path = "tmp", pattern="progress.txt", recursive = TRUE, full.names=TRUE))
|
||||
# file.remove(list.files(path = "tmp", pattern="preservedCalib.csv", recursive = TRUE, full.names=TRUE))
|
||||
|
||||
# ____ _ _ _ _
|
||||
# / ___|_ __ ___ __ _| |_ ___ | |_| |__ _ __ ___ __ _ __| |___
|
||||
# | | | '__/ _ \/ _` | __/ _ \ | __| '_ \| '__/ _ \/ _` |/ _` / __|
|
||||
# | |___| | | __/ (_| | || __/ | |_| | | | | | __/ (_| | (_| \__ \
|
||||
# \____|_| \___|\__,_|\__\___| \__|_| |_|_| \___|\__,_|\__,_|___/
|
||||
if(copyThread){
|
||||
unlink("tmp",recursive=TRUE)
|
||||
copyToThreadDirs2(iniSource=calTable$site_id, numCores=numCores, execPath=execPath)
|
||||
} else {
|
||||
print("copy skipped")
|
||||
file.remove(file.path(list.dirs("tmp",recursive=FALSE),"progress.txt"))
|
||||
file.remove(file.path(list.dirs("tmp", recursive=FALSE), "const_results.data"))
|
||||
}
|
||||
|
||||
# ____ _ _ _
|
||||
# | _ \ _ _ _ __ | |_| |__ _ __ ___ __ _ __| |___
|
||||
# | |_) | | | | '_ \ | __| '_ \| '__/ _ \/ _` |/ _` / __|
|
||||
# | _ <| |_| | | | | | |_| | | | | | __/ (_| | (_| \__ \
|
||||
# |_| \_\\__,_|_| |_| \__|_| |_|_| \___|\__,_|\__,_|___/
|
||||
|
||||
threadCount <- distributeCores(iterations, numCores)
|
||||
fut <- lapply(1:numCores, function(i) {
|
||||
future({
|
||||
tryCatch(
|
||||
|
||||
{
|
||||
multiSiteThread(measuredData = measurements, parameters = parameters, calTable=calTable,
|
||||
dataVar = dataVar, iterations = threadCount[i],
|
||||
likelihood = likelihood, threadNumber= i, constraints=constraints, th=th)
|
||||
# setwd("../")
|
||||
}
|
||||
|
||||
, error = function(e){
|
||||
saveRDS(e,"error.RDS")
|
||||
writeLines(as.character(iterations),"progress.txt")
|
||||
})
|
||||
})
|
||||
})
|
||||
|
||||
# _ _
|
||||
# __ ____ _| |_ ___| |__ _ __ _ __ ___ __ _ _ __ ___ ___ ___
|
||||
# \ \ /\ / / _` | __/ __| '_ \ | '_ \| '__/ _ \ / _` | '__/ _ \/ __/ __|
|
||||
# \ V V / (_| | || (__| | | | | |_) | | | (_) | (_| | | | __/\__ \__ \
|
||||
# \_/\_/ \__,_|\__\___|_| |_| | .__/|_| \___/ \__, |_| \___||___/___/
|
||||
# |_| |___/
|
||||
|
||||
getProgress <- function(){
|
||||
# threadfiles <- list.files(settings$inputLoc, pattern="progress.txt", recursive = TRUE)
|
||||
threadfiles <- list.files(pattern="progress.txt", recursive = TRUE)
|
||||
if(length(threadfiles)==0){
|
||||
return(0)
|
||||
} else {
|
||||
sum(sapply(threadfiles, function(x){
|
||||
partRes <- readLines(x)
|
||||
if(length(partRes)==0){
|
||||
return(0)
|
||||
} else {
|
||||
return(as.numeric(partRes))
|
||||
}
|
||||
|
||||
}))
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
progress <- 0
|
||||
while(progress < iterations){
|
||||
Sys.sleep(1)
|
||||
progress <- tryCatch(getProgress(), error=function(e){progress})
|
||||
if(is.null(pb)){
|
||||
pbUpdate(as.numeric(progress))
|
||||
} else {
|
||||
pbUpdate(pb,as.numeric(progress))
|
||||
}
|
||||
}
|
||||
if(!is.null(pb)){
|
||||
close(pb)
|
||||
}
|
||||
|
||||
# ____ _ _
|
||||
# / ___|___ _ __ ___ | |__ (_)_ __ ___
|
||||
# | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \
|
||||
# | |__| (_) | | | | | | |_) | | | | | __/
|
||||
# \____\___/|_| |_| |_|_.__/|_|_| |_|\___|
|
||||
if(!is.null(constraints)){
|
||||
constRes <- file.path(list.dirs("tmp", recursive=FALSE), "const_results.data")
|
||||
constRes <- lapply(constRes, function(f){read.csv(f, stringsAsFactors=FALSE, header=FALSE)})
|
||||
constRes <- do.call(rbind,constRes)
|
||||
write.csv(constRes, "constRes.csv")
|
||||
}
|
||||
resultFiles <- list.files(pattern="preservedCalib.*csv$",recursive=TRUE)
|
||||
res0 <- read.csv(grep("thread_1/",resultFiles, value=TRUE),stringsAsFactors=FALSE)
|
||||
resultFilesSans0 <- grep("thread_1/", resultFiles, value=TRUE, invert=TRUE)
|
||||
# results <- do.call(rbind,lapply(resultFilesSans0, function(f){read.csv(f, stringsAsFactors=FALSE)}))
|
||||
resultsSans0 <- lapply(resultFilesSans0, function(f){read.csv(f, stringsAsFactors=FALSE, header=FALSE)})
|
||||
resultsSans0 <- do.call(rbind,resultsSans0)
|
||||
colnames(resultsSans0) <- colnames(res0)
|
||||
results <- (rbind(res0,resultsSans0))
|
||||
write.csv(results,"result.csv")
|
||||
calibrationPar <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["calibrationPar"]]
|
||||
if(!is.null(constraints)){
|
||||
notForTree <- c(seq(from = (length(calibrationPar)+1), length.out=3))
|
||||
notForTree <- c(notForTree,which(sapply(seq_along(calibrationPar),function(i){sd(results[,i])==0})))
|
||||
treeData <- results[,-notForTree]
|
||||
treeData["failType"] <- as.factor(results$failType)
|
||||
if(ncol(treeData) > 4){
|
||||
rp <- rpart(failType ~ .,data=treeData,control=treeControl)
|
||||
svg("treeplot.svg")
|
||||
tryCatch(rpart.plot(rp), error = function(e){
|
||||
print(e)
|
||||
})
|
||||
dev.off()
|
||||
}
|
||||
}
|
||||
origModOut <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["origModOut"]]
|
||||
# Just single objective version TODO:Multiobjective
|
||||
results <- results[results[,"Const"] == 1,]
|
||||
if(nrow(results)==0){
|
||||
stop("No simulation suitable for constraints\n Please see treeplot.png for explanation, if you have more than four parameters.")
|
||||
}
|
||||
bestCase <- which.max(results[,length(calibrationPar) + 1])
|
||||
parameters <- results[bestCase,1:length(calibrationPar)] # the last two column is the (log) likelihood and the rmse
|
||||
#TODO: Have to put that before multiSiteThread, we should not have to calculate it at every iterations
|
||||
|
||||
firstDir <- list.dirs("tmp/thread_1",full.names=TRUE,recursive =FALSE)[1]
|
||||
epcFile <- list.files(firstDir, pattern = "\\.epc",full.names=TRUE)
|
||||
settingsProto <- setupMuso(inputLoc = firstDir,
|
||||
iniInput =rep(list.files(firstDir, pattern = "\\.ini",full.names=TRUE),2))
|
||||
alignIndexes <- commonIndexes(settingsProto, measurements)
|
||||
musoCodeToIndex <- sapply(dataVar,function(musoCode){
|
||||
settingsProto$dailyOutputTable[settingsProto$dailyOutputTable$code == musoCode,"index"]
|
||||
})
|
||||
|
||||
|
||||
setwd("tmp/thread_1")
|
||||
aposteriori<- spatialRun(settingsProto, calibrationPar, parameters, calTable)
|
||||
file.copy(list.files(list.dirs(full.names=TRUE, recursive=FALSE)[1], pattern=".*\\.epc", full.names=TRUE),
|
||||
"../../multiSiteOptim.epc", overwrite=TRUE)
|
||||
setwd("../../")
|
||||
#TODO: Have to put that before multiSiteThread, we should not have to calculate it at every iterations
|
||||
nameGroupTable <- calTable
|
||||
nameGroupTable[,1] <- tools::file_path_sans_ext(basename(nameGroupTable[,1]))
|
||||
res <- list()
|
||||
res[["calibrationPar"]] <- calibrationPar
|
||||
res[["parameters"]] <- parameters
|
||||
res[["comparison"]] <- compareCalibratedWithOriginal(key = "grainDM", modOld=origModOut, modNew=aposteriori, mes=measurements,
|
||||
likelihoods = likelihood,
|
||||
alignIndexes = alignIndexes,
|
||||
musoCodeToIndex = musoCodeToIndex,
|
||||
nameGroupTable = nameGroupTable, mean)
|
||||
res[["likelihood"]] <- results[bestCase,ncol(results)-2]
|
||||
comp <- res$comparison
|
||||
res[["originalMAE"]] <- mean(abs((comp[,1]-comp[,3])))
|
||||
res[["MAE"]] <- mean(abs((comp[,2]-comp[,3])))
|
||||
res[["RMSE"]] <- results[bestCase,ncol(results)-2]
|
||||
res[["originalRMSE"]] <- sqrt(mean((comp[,1]-comp[,3])^2))
|
||||
res[["originalR2"]] <- summary(lm(measured ~ original,data=res$comparison))$r.squared
|
||||
res[["R2"]] <- summary(lm(measured ~ calibrated, data=res$comparison))$r.squared
|
||||
saveRDS(res,"results.RDS")
|
||||
png("calibRes.png")
|
||||
opar <- par(mar=c(5,5,4,2)+0.1, xpd=FALSE)
|
||||
with(data=res$comparison, {
|
||||
plot(measured,original,
|
||||
ylim=c(min(c(measured,original,calibrated)),
|
||||
max(c(measured,original,calibrated))),
|
||||
xlim=c(min(c(measured,original,calibrated)),
|
||||
max(c(measured,original,calibrated))),
|
||||
xlab=expression("measured "~(kg[C]~m^-2)),
|
||||
ylab=expression("simulated "~(kg[C]~m^-2)),
|
||||
cex.lab=1.3,
|
||||
col="red",
|
||||
pch=19,
|
||||
pty="s"
|
||||
)
|
||||
points(measured,calibrated, pch=19, col="blue")
|
||||
abline(0,1)
|
||||
legend(x="top",
|
||||
pch=c(19,19),
|
||||
col=c("red","blue"),
|
||||
inset=c(0,-0.1),
|
||||
legend=c("original","calibrated"),
|
||||
ncol=2,
|
||||
box.lty=0,
|
||||
xpd=TRUE
|
||||
)
|
||||
})
|
||||
dev.off()
|
||||
return(res)
|
||||
}
|
||||
|
||||
#' multiSiteThread
|
||||
#'
|
||||
#' This is an
|
||||
#' @author Roland HOLLOS
|
||||
|
||||
|
||||
multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
|
||||
endDate = NULL, formatString = "%Y-%m-%d", calTable,
|
||||
dataVar, outLoc = "./calib",
|
||||
outVars = NULL, iterations = 300,
|
||||
skipSpinup = TRUE, plotName = "calib.jpg",
|
||||
modifyOriginal=TRUE, likelihood, uncertainity = NULL, burnin=NULL,
|
||||
naVal = NULL, postProcString = NULL, threadNumber, constraints=NULL,th=10) {
|
||||
|
||||
originalRun <- list()
|
||||
nameGroupTable <- calTable
|
||||
nameGroupTable[,1] <- tools::file_path_sans_ext(basename(nameGroupTable[,1]))
|
||||
setwd(paste0("tmp/thread_",threadNumber))
|
||||
firstDir <- list.dirs(full.names=FALSE,recursive =FALSE)[1]
|
||||
epcFile <- list.files(firstDir, pattern = "\\.epc",full.names=TRUE)
|
||||
settingsProto <- setupMuso(inputLoc = firstDir,
|
||||
iniInput =rep(list.files(firstDir, pattern = "\\.ini",full.names=TRUE),2))
|
||||
|
||||
# Exanding likelihood
|
||||
likelihoodFull <- as.list(rep(NA,length(dataVar)))
|
||||
names(likelihoodFull) <- names(dataVar)
|
||||
if(!missing(likelihood)) {
|
||||
lapply(names(likelihood),function(x){
|
||||
likelihoodFull[[x]] <<- likelihood[[x]]
|
||||
})
|
||||
}
|
||||
|
||||
defaultLikelihood <- which(is.na(likelihood))
|
||||
if(length(defaultLikelihood)>0){
|
||||
likelihoodFull[[defaultLikelihood]] <- (function(x, y){
|
||||
exp(-sqrt(mean((x-y)^2)))
|
||||
})
|
||||
}
|
||||
|
||||
mdata <- measuredData
|
||||
if(is.null(parameters)){
|
||||
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.")
|
||||
})
|
||||
} else {
|
||||
if((!is.list(parameters)) & (!is.matrix(parameters))){
|
||||
parameters <- tryCatch(read.csv(parameters, stringsAsFactor=FALSE), error = function (e){
|
||||
stop("Cannot find neither parameters file neither the parameters matrix")
|
||||
})
|
||||
}}
|
||||
|
||||
print("optiMuso is randomizing the epc parameters now...",quote = FALSE)
|
||||
randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations)
|
||||
|
||||
origEpc <- readValuesFromFile(epcFile, randVals[[1]])
|
||||
partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar) + 2)
|
||||
colN <- randVals[[1]]
|
||||
colN[match(parameters[,2],randVals[[1]])] <- parameters[,1]
|
||||
colN[match(parameters[,2], randVals[[1]])[!is.na(match(parameters[,2],randVals[[1]]))]] <- parameters[,1]
|
||||
colnames(partialResult) <- c(colN,sprintf("%s_likelihood",names(dataVar)),
|
||||
sprintf("%s_rmse",names(dataVar)),"Const", "failType")
|
||||
numParameters <- length(colN)
|
||||
partialResult[1:numParameters] <- origEpc
|
||||
## Prepare the preservedCalib matrix for the faster
|
||||
## run.
|
||||
musoCodeToIndex <- sapply(dataVar,function(musoCode){
|
||||
settingsProto$dailyOutputTable[settingsProto$dailyOutputTable$code == musoCode,"index"]
|
||||
})
|
||||
|
||||
resultRange <- (numParameters + 1):(ncol(partialResult))
|
||||
randValues <- randVals[[2]]
|
||||
|
||||
settingsProto$calibrationPar <- randVals[[1]]
|
||||
|
||||
if(!is.null(naVal)){
|
||||
measuredData <- as.data.frame(measuredData)
|
||||
measuredData[measuredData == naVal] <- NA
|
||||
}
|
||||
resIterate <- 1:nrow(calTable)
|
||||
names(resIterate) <- tools::file_path_sans_ext(basename(calTable[,1]))
|
||||
alignIndexes <- commonIndexes(settingsProto, measuredData)
|
||||
if(threadNumber == 1){
|
||||
originalRun[["calibrationPar"]] <- randVals[[1]]
|
||||
|
||||
origModOut <- lapply(resIterate, function(i){
|
||||
dirName <- tools::file_path_sans_ext(basename(calTable[i,1]))
|
||||
setwd(dirName)
|
||||
settings <- settingsProto
|
||||
settings$outputLoc <- settings$inputLoc <- "./"
|
||||
settings$iniInput <- settings$inputFiles <- rep(paste0(dirName,".ini"),2)
|
||||
settings$outputNames <- rep(dirName,2)
|
||||
settings$executable <- ifelse(Sys.info()[1]=="Linux","./muso","./muso.exe") # set default exe option at start wold be better
|
||||
res <- tryCatch(calibMuso(settings=settings,parameters =origEpc, silent = TRUE, skipSpinup = TRUE), error=function(e){NA})
|
||||
setwd("../")
|
||||
res
|
||||
})
|
||||
originalRun[["origModOut"]] <- origModOut
|
||||
|
||||
partialResult[,resultRange] <- calcLikelihoodsForGroups(dataVar=dataVar,
|
||||
mod=origModOut,
|
||||
mes=measuredData,
|
||||
likelihoods=likelihood,
|
||||
alignIndexes=alignIndexes,
|
||||
musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, groupFun=mean, constraints=constraints,th=th)
|
||||
|
||||
write.csv(x=randVals[[1]],"../randIndexes.csv")
|
||||
write.csv(x=partialResult, file="preservedCalib.csv",row.names=FALSE)
|
||||
}
|
||||
|
||||
print("Running the model with the random epc values...", quote = FALSE)
|
||||
for(i in 2:(iterations+1)){
|
||||
tmp <- lapply(resIterate, function(siteI){
|
||||
dirName <- tools::file_path_sans_ext(basename(calTable[siteI,1]))
|
||||
setwd(dirName)
|
||||
settings <- settingsProto
|
||||
settings$outputLoc <- settings$inputLoc <- "./"
|
||||
settings$iniInput <- settings$inputFiles <- rep(paste0(dirName,".ini"),2)
|
||||
settings$outputNames <- rep(dirName,2)
|
||||
settings$executable <- ifelse(Sys.info()[1]=="Linux","./muso","./muso.exe") # set default exe option at start wold be better
|
||||
|
||||
res <- tryCatch(calibMuso(settings=settings,parameters=randValues[(i-1),], silent = TRUE, skipSpinup = TRUE), error=function(e){NA})
|
||||
setwd("../")
|
||||
res
|
||||
})
|
||||
|
||||
if(is.null(tmp)){
|
||||
partialResult[,resultRange] <- NA
|
||||
} else {
|
||||
partialResult[,resultRange] <- calcLikelihoodsForGroups(dataVar=dataVar,
|
||||
mod=tmp,
|
||||
mes=measuredData,
|
||||
likelihoods=likelihood,
|
||||
alignIndexes=alignIndexes,
|
||||
musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, groupFun=mean, constraints = constraints, th=th)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
partialResult[1:numParameters] <- randValues[(i-1),]
|
||||
write.table(x=partialResult, file="preservedCalib.csv", append=TRUE, row.names=FALSE,
|
||||
sep=",", col.names=FALSE)
|
||||
# write.csv(x=tmp, file=paste0(pretag, (i+1),".csv"))
|
||||
writeLines(as.character(i-1),"progress.txt") #UNCOMMENT IMPORTANT
|
||||
}
|
||||
}
|
||||
if(threadNumber == 1){
|
||||
return(originalRun)
|
||||
}
|
||||
}
|
||||
distributeCores <- function(iterations, numCores){
|
||||
perProcess<- iterations %/% numCores
|
||||
numSimu <- rep(perProcess,numCores)
|
||||
gainers <- sample(1:numCores, iterations %% numCores)
|
||||
numSimu[gainers] <- numSimu[gainers] + 1
|
||||
numSimu
|
||||
}
|
||||
|
||||
prepareFromAgroMo <- function(fName){
|
||||
obs <- read.table(fName, stringsAsFactors=FALSE, sep = ";", header=T)
|
||||
obs <- reshape(obs, timevar="var_id", idvar = "date", direction = "wide")
|
||||
dateCols <- apply(do.call(rbind,(strsplit(obs$date, split = "-"))),2,as.numeric)
|
||||
colnames(dateCols) <- c("year", "month", "day")
|
||||
cbind.data.frame(dateCols, obs)
|
||||
}
|
||||
|
||||
calcLikelihoodsForGroups <- function(dataVar, mod, mes,
|
||||
likelihoods, alignIndexes, musoCodeToIndex,
|
||||
nameGroupTable, groupFun, constraints,
|
||||
th = 10){
|
||||
|
||||
if(!is.null(constraints)){
|
||||
constRes<- sapply(mod,function(m){
|
||||
compoVect(m,constraints)
|
||||
})
|
||||
|
||||
failType <- constMatToDec(constRes)
|
||||
}
|
||||
|
||||
likelihoodRMSE <- sapply(names(dataVar),function(key){
|
||||
modelled <- as.vector(unlist(sapply(sort(names(alignIndexes)),
|
||||
function(domain_id){
|
||||
apply(do.call(cbind,
|
||||
lapply(nameGroupTable[,1][nameGroupTable[,2] == domain_id],
|
||||
function(site){mod[[site]][alignIndexes[[domain_id]]$model,musoCodeToIndex[key]]
|
||||
})),1,groupFun)
|
||||
|
||||
|
||||
|
||||
|
||||
})))
|
||||
|
||||
|
||||
measuredGroups <- split(mes,mes$domain_id)
|
||||
measured <- do.call(rbind.data.frame, lapply(names(measuredGroups), function(domain_id){
|
||||
measuredGroups[[domain_id]][alignIndexes[[domain_id]]$meas,]
|
||||
}))
|
||||
measured <- measured[measured$var_id == key,]
|
||||
|
||||
res <- c(likelihoods[[key]](modelled, measured),
|
||||
sqrt(mean((modelled-measured$mean)^2))
|
||||
)
|
||||
print(abs(mean(modelled)-mean(measured$mean)))
|
||||
res
|
||||
})
|
||||
|
||||
likelihoodRMSE <- c(likelihoodRMSE[1,], likelihoodRMSE[2,],
|
||||
ifelse((100 * sum(apply(constRes, 2, prod)) / ncol(constRes)) >= th,
|
||||
1,0), failType)
|
||||
names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar), "Const", "failType")
|
||||
return(likelihoodRMSE)
|
||||
}
|
||||
|
||||
commonIndexes <- function (settings,measuredData) {
|
||||
# Have to fix for other starting points also
|
||||
modelDates <- seq(from= as.Date(sprintf("%s-01-01",settings$startYear)),
|
||||
by="days",
|
||||
to=as.Date(sprintf("%s-12-31",settings$startYear+settings$numYears-1)))
|
||||
modelDates <- grep("-02-29",modelDates,invert=TRUE, value=TRUE)
|
||||
|
||||
lapply(split(measuredData,measuredData$domain_id),function(x){
|
||||
measuredDates <- x$date
|
||||
modIndex <- match(as.Date(measuredDates), as.Date(modelDates))
|
||||
measIndex <- which(!is.na(modIndex))
|
||||
modIndex <- modIndex[!is.na(modIndex)]
|
||||
cbind.data.frame(model=modIndex,meas=measIndex)
|
||||
})
|
||||
}
|
||||
|
||||
agroLikelihood <- function(modVector,measured){
|
||||
mu <- measured[,grep("mean", colnames(measured))]
|
||||
stdev <- measured[,grep("^sd", colnames(measured))]
|
||||
ndata <- nrow(measured)
|
||||
sum(sapply(1:ndata, function(x){
|
||||
dnorm(modVector, mu[x], stdev[x], log = TRUE)
|
||||
}), na.rm=TRUE)
|
||||
}
|
||||
|
||||
|
||||
#' compareCalibratedWithOriginal
|
||||
#'
|
||||
#' This functions compareses the likelihood and the RMSE values of the simulations and the measurements
|
||||
#' @param key
|
||||
compareCalibratedWithOriginal <- function(key, modOld, modNew, mes,
|
||||
likelihoods, alignIndexes, musoCodeToIndex, nameGroupTable,
|
||||
groupFun){
|
||||
|
||||
original <- as.vector(unlist(sapply(sort(names(alignIndexes)),
|
||||
function(domain_id){
|
||||
apply(do.call(cbind,
|
||||
lapply(nameGroupTable$site_id[nameGroupTable$domain_id == domain_id],
|
||||
function(site){
|
||||
modOld[[site]][alignIndexes[[domain_id]]$model,musoCodeToIndex[key]]
|
||||
})),1,groupFun)
|
||||
})))
|
||||
calibrated <- as.vector(unlist(sapply(sort(names(alignIndexes)),
|
||||
function(domain_id){
|
||||
apply(do.call(cbind,
|
||||
lapply(nameGroupTable$site_id[nameGroupTable$domain_id == domain_id],
|
||||
function(site){
|
||||
modNew[[site]][alignIndexes[[domain_id]]$model,musoCodeToIndex[key]]
|
||||
})),1,groupFun)
|
||||
})))
|
||||
measuredGroups <- split(mes,mes$domain_id)
|
||||
measured <- do.call(rbind.data.frame, lapply(names(measuredGroups), function(domain_id){
|
||||
measuredGroups[[domain_id]][alignIndexes[[domain_id]]$meas,]
|
||||
}))
|
||||
measured <- measured[measured$var_id == key,]
|
||||
return(data.frame(original = original, calibrated = calibrated,measured=measured$mean))
|
||||
}
|
||||
|
||||
|
||||
spatialRun <- function(settingsProto,calibrationPar, parameters, calTable){
|
||||
resIterate <- 1:nrow(calTable)
|
||||
names(resIterate) <- tools::file_path_sans_ext(basename(calTable[,1]))
|
||||
modOut <- lapply(resIterate, function(i){
|
||||
dirName <- tools::file_path_sans_ext(basename(calTable[i,1]))
|
||||
setwd(dirName)
|
||||
settings <- settingsProto
|
||||
settings$outputLoc <- settings$inputLoc <- "./"
|
||||
settings$iniInput <- settings$inputFiles <- rep(paste0(dirName,".ini"),2)
|
||||
settings$outputNames <- rep(dirName,2)
|
||||
settings$calibrationPar <- calibrationPar
|
||||
settings$executable <- ifelse(Sys.info()[1]=="Linux","./muso","./muso.exe") # set default exe option at start wold be better
|
||||
res <- tryCatch(calibMuso(settings=settings,parameters =parameters, silent = TRUE, skipSpinup = TRUE), error=function(e){NA})
|
||||
setwd("../")
|
||||
res
|
||||
})
|
||||
modOut
|
||||
}
|
||||
@@ -8,7 +8,7 @@
|
||||
#' @importFrom limSolve xsample
|
||||
#' @export
|
||||
|
||||
musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL){
|
||||
musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL, burnin = NULL){
|
||||
if(is.null(constrains)){
|
||||
constMatrix <- constrains
|
||||
constMatrix <- getOption("RMuso_constMatrix")[[fileType]][[as.character(getOption("RMuso_version"))]]
|
||||
@@ -176,7 +176,7 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N
|
||||
E <- do.call(rbind,lapply(Ef,function(x){x$E}))
|
||||
f <- do.call(c,lapply(Ef,function(x){x$f}))
|
||||
# browser()
|
||||
randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,iter = iterations))$X
|
||||
randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,burninlength=burnin, iter = iterations))$X
|
||||
} else{
|
||||
Gh0<-genMat0(dependences)
|
||||
randVal <- suppressWarnings(xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
#' updateMusoMapping
|
||||
#'
|
||||
#' This function updates the Biome-BGCMuSo output code-variable matrix. Within Biome-BGCMuSo the state variables and fluxes are marked by integer numbers. In order to provide meaningful variable names (e.g. 3009 means Gross Primary Production in Biome-BGCMuSo v5) a conversion table is needed which is handled by this function.
|
||||
#' This function updates the Biome-BGCMuSo output code-variable matrix (creates a json file that is used internally by RBBGCMuso). Within Biome-BGCMuSo the output state variablesare marked by integer numbers (see the User's Guide). In order to provide meaningful variable names (e.g. 3009 means Gross Primary Production) a conversion table is needed which is handled by this function. The input Excel file must have the following column order: name, index, units, description (plus other optional columns line group). name refers to the abbreviation of the variable; index is the integer number of the output variable; unit is the unit of the variable; description is a meaningful text to explain the variable. The script will NOT work with other column order!
|
||||
#' @author Roland HOLLOS
|
||||
#' @param excelName Name of the excelfile which contains the parameters
|
||||
#' @importFrom openxlsx read.xlsx
|
||||
|
||||
@@ -1,37 +0,0 @@
|
||||
---
|
||||
title: "An easy grouping algorithm"
|
||||
author: "Hollós Roland"
|
||||
date: "10/1/2019"
|
||||
output: pdf_document
|
||||
---
|
||||
|
||||
|
||||
|
||||
## R Markdown
|
||||
|
||||
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.
|
||||
|
||||
When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
|
||||
|
||||
|
||||
```r
|
||||
summary(cars)
|
||||
```
|
||||
|
||||
```
|
||||
## speed dist
|
||||
## Min. : 4.0 Min. : 2.00
|
||||
## 1st Qu.:12.0 1st Qu.: 26.00
|
||||
## Median :15.0 Median : 36.00
|
||||
## Mean :15.4 Mean : 42.98
|
||||
## 3rd Qu.:19.0 3rd Qu.: 56.00
|
||||
## Max. :25.0 Max. :120.00
|
||||
```
|
||||
|
||||
## Including Plots
|
||||
|
||||
You can also embed plots, for example:
|
||||
|
||||
|
||||
|
||||
Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 5.4 KiB |
@@ -1,27 +0,0 @@
|
||||
#' getDailyOutputList
|
||||
#'
|
||||
#' bla bla
|
||||
#' @param settings bla
|
||||
#' @export
|
||||
|
||||
|
||||
getDailyOutputList <- function(settings=NULL){
|
||||
if(is.null(settings)){
|
||||
settings<- setupMuso()
|
||||
}
|
||||
settings$dailyOutputTable
|
||||
}
|
||||
|
||||
#' getAnnualOutputList
|
||||
#'
|
||||
#' bla bla
|
||||
#' @param settings bla
|
||||
#' @export
|
||||
|
||||
|
||||
getAnnualOutputList <- function(settings=NULL){
|
||||
if(is.null(settings)){
|
||||
settings<- setupMuso()
|
||||
}
|
||||
settings$annualOutputTable
|
||||
}
|
||||
@@ -0,0 +1,18 @@
|
||||
"child","parent","mod","name"
|
||||
"wth","ini",1,"weather"
|
||||
"endpoint","ini",1,"endpointIn"
|
||||
"endpoint","ini",2,"endpointOut"
|
||||
"txt","ini",1,"co2"
|
||||
"txt","ini",2,"nitrogen"
|
||||
"soi","ini",1,"soil"
|
||||
"epc","ini",1,"startEpc"
|
||||
"mgm","ini",1,"management"
|
||||
"plt","mgm",1,"planting"
|
||||
"thn","mgm",1,"thining"
|
||||
"mow","mgm",1,"mowing"
|
||||
"grz","mgm",1,"grazing"
|
||||
"hrv","mgm",1,"harvest"
|
||||
"cul","mgm",1,"cultivation"
|
||||
"frz","mgm",1,"fertilization"
|
||||
"irr","mgm",1,"irrigation"
|
||||
"epc","plt",0,"plantEpc"
|
||||
|
@@ -865,6 +865,7 @@
|
||||
"UNIT": "prop",
|
||||
"MIN": 0,
|
||||
"MAX": 0.1,
|
||||
"DEPENDENCE": 0,
|
||||
"GROUP": 0,
|
||||
"TYPE": 0
|
||||
},
|
||||
@@ -875,6 +876,7 @@
|
||||
"UNIT": "prop",
|
||||
"MIN": 0,
|
||||
"MAX": 0.1,
|
||||
"DEPENDENCE": 1,
|
||||
"GROUP": 0,
|
||||
"TYPE": 0
|
||||
},
|
||||
|
||||
+4874
-4635
File diff suppressed because it is too large
Load Diff
@@ -8,20 +8,20 @@ FLAGS
|
||||
PLANT FUNCTIONING PARAMETERS
|
||||
0 (yday) yearday to start new growth (when phenology flag = 0)
|
||||
364 (yday) yearday to end litterfall (when phenology flag = 0)
|
||||
1.0 (prop.) transfer growth period as fraction of growing season (when transferGDD_flag = 0)
|
||||
1.0 (prop.) litterfall as fraction of growing season (when transferGDD_flag = 0)
|
||||
0.5 (prop.) transfer growth period as fraction of growing season (when transferGDD_flag = 0)
|
||||
0.5 (prop.) litterfall as fraction of growing season (when transferGDD_flag = 0)
|
||||
0 (Celsius) base temperature
|
||||
-9999 (Celsius) minimum temperature for growth displayed on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) optimal1 temperature for growth displayed on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) optimal2 temperature for growth displayed on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) maxmimum temperature for growth displayed on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) minimum temperature for carbon assimilation on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) optimal1 temperature for carbon assimilation on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) optimal2 temperature for carbon assimilation on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) maxmimum temperature for carbon assimilation on current day (-9999: no T-dependence)
|
||||
-9999 (Celsius) minimum temperature for growth displayed on current day (-9999: no T-dependence of allocation)
|
||||
-9999 (Celsius) optimal1 temperature for growth displayed on current day (-9999: no T-dependence of allocation)
|
||||
-9999 (Celsius) optimal2 temperature for growth displayed on current day (-9999: no T-dependence of allocation)
|
||||
-9999 (Celsius) maxmimum temperature for growth displayed on current day (-9999: no T-dependence of allocation)
|
||||
-9999 (Celsius) minimum temperature for carbon assimilation displayed on current day (-9999: no limitation)
|
||||
-9999 (Celsius) optimal1 temperature for carbon assimilation displayed on current day (-9999: no limitation)
|
||||
-9999 (Celsius) optimal2 temperature for carbon assimilation displayed on current day (-9999: no limitation)
|
||||
-9999 (Celsius) maxmimum temperature for carbon assimilation displayed on current day (-9999: no limitation)
|
||||
1.0 (1/yr) annual leaf and fine root turnover fraction
|
||||
0.00 (1/yr) annual live wood turnover fraction
|
||||
0.0 (1/yr) annual fire mortality fraction
|
||||
0.03 (1/yr) annual fire mortality fraction
|
||||
0.01 (1/vegper) whole-plant mortality fraction in vegetation period
|
||||
36.6 (kgC/kgN) C:N of leaves
|
||||
45.0 (kgC/kgN) C:N of leaf litter, after retranslocation
|
||||
@@ -47,30 +47,30 @@ PLANT FUNCTIONING PARAMETERS
|
||||
0.23 (DIM) soft stem litter cellulose proportion
|
||||
0.00 *(DIM) dead wood cellulose proportion
|
||||
0.01 (1/LAI/d) canopy water interception coefficient
|
||||
0.7 (DIM) canopy light extinction coefficient
|
||||
0.6 (g/MJ) potential radiation use efficiency
|
||||
0.63 (DIM) canopy light extinction coefficient
|
||||
2.0 (g/MJ) potential radiation use efficiency
|
||||
0.781 (DIM) radiation parameter1 (Jiang et al.2015)
|
||||
-13.596 (DIM) radiation parameter2 (Jiang et al.2015)
|
||||
2.0 (DIM) all-sided to projected leaf area ratio
|
||||
2.0 (DIM) ratio of shaded SLA:sunlit SLA
|
||||
0.3 (DIM) fraction of leaf N in Rubisco
|
||||
0.14 (DIM) fraction of leaf N in Rubisco
|
||||
0.03 (DIM) fraction of leaf N in PEP Carboxylase
|
||||
0.002 (m/s) maximum stomatal conductance (projected area basis)
|
||||
0.004 (m/s) maximum stomatal conductance (projected area basis)
|
||||
0.00006 (m/s) cuticular conductance (projected area basis)
|
||||
0.039 (m/s) boundary layer conductance (projected area basis)
|
||||
1.0 (m) maximum height of plant
|
||||
0.04 (m/s) boundary layer conductance (projected area basis)
|
||||
1.5 (m) maximum height of plant
|
||||
0.8 (kgC) stem weight corresponding to maximum height
|
||||
0.5 (dimless) plant height function shape parameter (slope)
|
||||
1.0 (m) maximum depth of rooting zone
|
||||
2.67 (DIM) root distribution parameter
|
||||
4.0 (m) maximum depth of rooting zone
|
||||
3.67 (DIM) root distribution parameter
|
||||
0.4 (kgC) root weight corresponding to max root depth
|
||||
0.5 (dimless) root depth function shape parameter (slope)
|
||||
1000 (m/kg) root weight to rooth length conversion factor
|
||||
1000 (m/kg) root weight to root length conversion factor
|
||||
0.3 (prop.) growth resp per unit of C grown
|
||||
0.195 (kgC/kgN/d) maintenance respiration in kgC/day per kg of tissue N
|
||||
0.218 (kgC/kgN/d) maintenance respiration in kgC/day per kg of tissue N
|
||||
0.1 (DIM) theoretical maximum prop. of non-structural and structural carbohydrates
|
||||
0.24 (DIM) prop. of non-structural carbohydrates available for maintanance respiration
|
||||
0.0048 (kgN/m2/yr) symbiotic+asymbiotic fixation of N
|
||||
0.02 (kgN/m2/yr) symbiotic+asymbiotic fixation of N
|
||||
0 (day) time delay for temperature in photosynthesis acclimation
|
||||
----------------------------------------------------------------------------------------
|
||||
CROP SPECIFIC PARAMETERS
|
||||
@@ -93,21 +93,21 @@ CROP SPECIFIC PARAMETERS
|
||||
0.2 (prop.) theoretical maximum of flowering thermal stress mortality parameter
|
||||
----------------------------------------------------------------------------------------
|
||||
STRESS AND SENESCENCE PARAMETERS
|
||||
1.0 (prop) VWC ratio to calc. soil moisture limit 1 (prop. to FC-WP)
|
||||
0.95 (prop) VWC ratio to calc. soil moisture limit 2 (prop. to SAT-FC)
|
||||
0.98 (prop) VWC ratio to calc. soil moisture limit 1 (prop. to FC-WP)
|
||||
0.7 (prop) VWC ratio to calc. soil moisture limit 2 (prop. to SAT-FC)
|
||||
0.4 (prop) minimum of soil moisture limit2 multiplicator (full anoxic stress value)
|
||||
1000 (Pa) vapor pressure deficit: start of conductance reduction
|
||||
2800 (Pa) vapor pressure deficit: complete conductance reduction
|
||||
0.03 (prop.) maximum senescence mortality coefficient of aboveground plant material
|
||||
0.02 (prop.) maximum senescence mortality coefficient of belowground plant material
|
||||
0.01 (prop.) maximum senescence mortality coefficient of non-structured plant material
|
||||
4000 (Pa) vapor pressure deficit: complete conductance reduction
|
||||
0.003 (prop.) maximum senescence mortality coefficient of aboveground plant material
|
||||
0.001 (prop.) maximum senescence mortality coefficient of belowground plant material
|
||||
0.0 (prop.) maximum senescence mortality coefficient of non-structured plant material
|
||||
35 (Celsius) lower limit extreme high temperature effect on senescence mortality
|
||||
40 (Celsius) upper limit extreme high temperature effect on senescence mortality
|
||||
0.01 (prop.) turnover rate of wilted standing biomass to litter
|
||||
0.047 (prop.) turnover rate of non-woody cut-down biomass to litter
|
||||
0.01 (prop.) turnover rate of woody cut-down biomass to litter
|
||||
30 (nday) drought tolerance parameter (critical value of DSWS)
|
||||
0.2 (dimless) effect of soilstress factor on photosynthesis (1: full effect, 0: no effect)
|
||||
17 (nday) drought tolerance parameter (critical value of DSWS)
|
||||
0.3 (prop) soil water deficit effect on photosynthesis downregulation
|
||||
----------------------------------------------------------------------------------------
|
||||
GROWING SEASON PARAMETERS
|
||||
5 (kg/m2) crit. amount of snow limiting photosyn.
|
||||
@@ -125,9 +125,9 @@ GROWING SEASON PARAMETERS
|
||||
----------------------------------------------------------------------------------------
|
||||
PHENOLOGICAL (ALLOCATION) PARAMETERS (7 phenological phases)
|
||||
phase1 phase2 phase3 phase4 phase5 phase6 phase7 (text) name of the phenophase
|
||||
10000 200 500 200 400 200 100 (Celsius) length of phenophase (GDD)
|
||||
0.3 0.3 0.3 0.3 0.3 0.3 0.3 (ratio) leaf ALLOCATION
|
||||
0.5 0.5 0.5 0.5 0.5 0.5 0.5 (ratio) fine root ALLOCATION
|
||||
5000 200 500 200 400 200 100 (Celsius) length of phenophase (GDD)
|
||||
0.3 0.4 0.4 0.4 0.4 0.4 0.4 (ratio) leaf ALLOCATION
|
||||
0.5 0.4 0.4 0.4 0.4 0.4 0.4 (ratio) fine root ALLOCATION
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 (ratio) fruit ALLOCATION
|
||||
0.2 0.2 0.2 0.2 0.2 0.2 0.2 (ratio) soft stem ALLOCATION
|
||||
0 0 0 0 0 0 0 (ratio) live woody stem ALLOCATION
|
||||
@@ -136,4 +136,4 @@ phase1 phase2 phase3 phase4 phase5 phase6 phase7 (text) name of the pheno
|
||||
0 0 0 0 0 0 0 (ratio) dead coarse root ALLOCATION
|
||||
49 49 49 49 49 49 49 (m2/kgC) canopy average specific leaf area (projected area basis)
|
||||
0.37 0.37 0.37 0.37 0.37 0.37 0.37 (prop.) current growth proportion
|
||||
10000 10000 10000 10000 10000 10000 10000 (Celsius) maximal lifetime of plant tissue
|
||||
10000 10000 10000 10000 10000 10000 10000 (Celsius) maximal lifetime of plant tissue
|
||||
|
||||
@@ -6,7 +6,7 @@ NITROGEN AND DECOMPOSITION PARAMETERS
|
||||
0.1 (prop.) nitrification coefficient 2
|
||||
0.02 (prop.) coefficient of N2O emission of nitrification
|
||||
0.1 (prop.) NH4 mobilen proportion
|
||||
1.0 (prop.) NO3 mobilen proportion
|
||||
1.0 denitrification related N2/N2O ratio multiplier (soil texture effect)
|
||||
10 (m) e-folding depth of decomposition rate's depth scalar
|
||||
0.002 (prop.) fraction of dissolved part of SOIL1 organic matter
|
||||
0.002 (prop.) fraction of dissolved part of SOIL2 organic matter
|
||||
@@ -16,7 +16,7 @@ NITROGEN AND DECOMPOSITION PARAMETERS
|
||||
0.45 (prop.) lower optimum WFPS for scalar of nitrification calculation
|
||||
0.55 (prop.) higher optimum WFPS for scalar of nitrification calculation
|
||||
0.2 (prop.) minimum value for saturated WFPS scalar of nitrification calculation
|
||||
10 (ppm) critical value of dissolved N and C in bottom (inactive layer)
|
||||
10 (ppm) C:N ratio of recaltirant SOM (slowest)
|
||||
----------------------------------------------------------------------------------------
|
||||
RATE SCALARS
|
||||
0.39 (DIM) respiration fractions for fluxes between compartments (l1s1)
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -6,7 +6,7 @@
|
||||
\usage{
|
||||
calibrateMuso(
|
||||
measuredData,
|
||||
parameters = NULL,
|
||||
parameters = read.csv("parameters.csv", stringsAsFactor = FALSE),
|
||||
startDate = NULL,
|
||||
endDate = NULL,
|
||||
formatString = "\%Y-\%m-\%d",
|
||||
@@ -28,6 +28,7 @@ calibrateMuso(
|
||||
pb = txtProgressBar(min = 0, max = iterations, style = 3),
|
||||
maxLikelihoodEpc = TRUE,
|
||||
pbUpdate = setTxtProgressBar,
|
||||
outputLoc = "./",
|
||||
method = "GLUE",
|
||||
lg = FALSE,
|
||||
w = NULL,
|
||||
|
||||
@@ -0,0 +1,16 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/flat.R
|
||||
\name{checkFileSystem}
|
||||
\alias{checkFileSystem}
|
||||
\title{checkFileSystem}
|
||||
\usage{
|
||||
checkFileSystem(iniName, root = ".", depTree = options("RMuso_depTree")[[1]])
|
||||
}
|
||||
\arguments{
|
||||
\item{iniName}{The name of the ini file}
|
||||
|
||||
\item{depTree}{The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]}
|
||||
}
|
||||
\description{
|
||||
This function checks the MuSo file system, if it is correct
|
||||
}
|
||||
@@ -0,0 +1,21 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/multiSite.R
|
||||
\name{compareCalibratedWithOriginal}
|
||||
\alias{compareCalibratedWithOriginal}
|
||||
\title{compareCalibratedWithOriginal}
|
||||
\usage{
|
||||
compareCalibratedWithOriginal(
|
||||
key,
|
||||
modOld,
|
||||
modNew,
|
||||
mes,
|
||||
likelihoods,
|
||||
alignIndexes,
|
||||
musoCodeToIndex,
|
||||
nameGroupTable,
|
||||
groupFun
|
||||
)
|
||||
}
|
||||
\description{
|
||||
This functions compareses the likelihood and the RMSE values of the simulations and the measurements
|
||||
}
|
||||
@@ -0,0 +1,25 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/flat.R
|
||||
\name{flatMuso}
|
||||
\alias{flatMuso}
|
||||
\title{flatMuso}
|
||||
\usage{
|
||||
flatMuso(
|
||||
iniName,
|
||||
execPath = "./",
|
||||
depTree = options("RMuso_depTree")[[1]],
|
||||
directory = "flatdir",
|
||||
d = TRUE,
|
||||
outE = TRUE
|
||||
)
|
||||
}
|
||||
\arguments{
|
||||
\item{iniName}{The name of the ini file}
|
||||
|
||||
\item{depTree}{The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]}
|
||||
|
||||
\item{directory}{The destination directory for flattening. At default it will be flatdir}
|
||||
}
|
||||
\description{
|
||||
This function reads the ini file and creates a directory (named after the directory argument) with all the files the modell uses with this file. the directory will be flat.
|
||||
}
|
||||
@@ -0,0 +1,23 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/flat.R
|
||||
\name{getFilePath}
|
||||
\alias{getFilePath}
|
||||
\title{getFilePath}
|
||||
\usage{
|
||||
getFilePath(
|
||||
iniName,
|
||||
fileType,
|
||||
execPath = "./",
|
||||
depTree = options("RMuso_depTree")[[1]]
|
||||
)
|
||||
}
|
||||
\arguments{
|
||||
\item{iniName}{The name of the ini file}
|
||||
|
||||
\item{depTree}{The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]}
|
||||
|
||||
\item{filetype}{The type of the choosen file. For options see options("RMuso_depTree")[[1]]$name}
|
||||
}
|
||||
\description{
|
||||
This function reads the ini file and for a chosen fileType it gives you the filePath
|
||||
}
|
||||
@@ -0,0 +1,20 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/flat.R
|
||||
\name{getFilesFromIni}
|
||||
\alias{getFilesFromIni}
|
||||
\title{getFilesFromIni}
|
||||
\usage{
|
||||
getFilesFromIni(
|
||||
iniName,
|
||||
execPath = "./",
|
||||
depTree = options("RMuso_depTree")[[1]]
|
||||
)
|
||||
}
|
||||
\arguments{
|
||||
\item{iniName}{The name of the ini file}
|
||||
|
||||
\item{depTree}{The file dependency defining dataframe. At default it is: options("RMuso_depTree")[[1]]}
|
||||
}
|
||||
\description{
|
||||
This function reads the ini file and gives yout back the path of all file involved in model run
|
||||
}
|
||||
@@ -0,0 +1,64 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/multiSite.R
|
||||
\name{multiSiteCalib}
|
||||
\alias{multiSiteCalib}
|
||||
\title{multiSiteCalib}
|
||||
\usage{
|
||||
multiSiteCalib(
|
||||
measurements,
|
||||
calTable,
|
||||
parameters,
|
||||
dataVar,
|
||||
iterations = 100,
|
||||
burnin = ifelse(iterations < 3000, 3000, NULL),
|
||||
likelihood,
|
||||
execPath,
|
||||
thread_prefix = "thread",
|
||||
numCores = (parallel::detectCores() - 1),
|
||||
pb = txtProgressBar(min = 0, max = iterations, style = 3),
|
||||
pbUpdate = setTxtProgressBar,
|
||||
copyThread = TRUE,
|
||||
constraints = NULL,
|
||||
th = 10,
|
||||
treeControl = rpart.control()
|
||||
)
|
||||
}
|
||||
\arguments{
|
||||
\item{calTable}{A dataframe which contantains the ini file locations and the domains they belongs to}
|
||||
|
||||
\item{parameters}{A dataframe with the name, the minimum, and the maximum value for the parameters used in MonteCarlo experiment}
|
||||
|
||||
\item{dataVar}{A named vector where the elements are the MuSo variable codes and the names are the same as provided in measurements and likelihood}
|
||||
|
||||
\item{iterations}{The number of MonteCarlo experiments to be executed}
|
||||
|
||||
\item{burnin}{Currently not used, altought it is the length of burnin period of the MCMC sampling used to generate random parameters}
|
||||
|
||||
\item{likelihood}{A list of likelihood functions which names are linked to dataVar}
|
||||
|
||||
\item{execPath}{If you are running the calibration from different location than the MuSo executable, you have to provide the path}
|
||||
|
||||
\item{thread_prefix}{The prefix of thread directory names in the tmp directory created during the calibrational process}
|
||||
|
||||
\item{numCores}{The number of processes used during the calibration. At default it uses one less than the number of threads available}
|
||||
|
||||
\item{pb}{The progress bar function. If you use (web-)GUI you can provide a different function}
|
||||
|
||||
\item{pbUpdate}{The update function for pb (progress bar)}
|
||||
|
||||
\item{copyThread}{A boolean, recreate tmp directory for calibration or not (case of repeating the calibration)}
|
||||
|
||||
\item{th}{A trashold value for multisite calibration. What percentage of the site should satisfy the constraints.}
|
||||
|
||||
\item{treeControl}{A list which controls (maximal complexity, maximal depth) the details of the decession tree making.}
|
||||
|
||||
\item{measuremets}{The table which contains the measurements}
|
||||
|
||||
\item{contsraints}{A dataframe containing the constraints logic the minimum and a maximum value for the calibration.}
|
||||
}
|
||||
\description{
|
||||
This funtion uses the Monte Carlo technique to uniformly sample the parameter space from user defined parameters of the Biome-BGCMuSo model. The sampling algorithm ensures that the parameters are constrained by the model logic which means that parameter dependencies are fully taken into account (parameter dependency means that e.g leaf C:N ratio must be smaller than C:N ratio of litter; more complicated rules apply to the allocation parameters where the allocation fractions to different plant compartments must sum up 1). This function implements a mathematically correct solution to provide uniform distriution of the random parameters on convex polytopes.
|
||||
}
|
||||
\author{
|
||||
Roland HOLLOS
|
||||
}
|
||||
@@ -0,0 +1,36 @@
|
||||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/multiSite.R
|
||||
\name{multiSiteThread}
|
||||
\alias{multiSiteThread}
|
||||
\title{multiSiteThread}
|
||||
\usage{
|
||||
multiSiteThread(
|
||||
measuredData,
|
||||
parameters = NULL,
|
||||
startDate = NULL,
|
||||
endDate = NULL,
|
||||
formatString = "\%Y-\%m-\%d",
|
||||
calTable,
|
||||
dataVar,
|
||||
outLoc = "./calib",
|
||||
outVars = NULL,
|
||||
iterations = 300,
|
||||
skipSpinup = TRUE,
|
||||
plotName = "calib.jpg",
|
||||
modifyOriginal = TRUE,
|
||||
likelihood,
|
||||
uncertainity = NULL,
|
||||
burnin = NULL,
|
||||
naVal = NULL,
|
||||
postProcString = NULL,
|
||||
threadNumber,
|
||||
constraints = NULL,
|
||||
th = 10
|
||||
)
|
||||
}
|
||||
\description{
|
||||
This is an
|
||||
}
|
||||
\author{
|
||||
Roland HOLLOS
|
||||
}
|
||||
@@ -4,7 +4,13 @@
|
||||
\alias{musoRand}
|
||||
\title{musoRand}
|
||||
\usage{
|
||||
musoRand(parameters, iterations = 3000, fileType = "epc", constrains = NULL)
|
||||
musoRand(
|
||||
parameters,
|
||||
iterations = 3000,
|
||||
fileType = "epc",
|
||||
constrains = NULL,
|
||||
burnin = NULL
|
||||
)
|
||||
}
|
||||
\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.}
|
||||
|
||||
@@ -10,7 +10,7 @@ saveAllMusoPlots(
|
||||
silent = TRUE,
|
||||
type = "line",
|
||||
outFile = "annual.csv",
|
||||
colour = NULL,
|
||||
colour = "blue",
|
||||
skipSpinup = FALSE
|
||||
)
|
||||
}
|
||||
|
||||
@@ -13,7 +13,7 @@ updateMusoMapping(excelName, dest = "./", version = getOption("RMuso_version"))
|
||||
The output code-variable matrix, and also the function changes the global variable
|
||||
}
|
||||
\description{
|
||||
This function updates the Biome-BGCMuSo output code-variable matrix. Within Biome-BGCMuSo the state variables and fluxes are marked by integer numbers. In order to provide meaningful variable names (e.g. 3009 means Gross Primary Production in Biome-BGCMuSo v5) a conversion table is needed which is handled by this function.
|
||||
This function updates the Biome-BGCMuSo output code-variable matrix (creates a json file that is used internally by RBBGCMuso). Within Biome-BGCMuSo the output state variablesare marked by integer numbers (see the User's Guide). In order to provide meaningful variable names (e.g. 3009 means Gross Primary Production) a conversion table is needed which is handled by this function. The input Excel file must have the following column order: name, index, units, description (plus other optional columns line group). name refers to the abbreviation of the variable; index is the integer number of the output variable; unit is the unit of the variable; description is a meaningful text to explain the variable. The script will NOT work with other column order!
|
||||
}
|
||||
\author{
|
||||
Roland HOLLOS
|
||||
|
||||
Reference in New Issue
Block a user