diff --git a/RBBGCMuso/R/calibMuso.R b/RBBGCMuso/R/calibMuso.R index 400c5e8..6d1fc77 100644 --- a/RBBGCMuso/R/calibMuso.R +++ b/RBBGCMuso/R/calibMuso.R @@ -99,25 +99,31 @@ calibMuso <- function(settings=NULL, calibrationPar=NULL, - if(aggressive==TRUE){ - cleanupMuso(location=outputLoc,deep = TRUE) + if(aggressive == TRUE){ + cleanupMuso(location = outputLoc,deep = TRUE) } - toModif<-c(epc[2],iniInput[2]) + toModif <- c(epc[2],iniInput[2]) # if(!modifyOriginal & (!is.null(parameters) | !is.null(outVars))) # { toModif <- sapply(toModif, function (x){ - paste0(tools::file_path_sans_ext(basename(x)),"-tmp.",tools::file_ext(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) + changemulline(filePaths = c(epc[2], iniInput[2]), + calibrationPar = calibrationPar, + contents = parameters, + fileOut = toModif, + fileToChange = fileToChange, + modifyOriginal = modifyOriginal) } @@ -150,7 +156,9 @@ calibMuso <- function(settings=NULL, calibrationPar=NULL, # sapply(c(iniInput,epc),) # # } - + if(modifyOriginal){ + iniInput[2] <- toModif[2] + } if(!skipSpinup) { @@ -250,7 +258,6 @@ calibMuso <- function(settings=NULL, calibrationPar=NULL, 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){ diff --git a/RBBGCMuso/R/calibration.R b/RBBGCMuso/R/calibration.R index 460ae3e..62829af 100644 --- a/RBBGCMuso/R/calibration.R +++ b/RBBGCMuso/R/calibration.R @@ -173,13 +173,13 @@ optiMuso <- function(measuredData, parameters = NULL, startDate = NULL, preservedCalib<- read.csv("preservedCalib.csv") p<-list() preservedCalib <- preservedCalib[-1,] + preservedCalib <- preservedCalib[!is.na(preservedCalib$likelihood),] dontInclude <-c((ncol(preservedCalib)-1),ncol(preservedCalib)) # for(i in seq_along(colnames(preservedCalib)[-dontInclude])){ # p[[i]] <- ggplot(as.data.frame(preservedCalib),aes_string(colnames(preservedCalib)[i],"likelihood")) + # geom_point(shape='.',size=1,alpha=0.8) # } - unfilteredLikelihood <- preservedCalib$likelihood preservedCalib <- preservedCalib[preservedCalib$likelihood>quantile(preservedCalib$likelihood,0.95),] optRanges <-t(apply(preservedCalib,2,function(x) quantile(x,c(0.05,0.5,0.95)))) @@ -208,6 +208,8 @@ dev.off() # file.rename(tempEpc, "optimizedEpc.epc") # return(preservedCalib[maxLikelihoodPlace,]) write.csv(optRanges,"optRanges.csv") + # is.num <- sapply(head(optRanges,-2), is.numeric) + # optRanges[is.num] <- lapply(optRanges[is.num], round, 4) return(head(optRanges,n=-2)) } diff --git a/RBBGCMuso/R/quickeffect.R b/RBBGCMuso/R/quickeffect.R index 1b8b150..549a978 100644 --- a/RBBGCMuso/R/quickeffect.R +++ b/RBBGCMuso/R/quickeffect.R @@ -15,7 +15,7 @@ #' @importFrom tidyr separate #' @export -musoQuickEffect <- function(settings = NULL, calibrationPar = NULL, startVal, endVal, nSteps = 1, fileToChange="epc",modifyOriginal=TRUE, outVar, parName = "parVal"){ +musoQuickEffect <- function(settings = setupMuso(), calibrationPar = NULL, startVal, endVal, nSteps = 1, fileToChange="epc",modifyOriginal=TRUE, outVar, parName = "parVal", yearNum=1, year=(settings$startYear + yearNum -1)){ if(is.character(outVar)){ varNames <- as.data.frame(musoMappingFind(outVar)) @@ -33,9 +33,6 @@ musoQuickEffect <- function(settings = NULL, calibrationPar = NULL, startVal, e outVarIndex<-outVar } - if(is.null(settings)){ - settings <- setupMuso() - } if(is.null(calibrationPar)){ calibrationPar <- settings$calibrationPar } @@ -43,21 +40,43 @@ musoQuickEffect <- function(settings = NULL, calibrationPar = NULL, startVal, e parVals <- seq(startVal, endVal, length = (nSteps + 1)) parVals <- dynRound(startVal, endVal, seqLen = (nSteps + 1)) a <- do.call(rbind,lapply(parVals, function(parVal){ - calResult <- tryCatch(calibMuso(settings = settings,calibrationPar = calibrationPar,modifyOriginal=modifyOriginal, parameters = parVal, outVars = outVarIndex, silent = TRUE,fileToChange = fileToChange), error = function(e){NA}) - if(all(is.na(calResult))){ + calResult <- tryCatch(calibMuso(settings = settings,calibrationPar = calibrationPar, + modifyOriginal = modifyOriginal, + parameters = parVal, + outVars = outVarIndex, + silent = TRUE, + fileToChange = fileToChange), error = function(e){NULL}) + if(is.null(calResult)){ b <- cbind(rep(NA,365),parVal) - rownames(b) <- tail(musoDate(startYear = settings$startYear, numYears = settings$numYears),365) + rownames(b) <- musoDate(startYear = year, numYears = 1) colnames(b)[1] <- varNames return(b) } else { - return(cbind(tail(calResult,365), parVal)) + m <- as.data.frame(calResult[musoDate(startYear = year, numYears = 1),]) + colnames(m) <- colnames(calResult) + return(cbind(m, parVal)) } })) - a %<>% tbl_df %>% mutate(date=as.Date(rownames(a),"%d.%m.%Y")) %>% select(date,as.character(varNames),parVal) print(suppressWarnings(ggplot(data = a, aes_string(x= "date", y= varNames))+geom_line(aes(alpha = factor(parVal))) + labs(y=varNames, alpha = parName) + scale_alpha_discrete(range=c(0.25,1)))) } +# calma <- calibMuso(settings = settings,calibrationPar = calibrationPar, +# modifyOriginal = modifyOriginal, +# parameters = parVal, +# outVars = outVarIndex, +# silent = TRUE, +# fileToChange = fileToChange) +# plot(calma[,1]) +# calma <- calibMuso(settings = settings,calibrationPar = calibrationPar, +# modifyOriginal = modifyOriginal, +# parameters = parVal, +# silent = TRUE, +# fileToChange = fileToChange) +# calm <- calibMuso(calibrationPar=calibrationPar,parameters=parVal,modifyOriginal=TRUE) +# plot(x=as.Date(musoDate(2015,numYears=1),"%d.%m.%Y"),y=calm[musoDate(2015,numYears=1),"daily_gpp"],type="l") +# calibrationPar +# parVal