diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 6b6ef20..dd56fc7 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -55,6 +55,7 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s naVal, postProcString, i) , error = function(e){ # browser() + writeLines(as.character(e),"error.txt") writeLines(as.character(iterations),"progress.txt") }) @@ -132,7 +133,7 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s musoGlue(results, parameters=parameters,settings=settings, w=w, lg=lg) liks <- results[,sprintf("%s_likelihood",names(likelihood))] epcIndexes <- future::value(fut[[1]], stdout = FALSE, signal=FALSE) - if(ncol(liks) == 1){ + if(is.null(dim(liks)) || dim(liks)[2] == 1 ){ ml_place <- which.max(liks) } else { ml_place <- which.max(as.matrix(liks) %*% as.matrix(w)) @@ -249,6 +250,7 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, pretag <- file.path(outLoc,preTag) ##reading the original epc file at the specified ## row numbers + # browser() print("optiMuso is randomizing the epc parameters now...",quote = FALSE) if(iterations < 3000){ randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000,sourceFile=sourceFile) @@ -314,11 +316,21 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, # colNumb <- length(settings$dailyVarCodes) + 1 # } + singleVarCalib <- FALSE + if(is.null(dim(randValues))){ + singleVarCalib <- TRUE + } for(i in 2:(iterations+1)){ - + + if(!singleVarCalib){ + parameters <- randValues[(i-1),] + } else { + parameters <- randValues[i] + } + tmp <- tryCatch(calibMuso(settings = settings, - parameters = randValues[(i-1),], + parameters = parameters, silent= TRUE, skipSpinup = skipSpinup, modifyOriginal=modifyOriginal, postProcString = postProcString), error = function (e) NULL) if(is.null(tmp)){ @@ -332,7 +344,7 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, musoCodeToIndex = musoCodeToIndex, uncert = uncert) } - partialResult[1:numParameters] <- randValues[(i-1),] + partialResult[1:numParameters] <- parameters 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")) diff --git a/RBBGCMuso/R/calibration.R b/RBBGCMuso/R/calibration.R index db6e4af..0ad984b 100644 --- a/RBBGCMuso/R/calibration.R +++ b/RBBGCMuso/R/calibration.R @@ -248,7 +248,10 @@ musoGlue <- function(presCalFile, w, delta = 0.17, settings=setupMuso(), paramet preservedCalib <- preservedCalib[!is.na(preservedCalib$combined),] unfilteredLikelihood <- preservedCalib$combined top5points <- preservedCalib$combined>quantile(preservedCalib$combined,0.95) - preservedCalibtop5 <- preservedCalib[,] + preservedCalibtop5 <- tryCatch(preservedCalib[top5points,], error=function(e){ + warning("Too few simulation to calculate the top five percent, please increase the iteration number") + preservedCalib + }) optRanges <-t(apply(preservedCalibtop5,2,function(x) quantile(x,c(0.05,0.5,0.95)))) pdf("dotplot.pdf") if(lg){ diff --git a/RBBGCMuso/R/parametersweep.R b/RBBGCMuso/R/parametersweep.R index 5e95ce4..48d2365 100644 --- a/RBBGCMuso/R/parametersweep.R +++ b/RBBGCMuso/R/parametersweep.R @@ -53,7 +53,7 @@ You can download pandoc from here: 'https://pandoc.org/',\n or Rstudio from here rmdFile <- sprintf("---\ntitle: \"ParameterSweep basic\"\n---\n```{r setup, include=FALSE}\nknitr::opts_chunk$set(echo = TRUE)\n```\n```{r, echo=FALSE}\nsuppressWarnings(library(RBBGCMuso))\n```\n```{r, echo=FALSE}\nparameters <- read.csv(\"parameters.csv\")\n```\n```{r,fig.width=10, fig.height=3, echo=FALSE}\nnumPar\nfor(i in 1:numPar){\n suppressWarnings(musoQuickEffect(calibrationPar=parameters[i,2],startVal = parameters[i,3], endVal = parameters[i,4],\nnSteps = 9,\noutVar = \"daily_gpp\",\nparName = parameters[i,1],fixAlloc=%s))\n}\n```",fixAlloc) rmdVec <- unlist(strsplit(rmdFile,"\n")) - rmdVec[11] <- paste0("parameters <- read.csv(\"",parameters,"\", stringsAsFactor = FALSE)") + # rmdVec[11] <- paste0("parameters <- read.csv(\"",parameters,"\", stringsAsFactor = FALSE)") rmdVec[14] <- "numPar <- nrow(parameters)" rmdVec[17] <- paste0("nSteps = ", iterations - 1,",") rmdVec[18] <- paste0("outVar = \"",varNames,"\",")