upload binary package
This commit is contained in:
hollorol
parent
cb403ed7ae
commit
52b0f0553b
@ -1,4 +0,0 @@
|
||||
musoFilter <- function(dB=.,text){
|
||||
eval(parse(text = paste0("filter(",dB,",.,",text,")"))) %>%
|
||||
tbl_df
|
||||
}
|
||||
@ -1,172 +0,0 @@
|
||||
#' musoRand
|
||||
#'
|
||||
#' This funtion samples uniformly from the choosen parameters of the BiomeBGC-Muso model, which parameters are constrained by the model logic.
|
||||
#' @author Roland Hollos
|
||||
#' @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 input-file, the last two column consist the endpont of the parameter-ranges, where the parameters will be randomized.
|
||||
#' @param constrains This is a matrics wich specify the constrain rules for the sampling. Further informations coming son.
|
||||
#' @param iteration The number of sample-s. It is adviced to use at least 3000 iteration, because it is generally fast and it can be subsampled later at any time.
|
||||
#' @export
|
||||
|
||||
musoRand <- function(parameters, constrains = NULL, iterations=3000){
|
||||
|
||||
if(!is.null(constrains)){
|
||||
constMatrix <- constrains
|
||||
}
|
||||
|
||||
parameters <- parameters[,-1]
|
||||
constMatrix <- constMatrix[,-1]
|
||||
|
||||
depTableMaker <- function(constMatrix,parameters){
|
||||
parameters <- parameters[order(parameters[,1]),]
|
||||
constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],c(5,6)]<-parameters[,c(2,3)]
|
||||
logiConstrain <- (constMatrix[,"GROUP"] %in% constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],"GROUP"] &
|
||||
(constMatrix[,"GROUP"]!=0)) | ((constMatrix[,"INDEX"] %in% parameters[,1]) & (constMatrix[,"GROUP"] == 0))
|
||||
constMatrix<-constMatrix[logiConstrain,]
|
||||
constMatrix <- constMatrix[order(apply(constMatrix[,7:8],1,function(x){x[1]/10+abs(x[2])})),]
|
||||
constMatrix
|
||||
}
|
||||
|
||||
genMat0 <- function(dep){
|
||||
numberOfVariable <- nrow(dep)
|
||||
G <- rbind(diag(numberOfVariable), -1*diag(numberOfVariable))
|
||||
h <- c(dependences[,5], -1*dependences[,6])
|
||||
return(list(G=G,h=h))
|
||||
}
|
||||
|
||||
genMat1 <- function(dep, N){
|
||||
|
||||
## Range <- sapply(list(min,max),function(x){
|
||||
|
||||
## x(as.numeric(rownames(dep)))
|
||||
## }) It is more elegant, more general, but slower
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
numberOfVariables <- nrow(dep)
|
||||
G<- -1*diag(numberOfVariables)
|
||||
|
||||
for(i in 1:numberOfVariables){
|
||||
if(dep[i,4]!=0){
|
||||
G[i,dep[i,4]] <- 1
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
G<-G[dep[,4]!=0,]
|
||||
|
||||
if(Range[1]==1){
|
||||
G<-cbind(G,matrix(ncol=(N-Range[2]),nrow=nrow(G),data=0))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
G<-cbind(matrix(ncol=(Range[1]-1),nrow=nrow(G),data=0),G)
|
||||
} else {
|
||||
G <- cbind(matrix(ncol=(Range[1]-1),nrow=nrow(G),data=0),G,matrix(ncol=(N-Range[2]),nrow=nrow(G),data=0))
|
||||
}
|
||||
}
|
||||
return(list(G=G,h=rep(0,nrow(G))))
|
||||
}
|
||||
|
||||
genMat2 <- function(dep, N){
|
||||
G <- rep(1,nrow(dep))
|
||||
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
if(Range[1]==1){
|
||||
G<-c(G, numeric(N-Range[2]))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
G<-c(numeric(Range[1]-1), G)
|
||||
} else {
|
||||
G <- c(numeric(Range[1]-1), G, numeric(N-Range[2]))
|
||||
}
|
||||
}
|
||||
|
||||
G <- t(matrix(sign(dep[2,4])*G))
|
||||
h <- abs(dep[1,4])
|
||||
|
||||
return(list(G=G,h=h))
|
||||
}
|
||||
|
||||
genMat3 <- function(dep, N){
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
E <- rep(1,nrow(dep))
|
||||
|
||||
if(Range[1]==1){
|
||||
E<-c(E, numeric(N-Range[2]))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
E<-c(numeric(Range[1]-1), E)
|
||||
} else {
|
||||
E <- c(numeric(Range[1]-1), E, numeric(N-Range[2]))
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
E <- t(matrix(E))
|
||||
f <- dep[1,4]
|
||||
return(list(E=E,f=f))
|
||||
}
|
||||
|
||||
|
||||
applyRandTypeG <- function(dep,N){
|
||||
type <- unique(dep[,"TYPE"])
|
||||
minR <- min(dep[,"rowIndex"])
|
||||
maxR <- max(dep[,"rowIndex"])
|
||||
switch(type,
|
||||
invisible(Gh <- genMat1(dep, N)),
|
||||
invisible(Gh <- genMat2(dep, N)))
|
||||
return(Gh)
|
||||
}
|
||||
|
||||
applyRandTypeE <- function(dep,N){
|
||||
type <- unique(dep[,"TYPE"])
|
||||
minR <- min(dep[,"rowIndex"])
|
||||
maxR <- max(dep[,"rowIndex"])
|
||||
switch(-type,
|
||||
stop("Not implemented yet"),
|
||||
stop("Not implemented yet"),
|
||||
invisible(Ef <- genMat3(dep, N)))
|
||||
return(Ef)
|
||||
}
|
||||
|
||||
dependences <- depTableMaker(constMatrix, parameters)
|
||||
dependences <- cbind(dependences,1:nrow(dependences))
|
||||
colnames(dependences)[ncol(dependences)] <- "rowIndex"
|
||||
numberOfVariable <- nrow(dependences)
|
||||
nonZeroDeps<-dependences[dependences[,"TYPE"]!=0,]
|
||||
if(nrow(nonZeroDeps)!=0){
|
||||
splitedDeps<- split(nonZeroDeps,nonZeroDeps[,"GROUP"])
|
||||
Gh <- list()
|
||||
Ef <- list()
|
||||
|
||||
for(i in 1:length(splitedDeps)){
|
||||
print(splitedDeps[[i]][1,"TYPE"])
|
||||
if(splitedDeps[[i]][1,"TYPE"]>0){
|
||||
Gh[[i]]<-applyRandTypeG(splitedDeps[[i]],nrow(dependences))
|
||||
} else {
|
||||
Ef[[i]] <- applyRandTypeE(splitedDeps[[i]],nrow(dependences))
|
||||
}
|
||||
}
|
||||
|
||||
Gh0<- genMat0(dependences)
|
||||
G <- do.call(rbind,lapply(Gh,function(x){x$G}))
|
||||
G<- rbind(Gh0$G,G)
|
||||
h <- do.call(c,lapply(Gh,function(x){x$h}))
|
||||
h <- c(Gh0$h,h)
|
||||
E <- do.call(rbind,lapply(Ef,function(x){x$E}))
|
||||
f <- do.call(c,lapply(Ef,function(x){x$f}))
|
||||
randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,iter = iterations))$X
|
||||
} else{
|
||||
Gh0<-genMat0(dependences)
|
||||
randVal <- suppressWarnings(limSolve::xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X
|
||||
}
|
||||
|
||||
results <- list(INDEX =dependences$INDEX, randVal=randVal)
|
||||
return(results)
|
||||
}
|
||||
@ -1,348 +0,0 @@
|
||||
#' calibMuso
|
||||
#'
|
||||
#' This function changes the epc file and after that runs the BBGC-MuSo model and reads in its outputfile in a well-structured way.
|
||||
#'
|
||||
#' @author Roland Hollos
|
||||
#' @param settings You have to run the setupMuso function before calibMuso. It is its output which contains all of the necessary system variables. It sets the whole running environment
|
||||
#' @param timee The required timesteps in the modell output. It can be "d", if it is daily, "m", if it's monthly, "y", it it is yearly. I recommend to use daily data, the yearly and monthly data is not well-tested yet.
|
||||
#' @param debugging If it is TRUE, it copies the log file to a Log directory to store it, if it is stamplog it contatenate a number before the logfile, which is one more than the maximum of the represented ones in the LOG directory. If it is true or stamplog it collects the "wrong" logfiles
|
||||
#' @param keepEpc If TRUE, it keeps the epc file and stamp it, after these copies it to the EPCS directory. If debugging True or false, it copies the wrong epc files to the wrong epc directory.
|
||||
#' @param export if it is yes or you give a filename here, it converts the ouxtput to the specific extension. For example, if you set export to "example.csv", it converts the output to "csv", if you set it to "example.xls" it converts to example.xls with the xlsx package. If it is not installed it gives back a warning message and converts it to csv.
|
||||
#' @param silent If you set it TRUE all off the modells output to the screen will be suppressed. It can be usefull, because it increases the model-speed.
|
||||
#' @param aggressive It deletes every possible modell-outputs from the previous modell runs.
|
||||
#' @param parameters In the settings variable you have set the row indexes of the variables, you wish to change. In this parameter you can give an exact value for them in a vector like: c(1,2,3,4)
|
||||
#' @param logfilename If you want to set a specific name for your logfiles you can set this via logfile parameter
|
||||
#' @param leapYear Should the function do a leapyear correction on the outputdata? If TRUE, then the 31.12 day will be doubled.
|
||||
#' @param keepBinary In default RBBGCMuso to keep working area as clean as possible, deletes all the regular output files. The results are directly printed to the standard output, but you can redirect it, and save it to a variable, or you can export your results to the desired destination in a desired format. Whith this variable you can enable to keep the binary output files. If you want to set the location of the binary output, please take a look at the binaryPlace argument.
|
||||
#' @param binaryPlace The place of the binary output files.
|
||||
#' @param fileToChange You can change any line of the epc or the ini file, you just have to specify with this variable which file you van a change. Two options possible: "epc", "ini"
|
||||
#' @param skipSpinup If TRUE, calibMuso wont do spinup simulation
|
||||
#' @return No return, outputs are written to file
|
||||
#' @usage calibMuso(settings,parameters=NULL, timee="d", debugging=FALSE, logfilename=NULL,
|
||||
#' keepEpc=FALSE, export=FALSE, silent=FALSE, aggressive=FALSE, leapYear=FALSE)
|
||||
#' @import utils
|
||||
#' @export
|
||||
|
||||
calibMuso <- function(settings=NULL, calibrationPar=NULL,
|
||||
parameters=NULL, timee="d",
|
||||
debugging=FALSE, logfilename=NULL,
|
||||
keepEpc=FALSE, export=FALSE,
|
||||
silent=FALSE, aggressive=FALSE,
|
||||
leapYear=FALSE,keepBinary=FALSE,
|
||||
binaryPlace="./", fileToChange="epc",
|
||||
skipSpinup = TRUE, modifyOriginal =FALSE){
|
||||
##########################################################################
|
||||
###########################Set local variables and places########################
|
||||
########################################################################
|
||||
if(is.null(settings)){
|
||||
settings <- setupMuso()
|
||||
}
|
||||
|
||||
Linuxp <-(Sys.info()[1]=="Linux")
|
||||
##Copy the variables from settings
|
||||
inputLoc <- settings$inputLoc
|
||||
outputLoc <- settings$outputLoc
|
||||
outputNames <- settings$outputNames
|
||||
executable <- settings$executable
|
||||
iniInput <- settings$iniInput
|
||||
epc <- settings$epcInput
|
||||
calibrationPar <- settings$calibrationPar
|
||||
binaryPlace <- normalizePath(binaryPlace)
|
||||
whereAmI<-getwd()
|
||||
|
||||
|
||||
## Set the working directory to the inputLoc temporarly.
|
||||
setwd(inputLoc)
|
||||
|
||||
|
||||
if(debugging){#If debugging option turned on
|
||||
#If log or ERROR directory does not exists create it!
|
||||
dirName<-file.path(inputLoc,"LOG")
|
||||
dirERROR<-file.path(inputLoc,"ERROR")
|
||||
|
||||
if(!dir.exists(dirName)){
|
||||
dir.create(dirName)
|
||||
}
|
||||
|
||||
if(!dir.exists(dirERROR)){
|
||||
dir.create(dirERROR)
|
||||
}
|
||||
}
|
||||
|
||||
if(keepEpc) {
|
||||
epcdir <- dirname(epc[1])
|
||||
print(epcdir)
|
||||
WRONGEPC<-file.path(inputLoc,"WRONGEPC")
|
||||
EPCS<-file.path(inputLoc,"EPCS")
|
||||
|
||||
if(!dir.exists(WRONGEPC)){
|
||||
dir.create(WRONGEPC)
|
||||
}
|
||||
|
||||
if(!dir.exists(EPCS)){
|
||||
dir.create(EPCS)
|
||||
}
|
||||
}
|
||||
|
||||
#############################################################
|
||||
############################spinup run############################
|
||||
##########################################################
|
||||
|
||||
|
||||
|
||||
|
||||
if(aggressive==TRUE){
|
||||
cleanupMuso(location=outputLoc,deep = TRUE)
|
||||
}
|
||||
|
||||
toModif<-c(epc[2],iniInput[2])
|
||||
|
||||
if(!modifyOriginal & !is.null(parameters))
|
||||
{
|
||||
|
||||
toModif <- sapply(toModif, function (x){
|
||||
paste0(tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))
|
||||
})
|
||||
|
||||
}
|
||||
|
||||
##change the epc file if and only if there are given parameters
|
||||
if(!is.null(parameters)){
|
||||
changemulline(filePaths=c(epc[2],iniInput[2]), calibrationPar = calibrationPar,
|
||||
contents = parameters, fileOut=toModif, fileToChange=fileToChange, modifyOriginal=modifyOriginal)
|
||||
}
|
||||
|
||||
|
||||
##We change the working directory becase of the model, but we want to avoid sideeffects, so we save the current location and after that we will change everything to it.
|
||||
if(!modifyOriginal & !is.null(parameters)){
|
||||
epc[2]<-file.path(dirname(epc[2]),toModif[1]) # Writing back the lost path
|
||||
toModif[2]<-file.path(dirname(iniInput[2]),toModif[2]) #for the Initmp, also
|
||||
if(!file.exists(toModif[2]) & !modifyOriginal){
|
||||
file.copy(iniInput[2],toModif[2],overwrite = TRUE)
|
||||
}
|
||||
|
||||
iniInput[2] <- toModif[2]}
|
||||
|
||||
if(!is.null(parameters) & ((fileToChange == "epc") | (fileToChange == "both")) & !modifyOriginal){
|
||||
tmp<-readLines(iniInput[2])
|
||||
tmpInd<-grep("EPC_FILE",tmp)+1
|
||||
tmp[tmpInd]<-file.path(dirname(tmp[tmpInd]),basename(epc[2]))
|
||||
writeLines(tmp,iniInput[2])
|
||||
rm(list=c("tmp","tmpInd"))
|
||||
}
|
||||
|
||||
if(!skipSpinup) {
|
||||
|
||||
##Run the model for the spinup run.
|
||||
|
||||
if(silent){#silenc mode
|
||||
if(Linuxp){
|
||||
#In this case, in linux machines
|
||||
tryCatch(system(paste(executable,iniInput[1],"> /dev/null",sep=" ")),
|
||||
error= function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot run the modell-check the executable!")})
|
||||
} else {
|
||||
#In windows machines there is a show.output.on.console option
|
||||
tryCatch(system(paste(executable,iniInput[1],sep=" "),show.output.on.console = FALSE),
|
||||
error= function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot run the modell-check the executable!")})
|
||||
}
|
||||
|
||||
} else {
|
||||
system(paste(executable,iniInput[1],sep=" "))
|
||||
}
|
||||
|
||||
|
||||
logspinup <- getLogs(outputLoc,outputNames,type="spinup")
|
||||
## logspinup <- grep(paste0(outputNames[1],".log"), list.files(outputLoc),value = TRUE)
|
||||
## logspinup <- list.files(outputLoc)[grep("log$",list.files(outputLoc))]#load the logfiles
|
||||
|
||||
if(length(logspinup)==0){
|
||||
if(keepEpc){
|
||||
stampnum<-stamp(EPCS)
|
||||
lapply(epc,function (x) file.copy(from = x ,to=paste(EPCS,"/",(stampnum+1),"-", basename(x),sep="")))
|
||||
lapply(epc, function (x) file.copy(from = paste(EPCS,"/",(stampnum+1),"-",basename(x),sep=""), to=WRONGEPC))
|
||||
setwd(whereAmI)
|
||||
stop("Modell Failure")
|
||||
}
|
||||
setwd(whereAmI)
|
||||
stop("Modell Failure") #in that case the modell did not create even a logfile
|
||||
}
|
||||
|
||||
if(length(logspinup)>1){
|
||||
spincrash<-TRUE
|
||||
} else {
|
||||
if(identical(tail(readLines(paste(outputLoc,logspinup,sep="/"),-1),1),character(0))){
|
||||
spincrash<-TRUE
|
||||
} else {
|
||||
spincrash <- (tail(readLines(paste(outputLoc,logspinup,sep="/"),-1),1)!=1)
|
||||
}
|
||||
}
|
||||
} else {spincrash <- FALSE}
|
||||
#If the last line in the logfile is 0 There are mistakes so the spinup crashes
|
||||
|
||||
if(!spincrash){##If spinup did not crashed, run the normal run.
|
||||
|
||||
#####################################################################
|
||||
###########################normal run#########################
|
||||
#################################################################
|
||||
|
||||
##for the sake of safe we set the location again
|
||||
setwd(inputLoc)
|
||||
|
||||
if(silent){
|
||||
if(Linuxp){
|
||||
tryCatch(system(paste(executable,iniInput[2],"> /dev/null",sep=" ")),
|
||||
error =function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot run the modell-check the executable!")})
|
||||
} else {
|
||||
tryCatch(system(paste(executable,iniInput[2],sep=" "),show.output.on.console = FALSE),
|
||||
error =function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot run the modell-check the executable!")} )
|
||||
}
|
||||
|
||||
} else {
|
||||
tryCatch(system(paste(executable,iniInput[2],sep=" ")),
|
||||
error =function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot run the modell-check the executable!")})
|
||||
}
|
||||
|
||||
|
||||
##read the output
|
||||
|
||||
switch(timee,
|
||||
"d"=(Reva <- tryCatch(getdailyout(settings), #(:INSIDE: getOutput.R )
|
||||
error = function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot read binary output, please check if the output type is set 2 in the ini files!")})),
|
||||
"m"=(Reva <- tryCatch(getmonthlyout(settings), #(:INSIDE: getOutput.R )
|
||||
error = function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot read binary output, please check if the output type is set 2 in the ini files!")})),
|
||||
"y"=(Reva <- tryCatch(getyearlyout(settings), #(:INSIDE: getOutput.R )
|
||||
error = function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot read binary output, please check if the output type is set 2 in the ini files!")}))
|
||||
)
|
||||
|
||||
if(keepBinary){
|
||||
possibleNames <- tryCatch(getOutFiles(outputLoc = outputLoc,outputNames = outputNames),
|
||||
error=function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot find output files")})
|
||||
stampAndDir(outputLoc = outputLoc,names = possibleNames,stampDir=binaryPlace,type="output")
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if(skipSpinup){
|
||||
logfiles <- tryCatch(getLogs(outputLoc,outputNames,type="normal"),
|
||||
error = function (e){
|
||||
setwd(whereAmI)
|
||||
stop("Cannot find log files, something went wrong")})
|
||||
} else {
|
||||
logfiles <- tryCatch(getLogs(outputLoc,outputNames,type="both"),
|
||||
error = function (e){
|
||||
setwd(whereAmI)
|
||||
stop("Cannot find log files, something went wrong")})
|
||||
}
|
||||
## list.files(outputLoc)[grep("log$",list.files(outputLoc))]#creating a vector for logfilenames
|
||||
|
||||
###############################################
|
||||
#############LOG SECTION#######################
|
||||
###############################################
|
||||
|
||||
if(skipSpinup){
|
||||
errorsign <- readErrors(outputLoc=outputLoc,logfiles=logfiles,type="normal")
|
||||
} else {
|
||||
|
||||
perror <- readErrors(outputLoc=outputLoc,logfiles=logfiles) #vector of spinup and normalrun error
|
||||
|
||||
|
||||
##if errorsign is 1 there is error, if it is 0 everything ok
|
||||
perror[is.na(perror)]<-0
|
||||
if(length(perror)>sum(perror)){
|
||||
errorsign <- 1
|
||||
} else {
|
||||
if(length(perror)==1){
|
||||
errorsign <- 1
|
||||
} else {
|
||||
if(spincrash){
|
||||
errorsign <- 1
|
||||
} else {
|
||||
errorsign <- 0
|
||||
} }
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
if(keepEpc){#if keepepc option turned on
|
||||
|
||||
if(length(unique(dirname(epc)))>1){
|
||||
stop("Why are you playing with my nervs? Seriously? You hold your epc-s in different folders?")
|
||||
} else {
|
||||
if(skipSpinup){
|
||||
stampAndDir(stampDir=EPCS, wrongDir=WRONGEPC, names=epc[2], type="general", errorsign=errorsign, logfiles=logfiles)
|
||||
}
|
||||
stampAndDir(stampDir=EPCS, wrongDir=WRONGEPC, names=epc, type="general", errorsign=errorsign, logfiles=logfiles)
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(debugging){ #debugging is boolean
|
||||
logfiles <- file.path(outputLoc,logfiles)
|
||||
|
||||
stampAndDir(stampDir=dirName, wrongDir=dirERROR, names=logfiles, type="general",errorsign=errorsign,logfiles=logfiles)}
|
||||
|
||||
|
||||
#cleanupMuso(location=outputLoc,deep = FALSE)
|
||||
if(errorsign==1){
|
||||
return("Modell Failure")
|
||||
}
|
||||
|
||||
if(timee=="d"){
|
||||
colnames(Reva) <- unlist(settings$outputVars[[1]])
|
||||
} else {
|
||||
if(timee=="y")
|
||||
colnames(Reva) <- unlist(settings$outputVars[[2]])
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(leapYear){
|
||||
Reva <- corrigMuso(settings,Reva)
|
||||
rownames(Reva) <- musoDate(settings$startYear,settings$numYears)
|
||||
} else {
|
||||
rownames(Reva) <- musoDate(settings$startYear, settings$numYears, corrigated=FALSE)
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(export!=FALSE){
|
||||
setwd(whereAmI)
|
||||
|
||||
## switch(fextension(export),
|
||||
## "csv"=(write.csv(Reva,export)),
|
||||
## "xlsx"=(),
|
||||
## "odt"=
|
||||
|
||||
|
||||
## )
|
||||
write.csv(Reva,export)
|
||||
|
||||
} else{
|
||||
setwd(whereAmI)
|
||||
return(Reva)}
|
||||
}
|
||||
@ -1,50 +0,0 @@
|
||||
#' This is the function which is capable change multiple specific lines to other using their row numbers.
|
||||
#'
|
||||
#' he function uses the previous changspecline function to operate.
|
||||
##From now changespecline is in the forarcheologist file, because its no longer needed
|
||||
#'
|
||||
#' @author Roland Hollos
|
||||
#' @keywords internal
|
||||
#'
|
||||
|
||||
changemulline <- function(filePaths, calibrationPar, contents, fileOut=NULL, fileToChange, modifyOriginal=FALSE){
|
||||
selectFileToWrite <- function(filePaths, fileTochange){
|
||||
if(fileToChange == "epc"){
|
||||
return(filePaths[1])
|
||||
} else{
|
||||
return(filePaths[2])
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if(modifyOriginal){
|
||||
fileOut <- filePaths
|
||||
mutated <- TRUE
|
||||
} else {
|
||||
fileOut <- sapply(filePaths, function (x){
|
||||
paste0(tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(x))
|
||||
})
|
||||
mutated <- FALSE
|
||||
}
|
||||
|
||||
if(xor(is.list(calibrationPar), is.list(contents)) ){
|
||||
stop("If you change epc and ini files also, you have to use list for calibrationPar, and paramateters.")
|
||||
}
|
||||
|
||||
if(!is.element(fileToChange,c("ini","epc","both"))){
|
||||
stop("RBBGCMuso can only change ini or epc file, so fileToChange can be 'epc/ini/both'")
|
||||
}
|
||||
|
||||
if(fileToChange == "epc" | fileToChange == "ini"){
|
||||
changeMusoC(selectFileToWrite(filePaths, fileToChange),fileOut[1], cbind(calibrationPar, contents))
|
||||
}
|
||||
|
||||
if(fileToChange == "both"){
|
||||
changeMusoC(filePaths[1],fileOut[1], cbind(calibrationPar[[1]],contents[[1]]))
|
||||
changeMusoC(filePaths[2],fileOut[2], cbind(calibrationPar[[1]],contents[[2]]))
|
||||
}
|
||||
|
||||
return(mutated)
|
||||
}
|
||||
|
||||
|
||||
@ -1,200 +0,0 @@
|
||||
#' musoMont
|
||||
#'
|
||||
#' This function does monteCarlo on BiomeBGC-MuSo. It samples specified modell variables in given rangge from conditional multivariate uniform distribution, and runs the modell for each run.
|
||||
#' @author Roland Hollos
|
||||
#' @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 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 inputDir The location of the input directory, this directory must content a viable pack of all inputfiles and the executable file.
|
||||
#' @param iterations number of the monteCarlo run.
|
||||
#' @param preTag It will be the name of the output files. For example preTag-1.csv, pretag-2csv...
|
||||
#' @param outputType This parameter can be "oneCsv", "moreCsv", and "netCDF". If "oneCsv" is choosen the function create 1 big csv file for all of the runs, if "moreCsv" is choosen, every modell output goes to separate files, if netCDF is selected the outputs will be put in a netCDF file. The default value of the outputTypes is "moreCsv". netCDF is not implemented yet.
|
||||
#' @param fun If you select a variable from the possible outputs (with specify the varIndex parameter), you have to provide a function which maps to a subset of real numbers. The most frequent possibilities are: mean, min, max, var, but you can define any function for your need.
|
||||
#' @param varIndex This parameter specify which parameter of the output will be used. You can extract this information from the ini-files. At the output parameter specifications, the parameters order will determine this number. For example, if you have set these output parameters: 412, 874, 926, 888, and you want to use 926, you should address varIndex with 3.
|
||||
#' @param debugging If you set this parameter, you can save every logfile, and RBBGCMuso will select those which contains errors.
|
||||
#' @param keepEpc if you set keepEpc also true, it will save every selected epc file, and put the wrong ones in the WRONGEPC directory.
|
||||
#' @export
|
||||
|
||||
musoMont <- function(settings=NULL,
|
||||
parameters=NULL,
|
||||
inputDir = "./",
|
||||
outLoc = "./calib",
|
||||
iterations = 10,
|
||||
preTag = "mont-",
|
||||
outputType = "moreCsv",
|
||||
fun=mean,
|
||||
varIndex = 1,
|
||||
outVars = NULL,
|
||||
silent = TRUE,
|
||||
skipSpinup = TRUE,
|
||||
debugging = FALSE,
|
||||
keepEpc = FALSE,
|
||||
constrains = NULL,
|
||||
...){
|
||||
|
||||
|
||||
readValuesFromEpc <- function(epc, linums){
|
||||
rows <- numeric(2)
|
||||
values <- sapply(linums, function(x){
|
||||
rows[1] <- as.integer(x)
|
||||
rows[2] <- as.integer(round(100*x)) %% 10 + 1
|
||||
epcFile <- readLines(epc)
|
||||
|
||||
return(as.numeric(unlist(strsplit(epcFile[rows[1]], split= "[\t ]"))[rows[2]]))
|
||||
|
||||
})
|
||||
|
||||
return(values)
|
||||
}
|
||||
|
||||
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")
|
||||
})
|
||||
}}
|
||||
|
||||
outLocPlain <- basename(outLoc)
|
||||
currDir <- getwd()
|
||||
|
||||
if(!dir.exists(outLoc)){
|
||||
dir.create(outLoc)
|
||||
warning(paste(outLoc," is not exists, so it was created"))
|
||||
}
|
||||
|
||||
outLoc <- normalizePath(outLoc)
|
||||
|
||||
if(is.null(settings)){
|
||||
settings <- setupMuso()
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(is.null(outVars)){
|
||||
numVars <- length(settings$outputVars[[1]])
|
||||
outVarNames <- settings$outputVars[[1]]
|
||||
} else {
|
||||
numVars <- length(outVars)
|
||||
outVarNames <- sapply(outVars, musoMapping)
|
||||
}
|
||||
|
||||
parameterNames <- parameters[,1]
|
||||
# settings$calibrationPar <- A[,1] #:LATER:
|
||||
pretag <- file.path(outLoc,preTag)
|
||||
npar <- length(settings$calibrationPar)
|
||||
|
||||
##reading the original epc file at the specified
|
||||
## row numbers
|
||||
if(iterations < 3000){
|
||||
randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = 3000)
|
||||
randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),]
|
||||
} else {
|
||||
randVals <- musoRand(parameters = parameters,constrains = constrains, iterations = iterations)
|
||||
}
|
||||
|
||||
origEpc <- readValuesFromEpc(settings$epc[2],parameters[,2])
|
||||
|
||||
## Prepare the preservedEpc matrix for the faster
|
||||
## run.
|
||||
|
||||
pretag <- file.path(outLoc,preTag)
|
||||
|
||||
## Creating function for generating separate
|
||||
## csv files for each run
|
||||
|
||||
progBar <- txtProgressBar(1,iterations,style=3)
|
||||
|
||||
|
||||
moreCsv <- function(){
|
||||
randValues <- randVals[[2]]
|
||||
settings$calibrationPar <- randVals[[1]]
|
||||
## randValues <- randValues[,randVals[[1]] %in% parameters[,2]][,rank(parameters[,2])]
|
||||
modellOut <- matrix(ncol = numVars, nrow = iterations + 1)
|
||||
|
||||
origModellOut <- calibMuso(silent=TRUE)
|
||||
write.csv(x=origModellOut, file=paste0(pretag,1,".csv"))
|
||||
|
||||
if(!is.list(fun)){
|
||||
funct <- rep(list(fun), numVars)
|
||||
}
|
||||
|
||||
tmp2 <- numeric(numVars)
|
||||
|
||||
for(j in 1:numVars){
|
||||
tmp2[j]<-funct[[j]](origModellOut[,j])
|
||||
}
|
||||
modellOut[1,]<- tmp2
|
||||
|
||||
for(i in 2:(iterations+1)){
|
||||
tmp <- calibMuso(settings = settings,
|
||||
parameters = randValues[(i-1),],
|
||||
silent= TRUE,
|
||||
skipSpinup = skipSpinup,
|
||||
keepEpc = keepEpc,
|
||||
debugging = debugging,
|
||||
outVars = outVars)
|
||||
|
||||
|
||||
for(j in 1:numVars){
|
||||
tmp2[j]<-funct[[j]](tmp[,j])
|
||||
}
|
||||
|
||||
modellOut[i,]<- tmp2
|
||||
write.csv(x=tmp, file=paste0(pretag,(i+1),".csv"))
|
||||
setTxtProgressBar(progBar,i)
|
||||
}
|
||||
|
||||
|
||||
|
||||
preservedEpc <- cbind(rbind(parameters[,2], randValues[,randVals[[1]] %in% parameters[,2]][,rank(parameters[,2])]),
|
||||
modellOut)
|
||||
colnames(preservedEpc) <- c(parameterNames, sapply(outVarNames, function (x) paste0("mod.", x)))
|
||||
return(preservedEpc)
|
||||
}
|
||||
|
||||
## Creating function for generating one
|
||||
## csv files for each run
|
||||
|
||||
oneCsv <- function () {
|
||||
stop("This function is not implemented yet")
|
||||
## 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()))
|
||||
write.csv(a,"preservedEpc.csv")
|
||||
|
||||
setwd(currDir)
|
||||
return(a)
|
||||
}
|
||||
|
||||
@ -1,123 +0,0 @@
|
||||
musoMonte <- function(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,
|
||||
...){
|
||||
|
||||
currDir <- getwd()
|
||||
tmp <- file.path(outLoc,"tmp/")
|
||||
|
||||
if(!dir.exists(outLoc)){
|
||||
dir.create(outLoc)
|
||||
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)
|
||||
|
||||
|
||||
inputFiles <- file.path(inputDir,grep(basename(outLoc),list.files(inputDir),invert = TRUE,value = TRUE))
|
||||
|
||||
|
||||
for(i in inputFiles){
|
||||
file.copy(i,tmp)
|
||||
}
|
||||
setwd(tmp)
|
||||
|
||||
if(is.null(settings)){
|
||||
settings <- setupMuso()
|
||||
}
|
||||
|
||||
parameterNames <- parameters[,1]
|
||||
parReal <- parameters[,-1]
|
||||
|
||||
|
||||
|
||||
##reading the original epc file at the specified
|
||||
## row numbers
|
||||
|
||||
origEpcFile <- readLines(settings$epcinput[2])
|
||||
origEpc <- unlist( lapply(parameters[,2], 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 = nrow(parameters))
|
||||
preservedEpc[1,] <- origEpc
|
||||
colnames(preservedEpc) <- parameters[,1]
|
||||
|
||||
## Save the backupEpc, while change the settings
|
||||
## variable and set the output.
|
||||
file.copy(settings$epc[2],savedEpc,overwrite = TRUE) # do I need this?
|
||||
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)
|
||||
|
||||
## Creating function for generating separate
|
||||
## csv files for each run
|
||||
moreCsv <- function(){
|
||||
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,i,".csv")
|
||||
calibMuso(settings,debugging = "stamplog",
|
||||
parameters = parVar,export = exportName,
|
||||
keepEpc = TRUE)
|
||||
}
|
||||
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" = (preservedEpc <- oneCsv()),
|
||||
"moreCsv" = (preservedEpc <- moreCsv()),
|
||||
"netCDF" = (preservedEpc <- netCDF()))
|
||||
|
||||
## Change back the epc file to the original
|
||||
|
||||
file.copy(savedEpc,settings$epc[2],overwrite = TRUE)
|
||||
write.csv(preservedEpc,"preservedEpc.csv")
|
||||
}
|
||||
@ -1,165 +0,0 @@
|
||||
#' musoRand
|
||||
#'
|
||||
#' This funtion samples uniformly from the choosen parameters of the BiomeBGC-Muso model, which parameters are constrained by the model logic.
|
||||
#' @author Roland Hollos
|
||||
#' @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 input-file, the last two column consist the endpont of the parameter-ranges, where the parameters will be randomized.
|
||||
#' @param constrains This is a matrics wich specify the constrain rules for the sampling. Further informations coming son.
|
||||
#' @param iteration The number of sample-s. It is adviced to use at least 3000 iteration, because it is generally fast and it can be subsampled later at any time.
|
||||
#' @export
|
||||
|
||||
musoRand <- function(parameters, constrains = NULL, iterations=3000){
|
||||
|
||||
if(!is.null(constrains)){
|
||||
constMatrix <- constrains
|
||||
}
|
||||
|
||||
parameters <- parameters[,-1]
|
||||
constMatrix <- constMatrix[,-1]
|
||||
|
||||
depTableMaker <- function(constMatrix,parameters){
|
||||
parameters <- parameters[order(parameters[,1]),]
|
||||
constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],c(5,6)]<-parameters[,c(2,3)]
|
||||
logiConstrain <- (constMatrix[,"GROUP"] %in% constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],"GROUP"] &
|
||||
(constMatrix[,"GROUP"]!=0)) | ((constMatrix[,"INDEX"] %in% parameters[,1]) & (constMatrix[,"GROUP"] == 0))
|
||||
constMatrix<-constMatrix[logiConstrain,]
|
||||
constMatrix <- constMatrix[order(apply(constMatrix[,7:8],1,function(x){x[1]/10+abs(x[2])})),]
|
||||
constMatrix
|
||||
}
|
||||
|
||||
genMat0 <- function(dep){
|
||||
numberOfVariable <- nrow(dep)
|
||||
G <- rbind(diag(numberOfVariable), -1*diag(numberOfVariable))
|
||||
h <- c(dependences[,5], -1*dependences[,6])
|
||||
return(list(G=G,h=h))
|
||||
}
|
||||
|
||||
genMat1 <- function(dep, N){
|
||||
|
||||
## Range <- sapply(list(min,max),function(x){
|
||||
## x(as.numeric(rownames(dep)))
|
||||
## }) It is more elegant, more general, but slower
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
numberOfVariables <- nrow(dep)
|
||||
G<- -1*diag(numberOfVariables)
|
||||
|
||||
for(i in 1:numberOfVariables){
|
||||
if(dep[i,4]!=0){
|
||||
G[i,dep[i,4]] <- 1
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
G<-G[dep[,4]!=0,]
|
||||
|
||||
if(Range[1]==1){
|
||||
G<-cbind(G,matrix(ncol=(N-Range[2]),nrow=nrow(G),data=0))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
G<-cbind(matrix(ncol=(Range[1]-1),nrow=nrow(G),data=0),G)
|
||||
} else {
|
||||
G <- cbind(matrix(ncol=(Range[1]-1),nrow=nrow(G),data=0),G,matrix(ncol=(N-Range[2]),nrow=nrow(G),data=0))
|
||||
}
|
||||
}
|
||||
return(list(G=G,h=rep(0,nrow(G))))
|
||||
}
|
||||
|
||||
genMat2 <- function(dep, N){
|
||||
G <- rep(1,nrow(dep))
|
||||
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
if(Range[1]==1){
|
||||
G<-c(G, numeric(N-Range[2]))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
G<-c(numeric(Range[1]-1), G)
|
||||
} else {
|
||||
G <- c(numeric(Range[1]-1), G, numeric(N-Range[2]))
|
||||
}
|
||||
}
|
||||
|
||||
G <- t(matrix(sign(dep[2,4])*G))
|
||||
h <- abs(dep[1,4])
|
||||
|
||||
return(list(G=G,h=h))
|
||||
}
|
||||
|
||||
genMat3 <- function(dep, N){
|
||||
Range <- (function(x){
|
||||
c(min(x), max(x))
|
||||
})(as.numeric(dep[,"rowIndex"]))
|
||||
|
||||
E <- rep(1,nrow(dep))
|
||||
|
||||
if(Range[1]==1){
|
||||
E<-c(E, numeric(N-Range[2]))
|
||||
} else{
|
||||
if(Range[2]==N){
|
||||
E<-c(numeric(Range[1]-1), E)
|
||||
} else {
|
||||
E <- c(numeric(Range[1]-1), E, numeric(N-Range[2]))
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
E <- t(matrix(E))
|
||||
f <- dep[1,4]
|
||||
return(list(E=E,f=f))
|
||||
}
|
||||
|
||||
|
||||
applyRandTypeG <- function(dep,N){
|
||||
type <- unique(dep[,"TYPE"])
|
||||
minR <- min(dep[,"rowIndex"])
|
||||
maxR <- max(dep[,"rowIndex"])
|
||||
switch(type,
|
||||
invisible(Gh <- genMat1(dep, N)),
|
||||
invisible(Gh <- genMat2(dep, N)))
|
||||
return(Gh)
|
||||
}
|
||||
|
||||
applyRandTypeE <- function(dep,N){
|
||||
type <- unique(dep[,"TYPE"])
|
||||
minR <- min(dep[,"rowIndex"])
|
||||
maxR <- max(dep[,"rowIndex"])
|
||||
switch(-type,
|
||||
stop("Not implemented yet"),
|
||||
stop("Not implemented yet"),
|
||||
invisible(Ef <- genMat3(dep, N)))
|
||||
return(Ef)
|
||||
}
|
||||
|
||||
dependences <- depTableMaker(constMatrix, parameters)
|
||||
dependences <- cbind(dependences,1:nrow(dependences))
|
||||
colnames(dependences)[ncol(dependences)] <- "rowIndex"
|
||||
numberOfVariable <- nrow(dependences)
|
||||
nonZeroDeps<-dependences[dependences[,"TYPE"]!=0,]
|
||||
splitedDeps<- split(nonZeroDeps,nonZeroDeps[,"GROUP"])
|
||||
Gh <- list()
|
||||
Ef <- list()
|
||||
|
||||
for(i in 1:length(splitedDeps)){
|
||||
print(splitedDeps[[i]][1,"TYPE"])
|
||||
if(splitedDeps[[i]][1,"TYPE"]>0){
|
||||
Gh[[i]]<-applyRandTypeG(splitedDeps[[i]],nrow(dependences))
|
||||
} else {
|
||||
Ef[[i]] <- applyRandTypeE(splitedDeps[[i]],nrow(dependences))
|
||||
}
|
||||
}
|
||||
|
||||
Gh0<- genMat0(dependences)
|
||||
G <- do.call(rbind,lapply(Gh,function(x){x$G}))
|
||||
G<- rbind(Gh0$G,G)
|
||||
h <- do.call(c,lapply(Gh,function(x){x$h}))
|
||||
h <- c(Gh0$h,h)
|
||||
E <- do.call(rbind,lapply(Ef,function(x){x$E}))
|
||||
f <- do.call(c,lapply(Ef,function(x){x$f}))
|
||||
randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,iter = iteration))$X
|
||||
results <- list(INDEX =dependences$INDEX, randVal=randVal)
|
||||
return(results)
|
||||
}
|
||||
@ -1,108 +0,0 @@
|
||||
#' musoSensi
|
||||
#'
|
||||
#' This function does regression based sensitivity analysis based on the output of musoMonte.
|
||||
#' @author Roland Hollos
|
||||
#' @param monteCarloFile If you run musoMonte function previously, you did not have to rerun the monteCarlo, just provide the preservedEpc.csv file with its path. If you do not set this parameter, musoSensi will fun the musoMonte function to get all of the information.
|
||||
#' @param outputFile The filename in which the output of musoSensi function will be saved. It's default value is: "sensitivity.csv"
|
||||
#' @param plotName The name of the output barplot. It's default value is: "sensitivity.jpg"
|
||||
#' @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 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 inputDir The location of the input directory, this directory must content a viable pack of all inputfiles and the executable file.
|
||||
#' @param iterations number of the monteCarlo run.
|
||||
#' @param preTag It will be the name of the output files. For example preTag-1.csv, pretag-2csv...
|
||||
#' @param outputType This parameter can be "oneCsv", "moreCsv", and "netCDF". If "oneCsv" is choosen the function create 1 big csv file for all of the runs, if "moreCsv" is choosen, every modell output goes to separate files, if netCDF is selected the outputs will be put in a netCDF file. The default value of the outputTypes is "moreCsv". netCDF is not implemented yet.
|
||||
#' @param fun If you select a variable from the possible outputs (with specify the varIndex parameter), you have to provide a function which maps to a subset of real numbers. The most frequent possibilities are: mean, min, max, var, but you can define any function for your need.
|
||||
#' @param varIndex This parameter specify which parameter of the output will be used. You can extract this information from the ini-files. At the output parameter specifications, the parameters order will determine this number. For example, if you have set these output parameters: 412, 874, 926, 888, and you want to use 926, you should address varIndex with 3.
|
||||
#' @param skipSpinup With this parameter, you can turn of the spinup phase after the first spinup. I will decrease the time significantly.
|
||||
#' @import dplyr
|
||||
#' @import graphics
|
||||
#' @import grDevices
|
||||
#' @import ggplot2
|
||||
#' @export
|
||||
|
||||
musoSensi <- function(monteCarloFile = NULL,
|
||||
parameters = NULL,
|
||||
settings = NULL,
|
||||
inputDir = "./",
|
||||
outLoc = "./calib",
|
||||
iterations = 30,
|
||||
preTag = "mont-",
|
||||
outputType = "moreCsv",
|
||||
fun = mean,
|
||||
varIndex = 1,
|
||||
outputFile = "sensitivity.csv",
|
||||
plotName = "sensitivity.png",
|
||||
plotTitle = "Sensitivity",
|
||||
skipSpinup = TRUE,
|
||||
dpi=300){
|
||||
|
||||
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")
|
||||
})
|
||||
}}
|
||||
|
||||
doSensi <- function(M){
|
||||
npar <- ncol(M)-1
|
||||
M <- M[which(!is.na(M[,ncol(M)])),]
|
||||
y <- M[,(npar+1)]
|
||||
M <- M[,colnames(M) %in% parameters[,1]]
|
||||
npar <- ncol(M)
|
||||
M <- apply(M[,1:npar],2,function(x){x-mean(x)})
|
||||
varNames<- colnames(M)[1:npar]
|
||||
w <- lm(y~M)$coefficients[-1]
|
||||
Sv <- apply(M,2,var)
|
||||
overalVar <- sum(Sv*w^2)
|
||||
S=numeric(npar)
|
||||
|
||||
for(i in 1:npar){
|
||||
S[i] <- ((w[i]^2*Sv[i])/(overalVar))*100
|
||||
}
|
||||
|
||||
S <- round(S,digits=2)
|
||||
names(S)<-varNames
|
||||
write.csv(file = outputFile, x = S)
|
||||
|
||||
sensiPlot <- ggplot(data.frame(name=varNames,y=S/100),aes(x=name,y=y))+
|
||||
geom_bar(stat = 'identity')+
|
||||
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
|
||||
xlab(NULL)+
|
||||
ylab(NULL)+
|
||||
ggtitle("Sensitivity")+
|
||||
scale_y_continuous(labels = scales::percent,limits=c(0,1))
|
||||
print(sensiPlot)
|
||||
ggsave(plotName,dpi=dpi)
|
||||
return(S)
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(is.null(monteCarloFile)){
|
||||
M <- musoMont(parameters = parameters,
|
||||
settings = settings,
|
||||
inputDir = inputDir,
|
||||
outLoc = outLoc,
|
||||
iterations = iterations,
|
||||
preTag = preTag,
|
||||
outputType = outputType,
|
||||
fun = fun,
|
||||
varIndex = varIndex,
|
||||
skipSpinup = skipSpinup
|
||||
)
|
||||
yInd <- grep("mod.", colnames(M))[varIndex]
|
||||
parNames <- grep("mod.",colnames(M), invert=TRUE, value = TRUE)
|
||||
M <- M[,c(grep("mod.", colnames(M),invert=TRUE),yInd)]
|
||||
|
||||
return(doSensi(M))
|
||||
|
||||
} else {
|
||||
M <- read.csv(monteCarloFile)
|
||||
return(doSensi(M))
|
||||
}
|
||||
}
|
||||
@ -1,284 +0,0 @@
|
||||
#'plot the BBGCMuso output
|
||||
#'
|
||||
#' This function runs the BBGC-MuSo model and reads in its outputfile in a very structured way, and after that plot the results automaticly
|
||||
#'
|
||||
#' @author Roland Hollos, Dora Hidy
|
||||
#' @param settings You have to run the setupMuso function before rungetMuso. It is its output which contains all of the necessary system variables. It sets the whole environment
|
||||
#' @param timee The required timesteps in the modell output. It can be "d", if it is daily, "m", if it's monthly, "y", it it is yearly
|
||||
#' @param debugging If it is TRUE, it copies the log file to a Log directory to store it, if it is stamplog it contatenate a number before the logfile, which is one more than the maximum of the represented ones in the LOG directory. If it is true or stamplog it collects the "wrong" logfiles
|
||||
#' @param keepEpc If TRUE, it keeps the epc file and stamp it, after these copies it to the EPCS directory. If debugging True or false, it copies the wrong epc files to the wrong epc directory.
|
||||
#' @param export if it is yes or you give a filename here, it converts the output to the specific extension. For example, if you set export to "example.csv", it converts the output to "csv", if you set it to "example.xls" it converts to example.xls with the xlsx package. If it is not installed it gives back a warning message and converts it to csv.
|
||||
#' @param silent If you set it TRUE all off the modells output to the screen will be suppressed. It can be usefull, because it increases the model-speed.
|
||||
#' @param aggressive It deletes every possible modell-outputs from the previous modell runs.
|
||||
#' @param variable column number of the variable which should be plottedor "all" if you have less than 10 variables. In this case it will plot everything in a matrix layout
|
||||
#' @param leapYear Should the function do a leapyear correction on the outputdata? If TRUE, then the 31.12 day will be doubled.
|
||||
#' @param logfilename If you want to set a specific name for your logfiles you can set this via logfile parameter
|
||||
#' @param plotType There are two options now: continious time series("cts") or disctrete time series("dts")
|
||||
#' @param skipSpinup If TRUE, calibMuso wont do spinup simulation
|
||||
#' @return It depends on the export parameter. The function returns with a matrix with the modell output, or writes this in a file, which is given previously
|
||||
#' @usage plotMuso(settings, variable,
|
||||
#' timee="d", silent=TRUE,
|
||||
#' debugging=FALSE, keepEpc=FALSE,
|
||||
#' logfilename=NULL, aggressive=FALSE,
|
||||
#' leapYear=FALSE, export=FALSE)
|
||||
#' @import ggplot2
|
||||
#' @import dplyr
|
||||
#' @import tidyr
|
||||
#' @export
|
||||
|
||||
plotMuso <- function(settings=NULL,
|
||||
variable=1,
|
||||
##compare,
|
||||
##plotname,
|
||||
timee="d",
|
||||
silent=TRUE,
|
||||
debugging=FALSE,
|
||||
keepEpc=FALSE,
|
||||
logfilename=NULL,
|
||||
aggressive=FALSE,
|
||||
leapYear=FALSE,
|
||||
plotName=NULL,
|
||||
plotType="cts",
|
||||
layerPlot=FALSE,
|
||||
colour="blue",
|
||||
skipSpinup=TRUE,
|
||||
fromData=FALSE,
|
||||
timeFrame="day",
|
||||
groupFun=mean,
|
||||
dpi=300){
|
||||
|
||||
if( plotType!="cts" && plotType != "dts"){
|
||||
warning(paste0("The plotType ", plotType," is not implemented, plotType is set to cts"))
|
||||
plotType <- "cts"
|
||||
}
|
||||
|
||||
if(is.null(settings)){
|
||||
settings <- setupMuso()
|
||||
}
|
||||
|
||||
numberOfYears <- settings$numYears
|
||||
startYear <- settings$startYear
|
||||
## musoData <- rungetMuso(settings=settings,
|
||||
## silent=silent,
|
||||
## timee=timee,
|
||||
## debugging=debugging,
|
||||
## keepEpc=keepEpc,
|
||||
## logfilename=logfilename,
|
||||
## export=export)
|
||||
|
||||
groupByTimeFrame <- function(data, timeFrame, groupFun){
|
||||
datas <- data %>%
|
||||
group_by(year) %>%
|
||||
summarize(variable=groupFun(eval(parse(text=variable))))
|
||||
datas[,1]<-as.numeric(unlist(datas[,1]))
|
||||
colnames(datas) <- c("date",variable)
|
||||
datas
|
||||
}
|
||||
|
||||
if(fromData){
|
||||
Reva <- tryCatch(getdailyout(settings), #(:INSIDE: getOutput.R )
|
||||
error = function (e){
|
||||
setwd((whereAmI))
|
||||
stop("Cannot read binary output, please check if the output type is set 2 in the ini files!")})
|
||||
colnames(Reva) <- unlist(settings$outputVars[[1]])
|
||||
rownames(Reva) <- NULL
|
||||
musoData <- cbind(musoDate(startYear = startYear,numYears = numberOfYears,combined = TRUE,corrigated=FALSE),
|
||||
rep(1:365,numberOfYears),
|
||||
musoDate(startYear = startYear,numYears = numberOfYears,combined = FALSE,corrigated=FALSE),as.data.frame(Reva))
|
||||
colnames(musoData)[1:5]<-c("date","yearDay","year","day","month")
|
||||
musoData <-musoData %>%
|
||||
mutate(date=as.Date(as.character(date),"%d.%m.%Y"))
|
||||
} else {
|
||||
if(!is.element("cum_yieldC_HRV",unlist(settings$outputVars[[1]]))){
|
||||
musoData <- calibMuso(settings,silent = TRUE,skipSpinup=skipSpinup) %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("date") %>%
|
||||
mutate(date2=date,date=as.Date(date,"%d.%m.%Y")) %>%
|
||||
tidyr::separate(date2,c("day","month","year"),sep="\\.")
|
||||
if(timeFrame!="day"){
|
||||
musoData <- tryCatch(groupByTimeFrame(data=musoData, timeFrame = timeFrame, groupFun = groupFun),
|
||||
error=function(e){stop("The timeFrame or the gropFun is not found")})
|
||||
}} else {
|
||||
musoData <- calibMuso(settings,silent = TRUE,skipSpinup=skipSpinup) %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("date") %>%
|
||||
mutate(date2=date,date=as.Date(date,"%d.%m.%Y"),
|
||||
yearDay=rep(1:365,numberOfYears), cum_yieldC_HRV=cum_yieldC_HRV*22.22) %>%
|
||||
tidyr::separate(date2,c("day","month","year"),sep="\\.")
|
||||
if(timeFrame!="day"){
|
||||
musoData <- tryCatch(groupByTimeFrame(data=musoData, timeFrame = timeFrame, groupFun = groupFun),
|
||||
error=function(e){stop("The timeframe or the gropFun is not found")})
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
## numVari <- ncol(musoData)
|
||||
numVari <- ncol(musoData)-5
|
||||
|
||||
pointOrLineOrPlot <- function(musoData, variableName, plotType="cts", expandPlot=FALSE, plotName=NULL){
|
||||
if(!expandPlot){
|
||||
if(plotType=="cts"){
|
||||
if(length(variableName)==1){
|
||||
p <- ggplot(musoData,aes_string("date",variableName))+geom_line(colour=colour)+theme(axis.title.x=element_blank())
|
||||
if(!is.null(plotName)){
|
||||
ggsave(as.character(plotName), plot = p)
|
||||
p
|
||||
}
|
||||
p
|
||||
} else {
|
||||
p <- musoData %>%
|
||||
select(c("date", variableName))%>%
|
||||
gather(., key= outputs, value = bla, variableName) %>%
|
||||
# head %>%
|
||||
ggplot(aes(x=date,y=bla))+
|
||||
facet_wrap(~ outputs, scales = "free_y",ncol=1) +
|
||||
geom_line(colour=colour)+
|
||||
theme(
|
||||
axis.title.y = element_blank()
|
||||
)
|
||||
if(!is.null(plotName)){
|
||||
ggsave(as.character(plotName), plot = p)
|
||||
}
|
||||
p
|
||||
}
|
||||
} else {
|
||||
if(length(variableName)==1){
|
||||
p <- ggplot(musoData,aes_string("date",variableName))+geom_point(colour=colour)+theme(axis.title.x=element_blank())
|
||||
if(!is.null(plotName)){
|
||||
ggsave(as.character(plotName),p)
|
||||
}
|
||||
p
|
||||
} else{
|
||||
p <- musoData %>%
|
||||
select(c("date",variableName))%>%
|
||||
gather(., key= outputs, value = bla,variableName) %>%
|
||||
# head %>%
|
||||
ggplot(aes(x=date,y=bla))+
|
||||
facet_wrap(~ outputs, scales = "free_y",ncol=1) +
|
||||
geom_line(colour=colour)+
|
||||
theme(
|
||||
axis.title.y = element_blank()
|
||||
)
|
||||
if(!is.null(plotName)){
|
||||
ggsave(as.character(plotName),p)
|
||||
}
|
||||
p
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if(!is.null(plotName)){
|
||||
stop("Cannot save a single plot layer to a graphics device")
|
||||
}
|
||||
|
||||
if(plotType=="cts"){
|
||||
if(length(variableName)==1){
|
||||
geom_line(data=musoData, colour=colour, aes_string("date",variableName))
|
||||
|
||||
} else {
|
||||
stop("you cannot add layers for multiple plots")
|
||||
}
|
||||
} else {
|
||||
if(length(variableName)==1){
|
||||
geom_point(data=musoData, colour=colour, aes_string("date",variableName))
|
||||
} else{
|
||||
stop("you cannot add layers for multiple plots")
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
variableName <- as.character(settings$outputVars[[1]])[variable]
|
||||
if(is.character(variable)){
|
||||
if(identical(variable,"all")){
|
||||
variableName <- as.character(settings$outputVars[[1]])
|
||||
} else {
|
||||
if(identical(character(0),setdiff(variable,as.character(settings$outputVars[[1]])))){
|
||||
variableName <- variable
|
||||
} else {
|
||||
stop("The symmetric difference of the set of the output variables specified in the ini files and the set specified with your variable parameter is not the empty set.")
|
||||
}
|
||||
}
|
||||
|
||||
if(length(variableName)>8){
|
||||
warning("Too many variables to plot, the output quality can be poor")
|
||||
}
|
||||
|
||||
} else {
|
||||
|
||||
if(prod(sapply(variable,function(x){
|
||||
return(x >= 0 && x <= numVari)
|
||||
}))){
|
||||
variableName <- as.character(settings$outputVars[[1]])[variable]
|
||||
} else {
|
||||
stop("Not all members of the variable parameter are among the output variables")
|
||||
}}
|
||||
|
||||
pointOrLineOrPlot(musoData = musoData,
|
||||
variableName = variableName,
|
||||
plotType = plotType,
|
||||
expandPlot = layerPlot,
|
||||
plotName = plotName)
|
||||
}
|
||||
|
||||
#'plot the BBGCMuso output with data
|
||||
#'
|
||||
#' This function runs the BBGC-MuSo model and reads in its outputfile in a very structured way, and after that plot the results automaticly along with a given measurement
|
||||
#'
|
||||
#' @author Roland Hollos, Dora Hidy
|
||||
#' @param settings You have to run the setupMuso function before rungetMuso. It is its output which contains all of the necessary system variables. It sets the whole environment
|
||||
#' @param sep This is the separator used in the measurement file
|
||||
#' @param savePlot It it is specified, the plot will be saved in a format specified with the immanent extension
|
||||
#' @param variable The name of the output variable to plot
|
||||
#' @param NACHAR This is not implemented yet
|
||||
#' @param csvFile The file of the measurement. It must contain a header.
|
||||
#' @usage plotMuso(settings, variable,
|
||||
#' timee="d", silent=TRUE,
|
||||
#' debugging=FALSE, keepEpc=FALSE,
|
||||
#' logfilename=NULL, aggressive=FALSE,
|
||||
#' leapYear=FALSE, export=FALSE)
|
||||
#' @import ggplot2
|
||||
#' @export
|
||||
plotMusoWithData <- function(csvFile, variable, NACHAR=NA, settings=NULL, sep=",", savePlot=NULL,colour=c("black","blue")){
|
||||
if(!is.na(NACHAR)){
|
||||
warning("NACHAR is not implemented yet")
|
||||
}
|
||||
if(is.null(settings)){
|
||||
settings <- setupMuso()
|
||||
}
|
||||
|
||||
numberOfYears <- settings$numYears
|
||||
startYear <- settings$startYear
|
||||
yearVec <- seq(from = startYear, length=numberOfYears,by=1)
|
||||
|
||||
|
||||
data <- read.table(csvFile,header = TRUE, sep = ",") %>%
|
||||
select(variable)
|
||||
|
||||
baseData <- calibMuso(settings,silent=TRUE) %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("date") %>%
|
||||
dplyr::mutate(date2=date,date=as.Date(date,"%d.%m.%Y"),yearDay=rep(1:365,numberOfYears)) %>%
|
||||
tidyr::separate(date2,c("day","month","year"),sep="\\.")
|
||||
baseData <- cbind(baseData,data)
|
||||
colnames(baseData)[ncol(baseData)] <- "measuredData"
|
||||
|
||||
p <- baseData %>%
|
||||
ggplot(aes_string("date",variable)) +
|
||||
geom_line(colour=colour[1]) +
|
||||
geom_point(colour=colour[2], aes(date,measuredData)) +
|
||||
labs(y = paste0(variable,"_measured"))+
|
||||
theme(axis.title.x = element_blank())
|
||||
if(!is.null(savePlot)){
|
||||
ggsave(savePlot,p)
|
||||
return(p)
|
||||
} else {
|
||||
return(p)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -1,13 +0,0 @@
|
||||
#' putOutVars
|
||||
#'
|
||||
#'This function is for adding variables in the inifiles.
|
||||
#'
|
||||
#' @author Roland Hollos
|
||||
#' @param IniFile The name of the normal ini file.
|
||||
#' @param outputVars List of the output codes
|
||||
#' @keywords internal
|
||||
putOutVars <- function(IniFile,outputVars){
|
||||
ini <- readLines(IniFile)
|
||||
sapply(outputVars,musoMapping)
|
||||
}
|
||||
|
||||
@ -1,274 +0,0 @@
|
||||
#' setupMuso
|
||||
#'
|
||||
#' This funcion is fundamental for the BiomBGC-MuSo modell related functions like spinupMuso, normalMuso, rungetMuso, because it sets the modells environment.
|
||||
#'
|
||||
#' @author Roland Hollos
|
||||
#' @param parallel Do you want to run multiple modell paralelly, if yes, set this variable to TRUE
|
||||
#' @param executable This parameter 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 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 outputLoc 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 inputLoc Usually it is the root directory, where you put the iniFiles for the modell
|
||||
#' @param metInput Via metInput parameter, you can tell the modell where are the meteorological files. As default it reads this from the iniFiles.
|
||||
#' @param CO2Input Via CO2 parameter, you can tell the modell where are the CO2 data files. As default it reads this from the iniFiles.
|
||||
#' @param plantInput 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 thinInput 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 mowInput 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 grazInput 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 harvInput 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.
|
||||
#' @param plougInput Via ploughing parameter, you can tell the modell where are the data files, which contains the ploughing informations. As default it reads this from the iniFiles.
|
||||
#' @param fertInput Via fertilizing parameter, you can tell the modell where are the fertilizing data files, which contains the fertilizing informations. As default it reads this from the iniFiles.
|
||||
#' @param irrInput Via irrigation parameter, you can tell the modell where are the data files, which contains the irrigation informations. As default it reads this from the iniFiles.
|
||||
#' @param nitInput Via this parameter, you can tell the modell where are the NO2 data files. As default it reads this from the iniFiles.
|
||||
#' @param iniInput Via this parameter, you can tell the modell where are the ini files. As default it reads this from the iniFiles.
|
||||
#' @param epcInput Via this parameter, you can tell the modell where are the epc data files. As default it reads this from the iniFiles.
|
||||
#' @param modelOutputs This parameter is the list of the codes of outpu
|
||||
#' @usage setupMuso(executable=NULL, parallel = F, calibrationPar =c(1),
|
||||
#' outputLoc=NULL, inputLoc=NULL,
|
||||
#' metInput=NULL, CO2Input=NULL,
|
||||
#' plantInput=NULL, thinInput=NULL,
|
||||
#' mowInput=NULL, grazInput=NULL,
|
||||
#' harvInput=NULL, plougInput=NULL,
|
||||
#' fertInput=NULL, irrInput=NULL,
|
||||
#' nitInput=NULL, iniInput=NULL, epcInput=NULL)
|
||||
#' @return The output is a the modell setting list wich contains the following elements:
|
||||
#' executable, calibrationPar, outputLoc, outputName, inputLoc, iniInput, metInput, epcInput,thinInput,CO2Input, mowInput, grazInput, harvInput, plougInput, fertInput,rrInput, nitInput, inputFiles, numData, startyear, numYears, outputVars
|
||||
#' @export
|
||||
setupMuso <- function(executable=NULL,
|
||||
parallel = F,
|
||||
calibrationPar =c(1),
|
||||
outputLoc=NULL,
|
||||
modelOutputs=NULL,
|
||||
inputLoc=NULL,
|
||||
metInput=NULL,
|
||||
CO2Input=NULL,
|
||||
plantInput=NULL,
|
||||
thinInput=NULL,
|
||||
mowInput=NULL,
|
||||
grazInput=NULL,
|
||||
harvInput=NULL,
|
||||
plougInput=NULL,
|
||||
fertInput=NULL,
|
||||
irrInput=NULL,
|
||||
nitInput=NULL,
|
||||
iniInput=NULL,
|
||||
epcInput=NULL,
|
||||
mapData=NULL,
|
||||
leapYear=FALSE,
|
||||
version=5
|
||||
){
|
||||
|
||||
Linuxp <-(Sys.info()[1]=="Linux")
|
||||
writep <- 0
|
||||
|
||||
if(is.null(mapData)&version==4){
|
||||
mData <- mMapping4
|
||||
}
|
||||
|
||||
inputParser <- function(string,fileName,counter,value=TRUE){
|
||||
unlist(strsplit(grep(string,fileName,value=TRUE),"[\ \t]"))[counter]
|
||||
}
|
||||
|
||||
outMaker <- function(inputVar,grepString,filep){
|
||||
tempVar <- eval(parse(text=inputVar))
|
||||
if(is.null(tempVar)){
|
||||
writep <<- writep+1
|
||||
if(filep)
|
||||
{
|
||||
tempVar["spinup"] <- file.path(inputLoc,inputParser(string=grepString,fileName=iniFiles$spinup,counter=1,value=TRUE))
|
||||
tempVar["normal"] <- file.path(inputLoc,inputParser(string=grepString,fileName=iniFiles$normal,counter=1,value=TRUE))
|
||||
} else {
|
||||
tempVar["spinup"] <- inputParser(string=grepString,fileName=iniFiles$spinup,counter=1,value=TRUE)
|
||||
tempVar["normal"] <- inputParser(string=grepString,fileName=iniFiles$normal,counter=1,value=TRUE)
|
||||
|
||||
}
|
||||
|
||||
} else {
|
||||
iniFiles$spinup[grep(grepString,iniFiles$spinup)] <<- paste0(tempVar[1],"\t ",grepString)
|
||||
|
||||
if(length(tempVar)==2){
|
||||
iniFiles$normal[grep(" grepString",iniFiles$normal)] <<- paste0(tempVar[2],"\t ",grepString)
|
||||
}
|
||||
}
|
||||
return(tempVar)
|
||||
}
|
||||
|
||||
if(is.null(inputLoc)){
|
||||
inputLoc<- normalizePath("./")
|
||||
} else{
|
||||
inputLoc <- normalizePath(inputLoc)
|
||||
}
|
||||
|
||||
#iniChangedp <- FALSE
|
||||
|
||||
if(is.null(iniInput)){
|
||||
spinups<-grep("s.ini$",list.files(inputLoc),value=TRUE)
|
||||
normals<-grep("n.ini$",list.files(inputLoc),value=TRUE)
|
||||
|
||||
if(length(spinups)==1){
|
||||
iniInput[1]<-file.path(inputLoc,spinups)
|
||||
} else {stop("There are multiple or no spinup ini files, please choose")}
|
||||
|
||||
|
||||
if(length(normals)==1){
|
||||
iniInput[2]<-file.path(inputLoc,normals)
|
||||
} else {stop("There are multiple or no normal ini files, please choose")}
|
||||
|
||||
}
|
||||
|
||||
##read the ini files for the further changes
|
||||
|
||||
iniFiles<-lapply(iniInput, function (x) readLines(x,-1))
|
||||
iniFiles[[1]] <- gsub("\\","/", iniFiles[[1]],fixed=TRUE) #replacing \ to /
|
||||
iniFiles[[2]] <- gsub("\\","/", iniFiles[[2]],fixed=TRUE) #replacing \ to /
|
||||
names(iniFiles) <- c("spinup","normal")
|
||||
|
||||
|
||||
|
||||
inputs <- lapply(1:nrow(grepHelper), function (x) {
|
||||
|
||||
outMaker(grepHelper[x,1],grepHelper[x,2],grepHelper[x,3])
|
||||
|
||||
})
|
||||
names(inputs) <- grepHelper$inputVar
|
||||
## grepHelper is in sysdata.rda it is a table like this:
|
||||
##
|
||||
## inputVar string isFile
|
||||
## 1 epcInput EPC file name TRUE
|
||||
## 2 metInput met file name TRUE
|
||||
## 3 CO2Input CO2 file TRUE
|
||||
## 4 nitInput N-dep file TRUE
|
||||
## 5 thinInput do THINNING FALSE
|
||||
## 6 plantInput do PLANTING FALSE
|
||||
## 7 mowInput do MOWING FALSE
|
||||
## 8 grazInput do GRAZING FALSE
|
||||
## 9 harvInput do HARVESTING FALSE
|
||||
## 10 plougInput do PLOUGHING FALSE
|
||||
## 11 fertInput do FERTILIZING FALSE
|
||||
## 12 irrInput do IRRIGATION FALSE
|
||||
# return(inputs) debug element
|
||||
|
||||
|
||||
if(is.null(mapData)){
|
||||
|
||||
outIndex<-grep("DAILY_OUTPUT",iniFiles[[2]])+1
|
||||
numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]"))[1])
|
||||
dailyVarCodes<-tryCatch(iniFiles[[2]][(outIndex+1):(outIndex+numVar)],
|
||||
error = function(e){
|
||||
stop("Cannot read indexes of output variables from the normal ini file, please make sure you have not skiped a line after the flag: \"DAILY_OUTPUT\"")
|
||||
})
|
||||
dailyVarnames<-lapply(dailyVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1]))
|
||||
|
||||
outIndex<-grep("ANNUAL_OUTPUT",iniFiles[[2]])+1
|
||||
numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][outIndex],"[\ \t]"))[1])
|
||||
annualVarCodes<-iniFiles[[2]][(outIndex+1):(outIndex+numVar)]
|
||||
annualVarnames<-lapply(annualVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1]))
|
||||
outputVars<-list(dailyVarnames,annualVarnames)} else {
|
||||
|
||||
c<-grep("DAILY_OUTPUT",iniFiles[[2]])+1
|
||||
numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][c],"[\ \t]"))[1])
|
||||
dailyVarCodes<-iniFiles[[2]][(c+1):(c+numVar)]
|
||||
dailyVarnames<-lapply(dailyVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1],mapData))
|
||||
|
||||
c<-grep("ANNUAL_OUTPUT",iniFiles[[2]])+1
|
||||
numVar<-as.numeric(unlist(strsplit(iniFiles[[2]][c],"[\ \t]"))[1])
|
||||
annualVarCodes<-iniFiles[[2]][(c+1):(c+numVar)]
|
||||
annualVarnames<-lapply(annualVarCodes, function(x) musoMapping(unlist(strsplit(x,"[\ \t]"))[1],mapData))
|
||||
outputVars<-list(dailyVarnames,annualVarnames)
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
if(is.null(executable)){
|
||||
if(Linuxp){
|
||||
executable<-file.path(inputLoc,"muso")
|
||||
} else {
|
||||
executable<-file.path(inputLoc,"muso.exe")
|
||||
}
|
||||
} else {
|
||||
file.copy(executable,inputLoc)
|
||||
}
|
||||
outputName <- character(2)
|
||||
outputName[1] <- basename(unlist(strsplit(iniFiles[[1]][grep("OUTPUT_CONTROL",iniFiles[[1]])+1],"[\ \t]"))[1])
|
||||
outputName[2] <- basename(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]])+1],"[\ \t]"))[1])
|
||||
## outputName <- unlist(strsplit(grep("output",grep("prefix",iniFiles[[2]],value=TRUE),value=TRUE),"[\ \t]"))[1]
|
||||
##THIS IS AN UGLY SOLUTION, WHICH NEEDS AN UPGRADE!!! FiXED (2017.09.11)
|
||||
## outputName <- unlist(strsplit(grep("prefix for output files",iniFiles[[2]],value=TRUE),"[\ \t]"))[1]
|
||||
if(is.null(outputName)){
|
||||
stop("I cannot find outputName in your default ini file \n Please make sure that the line wich contains the name also contains the prefix and the output keywords!")
|
||||
|
||||
}
|
||||
## outputName<-unlist(read.table(iniInput[2],skip=93,nrows = 1))[1]
|
||||
|
||||
|
||||
if(is.null(outputLoc)){
|
||||
## outputLoc<-paste((rev(rev(unlist(strsplit(outputName,"/")))[-1])),collapse="/")
|
||||
outputLoc <- dirname(unlist(strsplit(iniFiles[[2]][grep("OUTPUT_CONTROL",iniFiles[[2]])+1],"[\ \t]"))[1])
|
||||
if(substr(outputLoc,start = 1,stop = 1)!="/"){
|
||||
##if the outputName is not absolute path make it absolute
|
||||
outputLoc <- file.path(inputLoc,outputLoc)
|
||||
}
|
||||
} else {
|
||||
outputLoc <- normalizePath(outputLoc)
|
||||
}
|
||||
|
||||
|
||||
|
||||
inputFiles<-c(iniInput,epcInput,metInput)
|
||||
numData<-rep(NA,3)
|
||||
numYears <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]])+2],"[\ \t]"))[1])
|
||||
## numYears<-unlist(read.table(iniInput[2],skip = 14,nrows = 1)[1])
|
||||
numValues <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("DAILY_OUTPUT",iniFiles[[2]])+1],"[\ \t]"))[1])
|
||||
## numValues will be replaced to numVar
|
||||
## numValues<-unlist(read.table(iniInput[2],skip=102,nrows = 1)[1])
|
||||
startYear <- as.numeric(unlist(strsplit(iniFiles[[2]][grep("TIME_DEFINE",iniFiles[[2]])+3],"[\ \t]"))[1])
|
||||
numData[1] <- numValues * sumDaysOfPeriod(startYear,numYears,corrigated=leapYear)
|
||||
|
||||
numData[2] <- numYears * numValues*12
|
||||
numData[3] <- numYears * numValues
|
||||
|
||||
##Writing out changed ini-file
|
||||
|
||||
writeLines(iniFiles[[1]],iniInput[1])
|
||||
writeLines(iniFiles[[2]],iniInput[2])
|
||||
|
||||
suppressWarnings(file.remove(paste0(file.path(outputLoc,outputName[1]),".log")))
|
||||
## I use file.path additionally because We do not know if outputLoc ends or not to "/"
|
||||
suppressWarnings(file.remove(paste0(file.path(outputLoc,outputName[2]),".log")))
|
||||
|
||||
settings = list(executable = executable,
|
||||
calibrationPar = calibrationPar,
|
||||
outputLoc=outputLoc,
|
||||
outputNames=outputName,
|
||||
inputLoc=inputLoc,
|
||||
iniInput=iniInput,
|
||||
metInput=inputs$metInput,
|
||||
epcInput=inputs$epcInput,
|
||||
thinInput=inputs$thinInput,
|
||||
CO2Input=inputs$CO2Input,
|
||||
mowInput=inputs$mowInput,
|
||||
grazInput=inputs$grazInput,
|
||||
harvInput=inputs$harvInput,
|
||||
plougInput=inputs$plougInput,
|
||||
fertInput=inputs$fertInput,
|
||||
irrInput=inputs$irrInput,
|
||||
nitInput=inputs$nitInput,
|
||||
inputFiles=inputFiles,
|
||||
numData=numData,
|
||||
startYear=startYear,
|
||||
numYears=numYears,
|
||||
outputVars=outputVars
|
||||
)
|
||||
|
||||
if(writep!=nrow(grepHelper)){
|
||||
writeLines(iniFiles[[1]],iniInput[[1]])
|
||||
if(inputs$epcInput[1]!=inputs$epc$Input[2]){ #Change need here
|
||||
writeLines(iniFiles[[2]],iniInput[[2]])
|
||||
}
|
||||
}
|
||||
return(settings)
|
||||
|
||||
}
|
||||
|
||||
Binary file not shown.
@ -3,3 +3,4 @@ basePackages <- c("dplyr","digest","ggplot2","shiny","latex2expr",
|
||||
install.packages(basePackages)
|
||||
install.packages("https://github.com/hollorol/RBBGCMuso/raw/master/RBBGCMuso_0.5.0.0-0.zip", repos = NULL, type = "win.binary")
|
||||
|
||||
copyMu
|
||||
Loading…
Reference in New Issue
Block a user