hollorol/RBBGCMuso
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

constrained calibration

Browse Source
This commit is contained in:
Hollos Roland 2021-06-02 14:48:36 +02:00
parent b585f7eb7a
commit 18595570d3
1 changed files with 77 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
RBBGCMuso/R/multiSite.R
View File

@ -6,17 +6,47 @@ annualAggregate <- function(x, aggFun){
tapply(x, rep(1:(length(x)/365), each=365), aggFun) tapply(x, rep(1:(length(x)/365), each=365), aggFun)
} }
SELECT <- function(x,selectPart){ SELECT <- function(x, selectPart){
index <- as.numeric(selectPart) if(!is.function(selectPart)){
if(!is.na(index)){ index <- as.numeric(selectPart)
tapply(x,rep(1:(length(x)/365),each=365), function(y){ tapply(x,rep(1:(length(x)/365),each=365), function(y){
y[index] y[index]
}) })
} else { } else {
tapply(x,rep(1:(length(x)/365),each=365), match.fun(selectPart)) tapply(x,rep(1:(length(x)/365),each=365), selectPart)
} }
} }
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){
with(as.data.frame(mod), {
nexpr <- nrow(constrTable)
filtered <- numeric(nexpr)
for(i in 1:nexpr){
val <- compose(constrTable[i,1])
filtered[i] <- (val <= constrTable[i,3]) &&
(val >= constrTable[i,2])
}
filtered
})
}
modCont <- function(expr, datf, interval, dumping_factor){ modCont <- function(expr, datf, interval, dumping_factor){
tryCatch({ tryCatch({
@ -62,7 +92,6 @@ copyToThreadDirs2 <- function(iniSource, thread_prefix = "thread", numCores, exe
#' @export #' @export
multiSiteCalib <- function(measurements, multiSiteCalib <- function(measurements,
calTable, calTable,
constTable,
parameters, parameters,
dataVar, dataVar,
iterations = 100, iterations = 100,
@ -73,7 +102,8 @@ multiSiteCalib <- function(measurements,
numCores = (parallel::detectCores()-1), numCores = (parallel::detectCores()-1),
pb = txtProgressBar(min=0, max=iterations, style=3), pb = txtProgressBar(min=0, max=iterations, style=3),
pbUpdate = setTxtProgressBar, pbUpdate = setTxtProgressBar,
copyThread = TRUE copyThread = TRUE,
constraints=NULL, th = 10
){ ){
future::plan(future::multisession) future::plan(future::multisession)
@ -102,17 +132,18 @@ multiSiteCalib <- function(measurements,
threadCount <- distributeCores(iterations, numCores) threadCount <- distributeCores(iterations, numCores)
fut <- lapply(1:numCores, function(i) { fut <- lapply(1:numCores, function(i) {
# browser() # browser()
# future({ future({
# tryCatch( tryCatch(
multiSiteThread(measuredData = measurements, parameters = parameters, calTable=calTable, multiSiteThread(measuredData = measurements, parameters = parameters, calTable=calTable,
dataVar = dataVar, iterations = threadCount[i], dataVar = dataVar, iterations = threadCount[i],
likelihood = likelihood, threadNumber= i, burnin=burnin) likelihood = likelihood, threadNumber= i, burnin=burnin,constraints=constraints, th=th)
browser() # setwd("../")
# , error = function(e){ , error = function(e){
# writeLines(as.character(iterations),"progress.txt") saveRDS(e,"error.RDS")
# }) writeLines(as.character(iterations),"progress.txt")
# }) })
})
}) })
# _ _ # _ _
@ -160,7 +191,6 @@ multiSiteCalib <- function(measurements,
# | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \ # | | / _ \| '_ ` _ \| '_ \| | '_ \ / _ \
# | |__| (_) | | | | | | |_) | | | | | __/ # | |__| (_) | | | | | | |_) | | | | | __/
# \____\___/|_| |_| |_|_.__/|_|_| |_|\___| # \____\___/|_| |_| |_|_.__/|_|_| |_|\___|
resultFiles <- list.files(pattern="preservedCalib.*csv$",recursive=TRUE) resultFiles <- list.files(pattern="preservedCalib.*csv$",recursive=TRUE)
res0 <- read.csv(grep("thread_1/",resultFiles, value=TRUE),stringsAsFactors=FALSE) res0 <- read.csv(grep("thread_1/",resultFiles, value=TRUE),stringsAsFactors=FALSE)
resultFilesSans0 <- grep("thread_1/", resultFiles, value=TRUE, invert=TRUE) resultFilesSans0 <- grep("thread_1/", resultFiles, value=TRUE, invert=TRUE)
@ -173,8 +203,12 @@ multiSiteCalib <- function(measurements,
calibrationPar <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["calibrationPar"]] calibrationPar <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["calibrationPar"]]
origModOut <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["origModOut"]] origModOut <- future::value(fut[[1]], stdout = FALSE, signal=FALSE)[["origModOut"]]
# Just single objective version TODO:Multiobjective # Just single objective version TODO:Multiobjective
bestCase <- which.max(results[,ncol(results)-1]) results <- results[results[,"Const"] == 1,]
parameters <- results[bestCase,1:(ncol(results)-2)] # the last two column is the (log) likelihood and the rmse if(nrow(results)==0){
stop("No simulation suitable for constraints")
}
bestCase <- which.max(results[,ncol(results)-2])
parameters <- results[bestCase,1:(ncol(results)-3)] # 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 #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] firstDir <- list.dirs("tmp/thread_1",full.names=TRUE,recursive =FALSE)[1]
@ -187,7 +221,7 @@ multiSiteCalib <- function(measurements,
}) })
setwd("tmp/thread_1") setwd("tmp/thread_1")
aposteriori<- spatialRun(settingsProto, calibrationPar, parameters, calTable) aposteriori<- spatialRun(settingsProto, calibrationPar, parameters, calTable)
file.copy(list.files(list.dirs(full.names=TRUE, recursive=FALSE)[1],pattern=".*\\.epc", full.names=TRUE), file.copy(list.files(list.dirs(full.names=TRUE, recursive=FALSE)[1], pattern=".*\\.epc", full.names=TRUE),
"../../multiSiteOptim.epc", overwrite=TRUE) "../../multiSiteOptim.epc", overwrite=TRUE)
setwd("../../") setwd("../../")
#TODO: Have to put that before multiSiteThread, we should not have to calculate it at every iterations #TODO: Have to put that before multiSiteThread, we should not have to calculate it at every iterations
@ -196,12 +230,12 @@ multiSiteCalib <- function(measurements,
res <- list() res <- list()
res[["calibrationPar"]] <- calibrationPar res[["calibrationPar"]] <- calibrationPar
res[["parameters"]] <- parameters res[["parameters"]] <- parameters
res[["comparison"]] <- compareCalibratedWithOriginal(key="grainDM", modOld=origModOut, modNew=aposteriori, mes=measurements, res[["comparison"]] <- compareCalibratedWithOriginal(key = "grainDM", modOld=origModOut, modNew=aposteriori, mes=measurements,
likelihoods=likelihood, likelihoods = likelihood,
alignIndexes=alignIndexes, alignIndexes = alignIndexes,
musoCodeToIndex = musoCodeToIndex, musoCodeToIndex = musoCodeToIndex,
nameGroupTable = nameGroupTable, mean) nameGroupTable = nameGroupTable, mean)
res[["likelihood"]] <- results[bestCase,ncol(results)-1] res[["likelihood"]] <- results[bestCase,ncol(results)-2]
comp <- res$comparison comp <- res$comparison
res[["originalMAE"]] <-mean(abs((comp[,1]-comp[,3]))) res[["originalMAE"]] <-mean(abs((comp[,1]-comp[,3])))
res[["MAE"]] <- mean(abs((comp[,2]-comp[,3]))) res[["MAE"]] <- mean(abs((comp[,2]-comp[,3])))
@ -247,7 +281,7 @@ multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
outVars = NULL, iterations = 300, outVars = NULL, iterations = 300,
skipSpinup = TRUE, plotName = "calib.jpg", skipSpinup = TRUE, plotName = "calib.jpg",
modifyOriginal=TRUE, likelihood, uncertainity = NULL, burnin=NULL, modifyOriginal=TRUE, likelihood, uncertainity = NULL, burnin=NULL,
naVal = NULL, postProcString = NULL, threadNumber) { naVal = NULL, postProcString = NULL, threadNumber, constraints=NULL,th=10) {
originalRun <- list() originalRun <- list()
nameGroupTable <- calTable nameGroupTable <- calTable
@ -290,12 +324,12 @@ multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations, burnin = burnin) randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations, burnin = burnin)
origEpc <- readValuesFromFile(epcFile, randVals[[1]]) origEpc <- readValuesFromFile(epcFile, randVals[[1]])
partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar)) partialResult <- matrix(ncol=length(randVals[[1]])+2*length(dataVar) + 1)
colN <- randVals[[1]] colN <- randVals[[1]]
colN[match(parameters[,2],randVals[[1]])] <- parameters[,1] colN[match(parameters[,2],randVals[[1]])] <- parameters[,1]
colN[match(parameters[,2], randVals[[1]])[!is.na(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)), colnames(partialResult) <- c(colN,sprintf("%s_likelihood",names(dataVar)),
sprintf("%s_rmse",names(dataVar))) sprintf("%s_rmse",names(dataVar)),"Const")
numParameters <- length(colN) numParameters <- length(colN)
partialResult[1:numParameters] <- origEpc partialResult[1:numParameters] <- origEpc
## Prepare the preservedCalib matrix for the faster ## Prepare the preservedCalib matrix for the faster
@ -338,7 +372,7 @@ multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
mes=measuredData, mes=measuredData,
likelihoods=likelihood, likelihoods=likelihood,
alignIndexes=alignIndexes, alignIndexes=alignIndexes,
musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, mean) musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, groupFun=mean, constraints=constraints,th=th)
write.csv(x=randVals[[1]],"../randIndexes.csv") write.csv(x=randVals[[1]],"../randIndexes.csv")
write.csv(x=partialResult, file="preservedCalib.csv",row.names=FALSE) write.csv(x=partialResult, file="preservedCalib.csv",row.names=FALSE)
@ -369,7 +403,7 @@ multiSiteThread <- function(measuredData, parameters = NULL, startDate = NULL,
mes=measuredData, mes=measuredData,
likelihoods=likelihood, likelihoods=likelihood,
alignIndexes=alignIndexes, alignIndexes=alignIndexes,
musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, mean) musoCodeToIndex = musoCodeToIndex,nameGroupTable = nameGroupTable, groupFun=mean, constraints = constraints, th=th)
@ -405,7 +439,15 @@ prepareFromAgroMo <- function(fName){
calcLikelihoodsForGroups <- function(dataVar, mod, mes, calcLikelihoodsForGroups <- function(dataVar, mod, mes,
likelihoods, alignIndexes, musoCodeToIndex, likelihoods, alignIndexes, musoCodeToIndex,
nameGroupTable, groupFun){ nameGroupTable, groupFun, constraints,
th = 10){
if(!is.null(constraints)){
constRes<- sapply(mod,function(m){
compoVect(m,constraints)
})
}
likelihoodRMSE <- sapply(names(dataVar),function(key){ likelihoodRMSE <- sapply(names(dataVar),function(key){
modelled <- as.vector(unlist(sapply(sort(names(alignIndexes)), modelled <- as.vector(unlist(sapply(sort(names(alignIndexes)),
function(domain_id){ function(domain_id){
@ -426,17 +468,18 @@ calcLikelihoodsForGroups <- function(dataVar, mod, mes,
})) }))
measured <- measured[measured$var_id == key,] measured <- measured[measured$var_id == key,]
# measured <- measured[!is.na(measured)]
res <- c(likelihoods[[key]](modelled, measured), res <- c(likelihoods[[key]](modelled, measured),
sqrt(mean((modelled-measured$mean)^2)) sqrt(mean((modelled-measured$mean)^2))
) )
print(abs(mean(modelled)-mean(measured$mean))) print(abs(mean(modelled)-mean(measured$mean)))
# browser()
res res
}) })
browser()
names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar)) likelihoodRMSE <- c(likelihoodRMSE[1,],likelihoodRMSE[2,],
return(c(likelihoodRMSE[1,],likelihoodRMSE[2,])) ifelse((100 * sum(apply(constRes,2,prod)) / ncol(constRes)) >= th,
1,0))
names(likelihoodRMSE) <- c(sprintf("%s_likelihood",dataVar), sprintf("%s_rmse",dataVar), "Const")
return(likelihoodRMSE)
} }
commonIndexes <- function (settings,measuredData) { commonIndexes <- function (settings,measuredData) {