CIRM for article

2022-04-04 16:08:49 +02:00
parent 7017705f43
commit 7d6d81fa96
7 changed files with 831 additions and 0 deletions
@@ -0,0 +1,72 @@
+This is a si
+
+## Preparations
+
+### Loading the RBBGCMuso package and the necessary functions
+
+```{r}
+library(RBBGCMuso)
+source("make_individual_trees.R") # The DT creation and update algorithms
+source("glue.R") # GLUE optimizer algorithms
+```
+
+The file containing the path to the observation files (Martonvasar_maize.obs), and the parameter intervals (Martonj) 
+
+### Reading the observations
+
+The mean yield had to be adjust. see in art.
+
+```{r}
+measureFile <- "Martonvasar_maize.obs"
+measurements <- read.csv2(measureFile, stringsAsFactors=FALSE)
+measurements$mean <- measurements$mean / 10000 * 0.85
+measurements$sd <- measurements$sd / 10000 * 0.85
+```
+
+### Define conditioning functions
+
+constraints.json
+
+```{json}
+{
+    "constraints": [
+
+        {
+            "Expression": "SELECT(harvest_index, max)|median",
+            "Min": 0.45,
+            "Max": 0.55
+        },
+        {
+            "Expression": "SELECT(proj_lai, max)|quantile(.,0.5)",
+            "Min": 2.7,
+            "Max": 5
+        },
+        {
+            "Expression": "SELECT(rootdepth5, max)|quantile(.,0.5)",
+            "Min": 1.40,
+            "Max": 1.80
+        },			
+        {
+            "Expression": "SELECT(flower_date, max)|quantile(.,0.5)",
+            "Min": 180,
+            "Max": 190 
+        }
+    ],
+
+    "treshold": 80
+}
+```
+
+```{r}
+constraints <- jsonlite::read_json("constraints.json",simplifyVector=TRUE)
+```
+
+### Cal file:
+
+```{verbatim}
+
+Martonvasar_maize.obs
+Martonvasar_maize.set
+site
+Martonvasar_maize;211
+```
@@ -0,0 +1,55 @@
+zero_var <- function(m){
+    apply(m,2, function(v){
+        var(v) != 0
+    })
+
+}
+
+
+glue <- function(results="result.csv",res_r="results.RDS",output ="gplot.pdf",epcname="maize_glue.epc"){
+    res <- read.csv(results)[-1]
+    res <- res[-1,]
+    colnames(res)
+    non_zero <- res[,1:(ncol(res)-4)]
+    colnames(non_zero) <- colnames(res)[1:(ncol(res)-4)]
+    impvars <- zero_var(non_zero)
+    nonzero <- non_zero[,impvars]
+    likelihoods <- res[,(ncol(res)-3)]
+    rmse <- res[,(ncol(res)-2)]
+    const <- res$Const
+    namess <- gsub("__.*","",colnames(nonzero))
+
+
+
+
+    likelihoods <- likelihoods[res$Const==1]
+    goods <- res[res$Const==1,]
+    medlik <- median(likelihoods[likelihoods >= quantile(likelihoods,0.95)])
+    medlik_place <- which.min(abs(likelihoods - medlik))
+    parameters <- readRDS(res_r)
+    glue_opt <- goods[medlik_place, 1:(ncol(res)-4)][impvars]
+    nonka <- goods[likelihoods >= quantile(likelihoods,0.95),1:(ncol(res)-4)]
+    med_opt <- apply(nonka,2,median)[impvars]
+    # med_opt <- apply(nonka,2,mean)[impvars]
+    ml_opt <- goods[which.max(likelihoods),1:(ncol(res)-4)][impvars]
+
+    calibrationPar <- parameters$calibrationPar[impvars]
+    changemulline(src="maize.epc", calibrationPar = calibrationPar, contents=glue_opt, outFiles = epcname)
+    changemulline(src="maize.epc", calibrationPar = calibrationPar, contents=med_opt, outFiles = "maize_median.epc")
+    changemulline(src="maize.epc", calibrationPar = calibrationPar, contents=ml_opt, outFiles = "maize_ml.epc")
+
+
+
+    print(output)
+    pdf(output)
+    for(i in 1:ncol(nonzero)){
+        plot(nonzero[,i],res[,(ncol(res)-3)],main="",col="lightgray", pch=20, cex=0.4, xlab=namess[i],ylab="logLikelihood")
+        points(nonzero[const==1,i],res[const==1,(ncol(res)-3)],pch=20, cex=0.6, col="red",type="p",xlab=namess[i],ylab="logLikelihood")
+        abline(v=glue_opt[i],col="green")
+        abline(v=med_opt[i],col="blue")
+        abline(v=ml_opt[i],col="black")
+    }
+
+    dev.off()
+}
+
@@ -0,0 +1,45 @@
+start_intervals <- read.csv("1/Martonvasar_maize.set",skip=1,stringsAsFactors=FALSE)
+
+indices <- which(start_intervals[,3] != start_intervals[,4])
+
+png("kichen_sink.png",width=30,height=30,res=600,units = "cm")
+
+par(mfrow=c(5,4))
+
+for(i in indices){
+    ranges <- start_intervals[i,3:4]
+    optimes <- numeric(10)
+    for(j in 1:10){
+        base_table <- read.csv(paste(j,"Martonvasar_maize_after_tree.set",sep="/"),
+                               skip=1, stringsAsFactors=FALSE)
+        ranges <- rbind(ranges,base_table[i,3:4])
+        optimes[j] <- unlist(readRDS(paste0(j,"/results.RDS"))$parameters[start_intervals[indices,1]][indices==i])
+    }
+    plot(ranges[,1],11:1,type="l",xlim=range(ranges),main=base_table[i,1],xlab="",ylab="iterations",yaxt="n")
+    axis(2,at=11:1,labels = 0:10)
+    points(optimes,10:1)
+    lines(ranges[,2],11:1,type="l")
+}
+dev.off()
+
+postscript("kichen_sink.eps",paper="a4")
+
+par(mfrow=c(5,4))
+
+for(i in indices){
+    ranges <- start_intervals[i,3:4]
+    optimes <- numeric(10)
+    for(j in 1:10){
+        base_table <- read.csv(paste(j,"Martonvasar_maize_after_tree.set",sep="/"),
+                               skip=1, stringsAsFactors=FALSE)
+        ranges <- rbind(ranges,base_table[i,3:4])
+        optimes[j] <- unlist(readRDS(paste0(j,"/results.RDS"))$parameters[start_intervals[indices,1]][indices==i])
+    }
+    plot(ranges[,1],11:1,type="l",xlim=range(ranges),main=base_table[i,1],xlab="",ylab="iterations",yaxt="n")
+    axis(2,at=11:1,labels = 0:10)
+    points(optimes,10:1)
+    lines(ranges[,2],11:1,type="l")
+}
+dev.off()
+
+
@@ -0,0 +1,131 @@
+library(rpart)
+library(rpart.plot)
+zero_var <- function(m){
+    apply(m,2, function(v){
+        var(v) != 0
+    })
+
+}
+
+decbin <- function(decnum){
+    if(decnum < 2){
+        return(decnum)
+    }
+    c(decbin((decnum %/% 2)),decnum %% 2)
+}
+
+decpad <- function(decnum,len){
+    binrep <- decbin(decnum)
+    c(rep(0,len-length(binrep)),binrep)
+}
+
+
+tree_per_const <- function(results="result.csv",output ="tree_per_const.pdf",
+                           parameters_file="Martonvasar_maize.set"){
+    varname <-readLines(parameters_file)[1]
+    parameters <- read.csv(parameters_file,skip=1,stringsAsFactors=FALSE)
+
+
+    results <- read.csv(results, stringsAsFactors=FALSE)
+    # likelihoods <- results[,ncol(results)-3]
+    # results <- results[likelihoods>=quantile(likelihoods,0.95),]
+    len <- round(log(max(results$failType),2))
+    failTypes <- do.call(rbind,lapply(results$failType,function(x){decpad(x,len)}))
+    pdf(output)
+    sapply(1:len, function(const){
+        nonzero <- results[,1:(ncol(results)-4)]
+        nonzero <- nonzero[,-1]
+        nonzero <- nonzero[,zero_var(nonzero)]
+        colnames(nonzero) <- gsub("__.*","",colnames(nonzero))
+        constraint <- failTypes[,const]
+        baseTable <- cbind.data.frame(nonzero,constraint = as.factor(constraint))
+        try({
+            rp = rpart(constraint ~ .,data = baseTable)
+        })
+        try({
+            parameters <<- update_parameters_based_on_tree(rp, parameters)
+        })
+        try({
+            rpart.plot(rp)
+        })
+
+
+    })
+    dev.off()
+    outname <- paste0(tools::file_path_sans_ext(parameters_file),"_after_tree.",
+                      tools::file_ext(parameters_file))
+    writeLines(varname,outname)
+    write.table(parameters,outname,row.names=FALSE,append=TRUE,sep=",",quote=FALSE)
+}
+update_parameters_based_on_tree <- function(rp, parameters){
+    frm <- rp$frame
+    nodes <- labels(rp)
+    names(nodes) <- row.names(frm)
+    node <- get_start_node(frm)
+    parameters <- change_parameters_on_node(nodes,node,parameters)
+    while(node !=1){
+       node <- get_parent_node(node) 
+        if(node == 1){
+            break()
+        }
+       parameters <- change_parameters_on_node(nodes,node,parameters)
+    }
+    parameters
+}
+
+
+
+
+parse_rule_row <- function(string){
+
+  rule_row  <- regmatches(string,regexec("([a-zA-Z_0-9]+)([>=< ]+)(.*)",string,perl=TRUE))[[1]][-1]
+    if(rule_row[2] == ">="){
+        rule_num <- 1
+    } else {
+        rule_num <- 2
+    }
+    rule <- list(c(rule_num,as.numeric(rule_row[3])))
+    names(rule) <- rule_row[1]
+    return(rule)
+}
+
+get_start_node <- function(frm){
+    nfrm <- frm[frm$yval == 2,]
+    nfrm <- nfrm[nfrm[,"var"] == "<leaf>",]
+    pot_start <- as.numeric(row.names(nfrm))[which.max(nfrm$n)]
+    pot_start
+}
+
+get_parent_node <- function(node_id){
+    as.integer(node_id/2)
+}
+
+change_parameters_on_node <- function(nodes,node,parameters2){
+    crule <- parse_rule_row(nodes[as.character(node)])
+    minmax <- unlist(parameters2[parameters2[,1] == names(crule),c(3,4)])
+
+    if(crule[[1]][1] == 1){
+        if(minmax[1]<=crule[[1]][2]){
+            minmax[1] <- crule[[1]][2]
+            if(minmax[1] <= minmax[2]){
+                parameters2[parameters2[,1] == names(crule),c(3,4)] <- minmax
+            } else {
+                write(sprintf("WARNING: %s's  minimum(%s) > maximum(%s)", parameters2[,1],
+                                                   minmax[1], minmax[2]),  "errorlog.txt", append=TRUE)
+            }
+        }
+    } else {
+        if(minmax[2]>=crule[[1]][2]){
+            minmax[2] <- crule[[1]][2]
+            if(minmax[1] <= minmax[2]){
+                parameters2[parameters2[,1] == names(crule),c(3,4)] <- minmax
+            } else {
+                write(sprintf("WARNING: %s's  minimum(%s) > maximum(%s)", parameters2[,1],
+                                                   minmax[1], minmax[2]),  "errorlog.txt", append=TRUE)
+            }
+     }
+    }
+    parameters2
+}
+
+
@@ -0,0 +1,76 @@
+library(RBBGCMuso)
+file.copy("../../start_set/maize.epc","./start_set/maize.epc",overwrite=TRUE)
+setwd("start_set/")
+
+rmse  <- function(modelled, measured){
+    sqrt(mean((modelled-measured)**2))
+}
+
+r2 <- function(modelled, measured){
+    summary(lm("mod ~ meas", data= data.frame(mod=modelled,meas=measured)))$r.squared
+}
+
+modeff <- function(modelled, measured){
+    1 - (sum((modelled-measured)**2) / sum((measured - mean(measured))**2))
+}
+
+bias <- function(modelled, measured){
+    mean(modelled) - mean(measured)
+}
+
+get_stats <- function(modelled,measured){
+    c(r2=r2(modelled,measured),
+      rmse=rmse(modelled,measured),
+      bias=bias(modelled,measured),
+      modeff=modeff(modelled,measured))
+}
+
+get_modelled <- function(obsTable,settings, ...){
+    simulation <- runMuso(settings, ...)
+    yield <- simulation[,"fruit_DM"]
+    modelled <- yield[match(as.Date(obsTable$date),as.Date(names(yield),"%d.%m.%Y"))] * 10
+    modelled 
+}
+
+obsTable <- read.csv2("Martonvasar_maize.obs",stringsAsFactors=FALSE)
+obsTable$mean <- obsTable$mean / 1000 * 0.85
+obsTable$sd <- obsTable$sd / 1000
+measured <- obsTable$mean
+
+# apriori
+settings <- setupMuso(iniInput=c("n.ini","n.ini"))
+modelled <- get_modelled(obsTable, settings)
+results <- matrix(ncol=4,nrow=11)
+colnames(results) <- c("r2","rmse","bias","modeff")
+row.names(results) <- 0:10
+results[1,] <- get_stats(modelled,measured)
+# max_likelihood_stats
+for(i in 1:10){
+    file.copy(sprintf("../%s/maize_ml.epc",i),"maize.epc",overwrite=TRUE)
+    settings <- setupMuso(iniInput=c("n.ini","n.ini"))
+    modelled <- get_modelled(obsTable, settings)
+    results[i+1,] <- get_stats(modelled, measured)
+}
+# median_stats
+results_med <- results
+
+for(i in 1:10){
+    file.copy(sprintf("../maize_median_step%02d_corrected.epc",i),"maize.epc",overwrite=TRUE)
+    settings <- setupMuso(iniInput=c("n.ini","n.ini"))
+    modelled <- get_modelled(obsTable, settings)
+    results_med[i+1,] <- get_stats(modelled, measured)
+}
+results_med
+
+succes_ratio <- numeric(10)
+for(i in 1:10){
+    succes_ratio[i] <- sum(read.csv(sprintf("../%s/result.csv",i),stringsAsFactors=FALSE)$Const[-1])/10000
+}
+names(succes_ratio) <- 1:10
+png("../success_rate.png",height=30,width=30,res=300, units="cm")
+barplot(succes_ratio,ylim=c(0,1),xlab="iteration number",ylab="Succes rate")
+dev.off()
+
+postscript("../success_rate.eps")
+barplot(succes_ratio,ylim=c(0,1),xlab="iteration number",ylab="Succes rate")
+dev.off()
@@ -0,0 +1,70 @@
+library(rpart)
+library(rpart.plot)
+
+accuracy <- function(x,rp){
+    # Accuracy = (TP + TN)/(TP + TN + FP + FP)
+    # TP: True Positive
+    # TN: True Negative
+    # FP: False Positive
+    # FN: False Negative
+
+    predicted <- rpart.predict(rp,type = "vector")
+    predicted[predicted==1] <- 0
+    predicted[predicted==2] <- 1
+    (sum(x*predicted) +  sum((x + predicted) == 0)) / length(x)
+}
+
+zero_var <- function(m){
+    apply(m,2, function(v){
+        var(v) != 0
+    })
+
+}
+
+decbin <- function(decnum){
+    if(decnum < 2){
+        return(decnum)
+    }
+    c(decbin((decnum %/% 2)),decnum %% 2)
+}
+
+decpad <- function(decnum,len){
+    binrep <- decbin(decnum)
+    c(rep(0,len-length(binrep)),binrep)
+}
+
+
+tree_per_const <- function(results="result.csv",output ="tree_per_const.pdf",
+                           parameters_file="Martonvasar_maize.set"){
+    varname <-readLines(parameters_file)[1]
+    parameters <- read.csv(parameters_file,skip=1,stringsAsFactors=FALSE)
+
+
+    results <- read.csv(results, stringsAsFactors=FALSE)
+    # likelihoods <- results[,ncol(results)-3]
+    # results <- results[likelihoods>=quantile(likelihoods,0.95),]
+    len <- round(log(max(results$failType),2))
+    failTypes <- do.call(rbind,lapply(results$failType,function(x){decpad(x,len)}))
+    sapply(1:len, function(const){
+        nonzero <- results[,1:(ncol(results)-4)]
+        nonzero <- nonzero[,-1]
+        nonzero <- nonzero[,zero_var(nonzero)]
+        colnames(nonzero) <- gsub("__.*","",colnames(nonzero))
+        constraint <- failTypes[,const]
+        baseTable <- cbind.data.frame(nonzero,constraint = as.factor(constraint))
+        tryCatch({
+            rp <- rpart(constraint ~ .,data = baseTable)
+            accuracy(constraint, rp)
+        }, error = function(e){NA})
+    })
+}
+
+results <- matrix(nrow=10,ncol=4)
+row.names(results) <- 1:10
+colnames(results) <- c("Harvest Index", "LAI", "Root depth in phen. 5", "Flowering date")
+for(i in 1:10){
+    setwd(as.character(i))
+    results[i,] <- tree_per_const(parameters_file="Martonvasar_maize_after_tree.set")
+    setwd("../")
+}
+results