From 0d3666a8f5ad1468cfa0ce34f5666d166c053356 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roland=20Holl=C3=B3s?= Date: Wed, 14 Jun 2023 17:53:28 +0200 Subject: [PATCH 1/5] solving the old dependency group problem --- RBBGCMuso/R/assistantFunctions.R | 1 + RBBGCMuso/R/musoRand.R | 229 ++++++++++++++++++------------- 2 files changed, 138 insertions(+), 92 deletions(-) diff --git a/RBBGCMuso/R/assistantFunctions.R b/RBBGCMuso/R/assistantFunctions.R index 7d442f9..ad32f2a 100644 --- a/RBBGCMuso/R/assistantFunctions.R +++ b/RBBGCMuso/R/assistantFunctions.R @@ -150,6 +150,7 @@ readValuesFromFile <- function(filename, linums){ return(values) } + #' readMeasuredMuso #' #' MuSo data reader diff --git a/RBBGCMuso/R/musoRand.R b/RBBGCMuso/R/musoRand.R index 1ed8c3a..9148af4 100644 --- a/RBBGCMuso/R/musoRand.R +++ b/RBBGCMuso/R/musoRand.R @@ -8,147 +8,182 @@ #' @importFrom limSolve xsample #' @export -musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = NULL, burnin = NULL){ +musoRand <- function(parameters, iterations=3000, fileType="epc", sourceFile=NULL, constrains = NULL, burnin = NULL){ + if(is.null(constrains)){ constMatrix <- constrains constMatrix <- getOption("RMuso_constMatrix")[[fileType]][[as.character(getOption("RMuso_version"))]] } else { constMatrix <- constrains } - + parameters <- parameters[,-1] constMatrix <- constMatrix[,-1] - + depTableMaker <- function(constMatrix,parameters){ - # browser() parameters <- parameters[order(parameters[,1]),] ## BUG!!! selectedRows <- constMatrix[,"INDEX"] %in% parameters[,1] + # constMatrix[constMatrix[,"INDEX"] %in% parameters[,1],] rankList <- rank(constMatrix[selectedRows,2]) constMatrix[selectedRows,c(5,6)] <- parameters[rankList,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 + + + + + + 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])})),] + + + + paragroups <- unique(constMatrix[constMatrix[,"GROUP"] != 0,"GROUP"]) + missingMembers <- list() + for(group in paragroups){ + groupMembers <- constMatrix[constMatrix[,"GROUP"] == group,"INDEX"] + missingMemberElems <- groupMembers[is.na(match(groupMembers,parameters[,1]))] + if(length(missingMemberElems) > 0){ + missingMembers$indices <- c(missingMembers$indices, + match(missingMemberElems,constMatrix[,"INDEX"])) + if(is.null(sourceFile)){ + stop(sprintf("All group members of the group (%s) have to be in parameters if sourceFile not(epc file)",group)) + } + + missingMembers$values <- c(missingMembers$values, + suppressWarnings(readValuesFromFile(sourceFile,missingMemberElems))) + } + } + + list(depTable=constMatrix,missingMembers=missingMembers) } # browser() - 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)) + genMat0 <- function(dep, missingIndices){ + numberOfVariable <- nrow(dep) + if(length(missingMembers) != 0){ + dep <- dep[-missingIndices,] + } + 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"])) + ## 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) + numberOfVariables <- nrow(dep) + G<- -1*diag(numberOfVariables) - for(i in 1:numberOfVariables){ - if(dep[i,4]!=0){ - G[i,dep[i,4]] <- 1 - } + for(i in 1:numberOfVariables){ + if(dep[i,4]!=0){ + G[i,dep[i,4]] <- 1 + } + + } + # browser() + G<-G[dep[,4]!=0,] - } -# browser() - G<-G[dep[,4]!=0,] - if(is.null(nrow(G))){ G<-t(as.matrix(G)) } numRowsInG <- nrow(G) - if(Range[1]==1){ - G<-cbind(G,matrix(ncol=(N-Range[2]),nrow=numRowsInG,data=0)) - } else{ - if(Range[2]==N){ - G<-cbind(matrix(ncol=(Range[1]-1),nrow=numRowsInG,data=0),G) - } else { - G <- cbind(matrix(ncol=(Range[1]-1),nrow=numRowsInG,data=0),G,matrix(ncol=(N-Range[2]),nrow=numRowsInG,data=0)) - } - } - return(list(G=-1*G,h=-1*rep(0,nrow(G)))) + if(Range[1]==1){ + G<-cbind(G,matrix(ncol=(N-Range[2]),nrow=numRowsInG,data=0)) + } else{ + if(Range[2]==N){ + G<-cbind(matrix(ncol=(Range[1]-1),nrow=numRowsInG,data=0),G) + } else { + G <- cbind(matrix(ncol=(Range[1]-1), + nrow=numRowsInG,data=0), + G, + matrix(ncol=(N-Range[2]), + nrow=numRowsInG,data=0)) + } + } + return(list(G=-1*G,h=-1*rep(0,nrow(G)))) } genMat2 <- function(dep, N){ - G <- rep(1,nrow(dep)) + G <- rep(1,nrow(dep)) - Range <- (function(x){ - c(min(x), max(x)) - })(as.numeric(dep[,"rowIndex"])) + 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])) - } - } + 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]) + G <- t(matrix(sign(dep[2,4])*G)) + h <- abs(dep[1,4]) if(dep[1,"TYPE"]==2){ # This is not needed, I'll have to remove the if part, and keep the content G <- G*(-1) h <- h*(-1) } - return(list(G=G,h=h)) + return(list(G=G,h=h)) } genMat3 <- function(dep, N){ - Range <- (function(x){ - c(min(x), max(x)) - })(as.numeric(dep[,"rowIndex"])) + Range <- (function(x){ + c(min(x), max(x)) + })(as.numeric(dep[,"rowIndex"])) - E <- rep(1,nrow(dep)) + 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])) - } - } + 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)) + 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) + 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) + 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) + missingMembers <- dependences$missingMembers + dependences <- dependences$depTable dependences <- cbind(dependences,1:nrow(dependences)) colnames(dependences)[ncol(dependences)] <- "rowIndex" # browser() @@ -167,21 +202,31 @@ musoRand <- function(parameters, iterations=3000, fileType="epc", constrains = N Ef[[i]] <- applyRandTypeE(splitedDeps[[i]],nrow(dependences)) } } - - Gh0<- genMat0(dependences) + + Gh0<- genMat0(dependences, missingMembers$indices) 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})) - # browser() + # browser() + if(length(missingMembers$indices)!=0){ + Ep <- matrix(data=0,ncol=numberOfVariable,nrow=length(missingMembers$indices)) + Ep[1:length(missingMembers$indices),missingMembers$indices] <- 1 + E <- rbind(E,Ep) + f <- c(f,missingMembers$values) + } + randVal <- suppressWarnings(limSolve::xsample(G=G,H=h,E=E,F=f,burninlength=burnin, iter = iterations))$X } else{ - Gh0<-genMat0(dependences) + Gh0<-genMat0(dependences,NULL) randVal <- suppressWarnings(xsample(G=Gh0$G,H=Gh0$h, iter = iterations))$X } - + if(length(missingMembers$indices)!=0){ + return(list(INDEX = dependences$INDEX[-missingMembers$indices], + randVal=randVal[,-missingMembers$indices])) + } results <- list(INDEX =dependences$INDEX, randVal=randVal) return(results) } From 216f2109fe7f0d2b285725c5c9e10765a50e2609 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roland=20Holl=C3=B3s?= Date: Fri, 16 Jun 2023 09:05:20 +0200 Subject: [PATCH 2/5] delete test.R file --- test.R | 382 --------------------------------------------------------- 1 file changed, 382 deletions(-) delete mode 100644 test.R diff --git a/test.R b/test.R deleted file mode 100644 index 88e85f5..0000000 --- a/test.R +++ /dev/null @@ -1,382 +0,0 @@ -parameters <- -getOption("RMuso_constMatrix")[["epc"]][["6"]] - NAME - yearday to start new growth - yearday to end new growth - transfer growth period as fraction of growing season - litterfall as fraction of growing season - base temperature - minimum temperature for growth displayed on current day - optimal1 temperature for growth displayed on current day - optimal2 temperature for growth displayed on current day - maxmimum temperature for growth displayed on current day - minimum temperature for carbon assimilation displayed on current day - optimal1 temperature for carbon assimilation displayed on current day - optimal2 temperature for carbon assimilation displayed on current day - maxmimum temperature for carbon assimilation displayed on current day - annual leaf and fine root turnover fraction - annual live wood turnover fraction - annual fire mortality fraction - whole-plant mortality paramter for vegetation period - C:N of leaves - C:N of leaf litter - C:N of fine roots - C:N of fruit - C:N of softstem - C:N of live wood - C:N of dead wood - dry matter content of leaves - dry matter content of leaf litter - dry matter content of fine roots - dry matter content of fruit - dry matter content of softstem - dry matter content of live wood - dry matter content of dead wood - leaf litter labile proportion - leaf litter cellulose proportion - fine root labile proportion - fine root cellulose proportion - fruit labile proportion - fruit cellulose proportion - softstem labile proportion - softstem cellulose proportion - dead wood cellulose proportion - canopy water interception coefficient - canopy light extinction coefficient - potential radiation use efficiency - radiation parameter1 (Jiang et al.2015) - radiation parameter1 (Jiang et al.2015) - all-sided to projected leaf area ratio - ratio of shaded SLA:sunlit SLA - fraction of leaf N in Rubisco - fraction of leaf N in PeP - maximum stomatal conductance - cuticular conductance - boundary layer conductance - maximum height of plant - stem weight corresponding to maximum height - plant height function shape parameter (slope) - maximum depth of rooting zone - root distribution parameter - root weight corresponding to max root depth - root depth function shape parameter (slope) - root weight to rooth length conversion factor - growth resp per unit of C grown - maintenance respiration in kgC/day per kg of tissue N - theoretical maximum prop. of non-structural and structural carbohydrates - prop. of non-structural carbohydrates available for maintanance resp - symbiotic+asymbiotic fixation of N - time delay for temperature in photosynthesis acclimation - critical VWCratio (prop. to FC-WP) in germination - critical photoslow daylength - slope of relative photoslow development rate - critical vernalization temperature 1 - critical vernalization temperature 2 - critical vernalization temperature 3 - critical vernalization temperature 4 - slope of relative vernalization development rate - required vernalization days (in vernalization development rate) - critical flowering heat stress temperature 1 - critical flowering heat stress temperature 2 - theoretical maximum of flowering thermal stress mortality - VWC ratio to calc. soil moisture limit 1 (prop. to FC-WP) - VWC ratio to calc. soil moisture limit 2 (prop. to SAT-FC) - minimum of soil moisture limit2 multiplicator (full anoxic stress value) - vapor pressure deficit: start of conductance reduction - vapor pressure deficit: complete conductance reduction - maximum senescence mortality coefficient of aboveground plant material - maximum senescence mortality coefficient of belowground plant material - maximum senescence mortality coefficient of non-structured plant material - lower limit extreme high temperature effect on senescence mortality - upper limit extreme high temperature effect on senescence mortality - turnover rate of wilted standing biomass to litter - turnover rate of cut-down non-woody biomass to litter - turnover rate of cut-down woody biomass to litter - drought tolerance parameter (critical value of day since water stress) - effect of soilstress factor on photosynthesis - crit. amount of snow limiting photosyn. - limit1 (under:full constrained) of HEATSUM index - limit2 (above:unconstrained) of HEATSUM index - limit1 (under:full constrained) of TMIN index - limit2 (above:unconstrained) of TMIN index - limit1 (above:full constrained) of VPD index - limit2 (under:unconstrained) of VPD index - limit1 (under:full constrained) of DAYLENGTH index - limit2 (above:unconstrained) of DAYLENGTH index - moving average (to avoid the effects of extreme events) - GSI limit1 (greater that limit -> start of vegper) - GSI limit2 (less that limit -> end of vegper) - length of phenophase (GDD)-0 - leaf ALLOCATION -0 - fine root ALLOCATION-0 - fruit ALLOCATION -0 - soft stem ALLOCATION-0 - live woody stem ALLOCATION -0 - dead woody stem ALLOCATION -0 - live coarse root ALLOCATION-0 - dead coarse root ALLOCATION -0 - canopy average specific leaf area-0 - current growth proportion-0 - maximal lifetime of plant tissue-0 - length of phenophase (GDD)-1 - leaf ALLOCATION -1 - fine root ALLOCATION-1 - fruit ALLOCATION -1 - soft stem ALLOCATION-1 - live woody stem ALLOCATION -1 - dead woody stem ALLOCATION -1 - live coarse root ALLOCATION-1 - dead coarse root ALLOCATION -1 - canopy average specific leaf area-1 - current growth proportion-1 - maximal lifetime of plant tissue-1 - length of phenophase (GDD)-2 - leaf ALLOCATION -2 - fine root ALLOCATION-2 - fruit ALLOCATION -2 - soft stem ALLOCATION-2 - live woody stem ALLOCATION -2 - dead woody stem ALLOCATION -2 - live coarse root ALLOCATION-2 - dead coarse root ALLOCATION -2 - canopy average specific leaf area-2 - current growth proportion-2 - maximal lifetime of plant tissue-2 - length of phenophase (GDD)-3 - leaf ALLOCATION -3 - fine root ALLOCATION-3 - fruit ALLOCATION -3 - soft stem ALLOCATION-3 - live woody stem ALLOCATION -3 - dead woody stem ALLOCATION -3 - live coarse root ALLOCATION-3 - dead coarse root ALLOCATION -3 - canopy average specific leaf area-3 - current growth proportion-3 - maximal lifetime of plant tissue-3 - length of phenophase (GDD)-4 - leaf ALLOCATION -4 - fine root ALLOCATION-4 - fruit ALLOCATION -4 - soft stem ALLOCATION-4 - live woody stem ALLOCATION -4 - dead woody stem ALLOCATION -4 - live coarse root ALLOCATION-4 - dead coarse root ALLOCATION -4 - canopy average specific leaf area-4 - current growth proportion-4 - maximal lifetime of plant tissue-4 - length of phenophase (GDD)-5 - leaf ALLOCATION -5 - fine root ALLOCATION-5 - fruit ALLOCATION -5 - soft stem ALLOCATION-5 - live woody stem ALLOCATION -5 - dead woody stem ALLOCATION -5 - live coarse root ALLOCATION-5 - dead coarse root ALLOCATION -5 - canopy average specific leaf area-5 - current growth proportion-5 - maximal lifetime of plant tissue-5 - length of phenophase (GDD)-6 - leaf ALLOCATION -6 - fine root ALLOCATION-6 - fruit ALLOCATION -6 - soft stem ALLOCATION-6 - live woody stem ALLOCATION -6 - dead woody stem ALLOCATION -6 - live coarse root ALLOCATION-6 - dead coarse root ALLOCATION -6 - canopy average specific leaf area-6 - current growth proportion-6 - maximal lifetime of plant tissue-6 - INDEX UNIT DEPENDENCE MIN MAX GROUP TYPE - 9.00 yday NA 0.00000 364.0000 0 0 - 10.00 yday NA 0.00000 364.0000 0 0 - 11.00 prop NA 0.00000 1.0000 0 0 - 12.00 prop NA 0.00000 1.0000 0 0 - 13.00 Celsius NA 0.00000 12.0000 0 0 - 14.00 Celsius 0 0.00000 10.0000 1 1 - 15.00 Celsius 1 10.00000 20.0000 1 1 - 16.00 Celsius 2 20.00000 40.0000 1 1 - 17.00 Celsius 3 30.00000 50.0000 1 1 - 18.00 Celsius 0 0.00000 10.0000 2 1 - 19.00 Celsius 1 10.00000 20.0000 2 1 - 20.00 Celsius 2 20.00000 40.0000 2 1 - 21.00 Celsius 3 30.00000 50.0000 2 1 - 22.00 1/yr NA 0.10000 0.4000 0 0 - 23.00 1/yr NA 0.50000 1.0000 0 0 - 24.00 1/yr NA 0.00000 1.0000 0 0 - 25.00 1/vegper NA 0.00000 0.5000 0 0 - 26.00 kgC/kgN 0 10.00000 100.0000 3 1 - 27.00 kgC/kgN 1 10.00000 60.0000 3 1 - 28.00 kgC/kgN 1 10.00000 60.0000 3 1 - 29.00 kgC/kgN 1 10.00000 60.0000 3 1 - 30.00 kgC/kgN 1 10.00000 60.0000 3 1 - 31.00 kgC/kgN 0 50.00000 100.0000 4 1 - 32.00 kgC/kgN 1 300.00000 800.0000 4 1 - 33.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 34.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 35.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 36.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 37.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 38.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 39.00 kgC/kgDM NA 0.20000 0.6000 0 0 - 40.00 prop 1 0.10000 0.6000 5 2 - 41.00 prop 1 0.10000 0.6000 5 2 - 42.00 prop 1 0.10000 0.6000 6 2 - 43.00 prop 1 0.10000 0.6000 6 2 - 44.00 prop 1 0.10000 0.6000 7 2 - 45.00 prop 1 0.10000 0.6000 7 2 - 46.00 prop 1 0.10000 0.6000 8 2 - 47.00 prop 1 0.10000 0.6000 8 2 - 48.00 prop NA 0.50000 0.9000 0 0 - 49.00 1/LAI/d NA 0.01000 0.1000 0 0 - 50.00 dimless NA 0.20000 0.8000 0 0 - 51.00 g/MJ NA 2.00000 2.0000 0 0 - 52.00 dimless NA 0.78100 0.7810 0 0 - 53.00 dimless NA -13.59600 -13.5960 0 0 - 54.00 dimless NA 2.00000 2.0000 0 0 - 55.00 dimless NA 2.00000 2.0000 0 0 - 56.00 dimless NA 0.01000 0.2000 0 0 - 57.00 dimless NA 0.04240 0.0424 0 0 - 58.00 m/s NA 0.00100 0.1000 0 0 - 59.00 m/s NA 0.00001 0.0001 0 0 - 60.00 m/s NA 0.01000 0.0900 0 0 - 61.00 m NA 0.10000 10.0000 0 0 - 62.00 kgC NA 0.10000 100.0000 0 0 - 63.00 dimless NA 0.50000 0.5000 0 0 - 64.00 m NA 0.10000 10.0000 0 0 - 65.00 prop NA 3.67000 3.6700 0 0 - 66.00 kgC/m2 NA 0.40000 0.4000 0 0 - 67.00 prop NA 0.50000 0.5000 0 0 - 68.00 m/kg NA 1000.00000 1000.0000 0 0 - 69.00 prop NA 0.10000 0.5000 0 0 - 70.00 kgC/kgN/d NA 0.10000 0.5000 0 0 - 71.00 dimless NA 0.00000 1.0000 0 0 - 72.00 dimless NA 0.00000 1.0000 0 0 - 73.00 kgN/m2/yr NA 0.00000 0.0010 0 0 - 74.00 day NA 0.00000 50.0000 0 0 - 79.00 prop NA 0.00000 1.0000 0 0 - 81.00 hour NA 14.00000 18.0000 0 0 - 82.00 dimless NA 0.00500 0.0050 0 0 - 84.00 Celsius 0 -5.00000 5.0000 9 1 - 85.00 Celsius 1 0.00000 10.0000 9 1 - 86.00 Celsius 2 5.00000 15.0000 9 1 - 87.00 Celsius 3 10.00000 20.0000 9 1 - 88.00 dimless NA 0.04000 0.0400 0 0 - 89.00 dimless NA 30.00000 70.0000 0 0 - 91.00 Celsius 0 30.00000 40.0000 10 1 - 92.00 Celsius 1 30.00000 50.0000 10 1 - 93.00 prop NA 0.00000 0.4000 0 0 - 96.00 prop NA 0.50000 1.0000 0 0 - 97.00 prop NA 0.50000 1.0000 0 0 - 98.00 prop NA 0.00000 1.0000 0 0 - 99.00 Pa NA 500.00000 1500.0000 0 0 -100.00 Pa NA 1500.00000 3500.0000 0 0 -101.00 prop 0 0.00000 0.1000 0 0 -102.00 prop 1 0.00000 0.1000 0 0 -103.00 prop NA 0.00000 0.1000 0 0 -104.00 Celsius NA 30.00000 40.0000 0 0 -105.00 Celsius NA 30.00000 50.0000 0 0 -106.00 prop NA 0.00000 0.1000 0 0 -107.00 prop NA 0.00000 0.1000 0 0 -108.00 prop NA 0.00000 0.1000 0 0 -109.00 n_day NA 0.00000 100.0000 0 0 -110.00 dimless NA 0.00000 1.0000 0 0 -113.00 kg/m2 NA 0.00000 20.0000 0 0 -114.00 Celsius 0 0.00000 50.0000 11 1 -115.00 Celsius 1 0.00000 100.0000 11 1 -116.00 Celsius 0 -5.00000 5.0000 12 1 -117.00 Celsius 1 0.00000 10.0000 12 1 -118.00 Pa 0 2000.00000 600.0000 13 1 -119.00 Pa 1 500.00000 1500.0000 13 1 -120.00 s 0 0.00000 0.0000 14 1 -121.00 s 1 0.00000 0.0000 14 1 -122.00 n_day NA 2.00000 20.0000 0 0 -123.00 dimless NA 0.00000 0.2000 0 0 -124.00 dimless NA 0.00000 0.1000 0 0 -128.60 Celsius NA 0.00000 10000.0000 0 0 -129.60 prop 1 0.00000 1.0000 15 -3 -130.60 prop 1 0.00000 1.0000 15 -3 -131.60 prop 1 0.00000 1.0000 15 -3 -132.60 prop 1 0.00000 1.0000 15 -3 -133.60 prop 1 0.00000 1.0000 15 -3 -134.60 prop 1 0.00000 1.0000 15 -3 -135.60 prop 1 0.00000 1.0000 15 -3 -136.60 prop 1 0.00000 1.0000 15 -3 -137.60 m2/kg NA 0.00000 2.0000 0 0 -138.60 prop NA 0.00000 0.0000 0 0 -139.60 Celsius NA 1.00000 20000.0000 0 0 -128.61 Celsius NA 0.00000 10000.0000 0 0 -129.61 prop 1 0.00000 1.0000 16 -3 -130.61 prop 1 0.00000 1.0000 16 -3 -131.61 prop 1 0.00000 1.0000 16 -3 -132.61 prop 1 0.00000 1.0000 16 -3 -133.61 prop 1 0.00000 1.0000 16 -3 -134.61 prop 1 0.00000 1.0000 16 -3 -135.61 prop 1 0.00000 1.0000 16 -3 -136.61 prop 1 0.00000 1.0000 16 -3 -137.61 m2/kg NA 0.00000 2.0000 0 0 -138.61 prop NA 0.00000 0.0000 0 0 -139.61 Celsius NA 1.00000 20000.0000 0 0 -128.62 Celsius NA 0.00000 10000.0000 0 0 -129.62 prop 1 0.00000 1.0000 17 -3 -130.62 prop 1 0.00000 1.0000 17 -3 -131.62 prop 1 0.00000 1.0000 17 -3 -132.62 prop 1 0.00000 1.0000 17 -3 -133.62 prop 1 0.00000 1.0000 17 -3 -134.62 prop 1 0.00000 1.0000 17 -3 -135.62 prop 1 0.00000 1.0000 17 -3 -136.62 prop 1 0.00000 1.0000 17 -3 -137.62 m2/kg NA 0.00000 2.0000 0 0 -138.62 prop NA 0.00000 0.0000 0 0 -139.62 Celsius NA 1.00000 20000.0000 0 0 -128.63 Celsius NA 0.00000 10000.0000 0 0 -129.63 prop 1 0.00000 1.0000 18 -3 -130.63 prop 1 0.00000 1.0000 18 -3 -131.63 prop 1 0.00000 1.0000 18 -3 -132.63 prop 1 0.00000 1.0000 18 -3 -133.63 prop 1 0.00000 1.0000 18 -3 -134.63 prop 1 0.00000 1.0000 18 -3 -135.63 prop 1 0.00000 1.0000 18 -3 -136.63 prop 1 0.00000 1.0000 18 -3 -137.63 m2/kg NA 0.00000 2.0000 0 0 -138.63 prop NA 0.00000 0.0000 0 0 -139.63 Celsius NA 1.00000 20000.0000 0 0 -128.64 Celsius NA 0.00000 10000.0000 0 0 -129.64 prop 1 0.00000 1.0000 19 -3 -130.64 prop 1 0.00000 1.0000 19 -3 -131.64 prop 1 0.00000 1.0000 19 -3 -132.64 prop 1 0.00000 1.0000 19 -3 -133.64 prop 1 0.00000 1.0000 19 -3 -134.64 prop 1 0.00000 1.0000 19 -3 -135.64 prop 1 0.00000 1.0000 19 -3 -136.64 prop 1 0.00000 1.0000 19 -3 -137.64 m2/kg NA 0.00000 2.0000 0 0 -138.64 prop NA 0.00000 0.0000 0 0 -139.64 Celsius NA 1.00000 20000.0000 0 0 -128.65 Celsius NA 0.00000 10000.0000 0 0 -129.65 prop 1 0.00000 1.0000 20 -3 -130.65 prop 1 0.00000 1.0000 20 -3 -131.65 prop 1 0.00000 1.0000 20 -3 -132.65 prop 1 0.00000 1.0000 20 -3 -133.65 prop 1 0.00000 1.0000 20 -3 -134.65 prop 1 0.00000 1.0000 20 -3 -135.65 prop 1 0.00000 1.0000 20 -3 -136.65 prop 1 0.00000 1.0000 20 -3 -137.65 m2/kg NA 0.00000 2.0000 0 0 -138.65 prop NA 0.00000 0.0000 0 0 -139.65 Celsius NA 1.00000 20000.0000 0 0 -128.66 Celsius NA 0.00000 10000.0000 0 0 -129.66 prop 1 0.00000 1.0000 21 -3 -130.66 prop 1 0.00000 1.0000 21 -3 -131.66 prop 1 0.00000 1.0000 21 -3 -132.66 prop 1 0.00000 1.0000 21 -3 -133.66 prop 1 0.00000 1.0000 21 -3 -134.66 prop 1 0.00000 1.0000 21 -3 -135.66 prop 1 0.00000 1.0000 21 -3 -136.66 prop 1 0.00000 1.0000 21 -3 -137.66 m2/kg NA 0.00000 2.0000 0 0 -138.66 prop NA 0.00000 0.0000 0 0 -139.66 Celsius NA 1.00000 20000.0000 0 0 From b9e829951e4d542aa12e245fc9ad6f773b5ce822 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roland=20Holl=C3=B3s?= Date: Mon, 19 Jun 2023 10:09:54 +0200 Subject: [PATCH 3/5] add randomization feature to sensitivity and calibration with musoMonte --- RBBGCMuso/R/calibrateMuso.R | 11 ++++++++++- RBBGCMuso/R/musoMonte.R | 3 ++- RBBGCMuso/R/musoSensi.R | 2 ++ 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 87a9ee5..81fe884 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -12,6 +12,7 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s skipSpinup = TRUE, plotName = "calib.jpg", modifyOriginal=TRUE, likelihood, uncertainity = NULL, naVal = NULL, postProcString = NULL, + sourceFile=NULL, # bases for musoRand if dependecy group is not fully defined by parameters.csv thread_prefix="thread", numCores = max(c(parallel::detectCores()-1,1)), pb = txtProgressBar(min=0, max=iterations, style=3), maxLikelihoodEpc=TRUE, pbUpdate = setTxtProgressBar, outputLoc="./", method="GLUE",lg = FALSE, w=NULL, ...){ @@ -127,6 +128,14 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s switch(method, "GLUE"={ 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) + epcVals <- results[which.max(liks),1:length(epcIndexes)] + epcPlace <- file.path(dirname(settings$inputFiles),settings$epc)[2] + changemulline(filePaths= epcPlace, epcIndexes, + epcVals, src =epcPlace,# settings$epcInput[2], + outFiles = file.path(outputLoc, "maxLikelihood_epc.epc")) + names(epcVals) <- epcIndexes }, "agromo"={ liks <- results[,sprintf("%s_likelihood",names(likelihood))] @@ -234,7 +243,7 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, ## row numbers print("optiMuso is randomizing the epc parameters now...",quote = FALSE) if(iterations < 3000){ - randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000) + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000,sourceFile=sourceFile) randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] } else { randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations) diff --git a/RBBGCMuso/R/musoMonte.R b/RBBGCMuso/R/musoMonte.R index d94abeb..5ce4756 100644 --- a/RBBGCMuso/R/musoMonte.R +++ b/RBBGCMuso/R/musoMonte.R @@ -18,6 +18,7 @@ musoMonte <- function(settings=NULL, parameters=NULL, + sourceFile=NULL, inputDir = "./", outLoc = "./calib", iterations = 10, @@ -100,7 +101,7 @@ musoMonte <- function(settings=NULL, ##reading the original epc file at the specified ## row numbers if(iterations < 3000){ - randVals <- musoRand(parameters = parameters,fileType="epc", iterations = 3000) + randVals <- musoRand(parameters = parameters,fileType="epc", iterations = 3000,sourceFile=sourceFile) randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] } else { randVals <- musoRand(parameters = parameters,fileType="epc", iterations = iterations) diff --git a/RBBGCMuso/R/musoSensi.R b/RBBGCMuso/R/musoSensi.R index 4f4130b..33dec85 100644 --- a/RBBGCMuso/R/musoSensi.R +++ b/RBBGCMuso/R/musoSensi.R @@ -31,6 +31,7 @@ musoSensi <- function(monteCarloFile = NULL, plotTitle = "Sensitivity", skipSpinup = TRUE, skipZero = TRUE, + sourceFile=NULL, postProcString=NULL, modifyOut=TRUE, dpi=300){ @@ -87,6 +88,7 @@ musoSensi <- function(monteCarloFile = NULL, if(is.null(monteCarloFile)){ M <- musoMonte(parameters = parameters, settings = settings, + sourceFile=NULL, inputDir = inputDir, outLoc = outLoc, iterations = iterations, From 47b9a0eb09cd658c00759f7ee868749a88257d24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roland=20Holl=C3=B3s?= Date: Mon, 19 Jun 2023 10:47:37 +0200 Subject: [PATCH 4/5] calibration bugfix --- RBBGCMuso/R/calibrateMuso.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 81fe884..062bbd6 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -244,9 +244,9 @@ musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, print("optiMuso is randomizing the epc parameters now...",quote = FALSE) if(iterations < 3000){ randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = 3000,sourceFile=sourceFile) - randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] + randVals[[2]]<- randVals[[2]][sample(1:3000,iterations),] # TODO: last not random } else { - randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations) + randVals <- musoRand(parameters = parameters,constrains = NULL, iterations = iterations,sourceFile=sourceFile) } origEpc <- readValuesFromFile(settings$epc[2],randVals[[1]]) From 26db1a93fd89c16547e8a95da9a21524ad90de99 Mon Sep 17 00:00:00 2001 From: Hollos Roland Date: Tue, 20 Jun 2023 15:13:41 +0200 Subject: [PATCH 5/5] fixing calibration bug --- RBBGCMuso/R/calibrateMuso.R | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/RBBGCMuso/R/calibrateMuso.R b/RBBGCMuso/R/calibrateMuso.R index 062bbd6..6b6ef20 100644 --- a/RBBGCMuso/R/calibrateMuso.R +++ b/RBBGCMuso/R/calibrateMuso.R @@ -45,6 +45,7 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s future({ tryCatch( musoSingleThread(measuredData, parameters, startDate, + sourceFile=settings$epc[2], # EPC SPECIFIC endDate, formatString, dataVar, outLoc, preTag, settings, @@ -53,6 +54,7 @@ calibrateMuso <- function(measuredData, parameters =read.csv("parameters.csv", s modifyOriginal, likelihood, uncertainity, naVal, postProcString, i) , error = function(e){ + # browser() writeLines(as.character(iterations),"progress.txt") }) @@ -130,7 +132,12 @@ 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) - epcVals <- results[which.max(liks),1:length(epcIndexes)] + if(ncol(liks) == 1){ + ml_place <- which.max(liks) + } else { + ml_place <- which.max(as.matrix(liks) %*% as.matrix(w)) + } + epcVals <- results[ml_place,1:length(epcIndexes)] epcPlace <- file.path(dirname(settings$inputFiles),settings$epc)[2] changemulline(filePaths= epcPlace, epcIndexes, epcVals, src =epcPlace,# settings$epcInput[2], @@ -183,6 +190,7 @@ copyToThreadDirs <- function(prefix="thread", numcores=parallel::detectCores()-1 } musoSingleThread <- function(measuredData, parameters = NULL, startDate = NULL, + sourceFile=NULL, endDate = NULL, formatString = "%Y-%m-%d", dataVar, outLoc = "./calib", preTag = "cal-", settings = setupMuso(),