#  Copyright (C) 2011-2012 Sandra Ligges
#
#  This program is free software; you can redistribute it and/or modify it under the terms
#  of the GNU General Public License as published by the Free Software Foundation;
#  either version 2 of the License, or (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
#  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#  See the GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available in the same folder as this file.


#######################################################################################################################
### Funktion zur Erstellung der Balkendiagramme für die Ergebnisse der Punktschätzverfahren
#######################################################################################################################
SimEstHR.Balkendiagramme <- function(file.erg1, file.erg2) {
    
    ### Vorbereitungen
    ## Definition der betrachteten (logarithmierten) Hazard-Ratios und eines Indizes für die Szenarien
    haz.rat <- c(1/3, 1/2, 2/3, 1/1.2, 1, 1.2, 1.5, 2, 3) 
    log.haz.rat <- log(haz.rat)
    idx.sz <- 1:50
    
    ## Definition der Namen und Farben der Punktschätzer
    name <- c(expression(S[C]), expression(S[B]), expression(S[W]), 
        expression(S[W]^{B[1]}), expression(S[W]^{B[2]}), expression(S[L]))
    name.B <- c(expression(S[C]^B), expression(S[B]^B), expression(S[W]^B), 
        expression(S[W]^list(B[1],B)), expression(S[W]^list(B[2],B)), expression(S[L]^B))
    name.exp <- c(expression(exp(S[C])), expression(exp(S[B])), expression(exp(S[W])), 
        expression(exp(S[W]^{B[1]})), expression(exp(S[W]^{B[2]})), expression(exp(S[L])))
    name.expB <- c(expression(exp(S[C]^B)), expression(exp(S[B]^B)), expression(exp(S[W]^B)), 
        expression(exp(S[W]^list(B[1],B))), expression(exp(S[W]^list(B[2],B))), expression(exp(S[L]^B)))

    mycols <- c(rgb(0.67, 0.8, 1), rgb(0.67, 1, 0.67), rgb(1, 0.67, 0.67), rgb(0.79, 0.67, 0.9), 
        rgb(1, 0.8, 0.5), rgb(1, 0.92, 0.67))
        
    ## Laden des Versuchsplans
    load(file = paste(file.erg1, "/SimEstHR_Info_Vplan.RData", sep = ""))
    
    ## Spezifikationen, die für den Titel der Abbildungen benötig werden
    Vplan.Zens <- ifelse(Vplan[, 6] == "keine", "nein", "ja")
    Vplan.Bdg <- ifelse(Vplan[, 7] == "keine", "nein", "ja")
    Vplan.Stoer <- ifelse(Vplan[, 8] == "keiner", "nein", "ja")
    
    ## Folgende Szenarien der alten Nummerierung bilden jeweils das erste Szenario in den zu zeichnenden Abbildungen
    sz.reihf <- c(3, 4, 1, 2, 7, 8, 5, 6, 11, 12, 9, 10, 15, 16, 13, 14, 49)    
    
    
    ### Zeichnung der Grafiken für je zwei bis drei der Szenarien in einer getrennten Datei ###########################
    ## Setzen eines Zählers, der die Nummern der jeweils ersten Szenarios in den Abbildungen durchläuft
    z <- -2
    
    ## Durchlauf der einzelnen Abbildungen
    for(i in sz.reihf) {
        
        ## Erstellung der Szenarioauswahl für eine Abbildung sowie
        ## Erzeugen der Überschriften der zwei bis drei Grafiken in einer Abbildung
        if(i < 49) {
            sz.ausw <- i + c(0, 16, 32)
            main.text <- c("Exp(variierend)/Exp(0.0693) (Gruppe 0/Gruppe 1)", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1)", 
                "Gomp(0.3, variierend)/Gomp(0.3, 0.0109) (Gruppe 0/Gruppe 1)")
        }else{ 
            sz.ausw <- c(49, 50)
            main.text <- c("Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 2", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 4")
        }
        
        z <- z + 3
        
        ## Öffnen eines pdfs und allgemeine Einstellungen für das Layout der Grafiken
        pdf(paste(file.erg2, "/SimEstHR_Balkendiagramme_", z, ".pdf", sep = ""), width = 21/2.54, height = 29.7/2.54)
        par(mfrow = c(3, 1), lwd = 1.1, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.45, 
            mar = c(4.5, 5, 3, 1), oma = c(0.1, 0.1, 3.5, 0.1))    
        
        ## Setzen eines weiteren Zählers für die betrachteten Szenarien innerhalb einer Abbildung
        zaehler <- 0    
        ## Durchlauf der zwei bzw. drei Szenarien einer Szenarioauswahl für eine Abbildung
        for(j in sz.ausw) {
            # Laden der Datei des entsprechenden Szenarios
            load(paste(file.erg1, "/SimEstHR_Sz", idx.sz[j], ".RData", sep = ""))
            
            zaehler <- zaehler + 1 
            
            # Extrahierung der Schätzwerte und Berechnung der MSE-Schätzungen
            MSE.mat <- matrix(0, nrow = 6, ncol = 9)   
            for(l in seq(along = haz.rat)) {
                est.temp <- sim.erg[[l]][, 1:6]
                MSE.mat[, l] <- apply(est.temp, 2, function(x) var(x, na.rm = TRUE)) + 
                    (apply(est.temp, 2, function(x) mean(x, na.rm = TRUE)) - log.haz.rat[l])^2        
            }
            
            # Zeichnen des Balkendiagramms für ein Szenario der aktuellen Szenarioauswahl
            max.bar <- max(apply(MSE.mat, 2, sum))
            for(l in seq(along = haz.rat)) {        
                MSE.mat.plot <- matrix(0, nrow = 6, ncol = 9)
                so.idx <- order(MSE.mat[, l])
                MSE.mat.plot[, l] <- MSE.mat[so.idx, l]
                            
                if(l == 1) {
                    xachse <- barplot(MSE.mat.plot, col = mycols[so.idx], 
                        ylim = c(0, c(1.2 * max.bar, 1.05 * max.bar, 1.05 * max.bar)[zaehler]),
                        xlab = expression(Wahres ~~ ln(omega)), ylab = "MSE", 
                        main = main.text[zaehler])
                    axis(1, at = xachse, labels = round(log.haz.rat, 2))
                    if(j <= 16 | j == 49) {
                        if(Vplan$Bdg[i] == "keine") {
                            legend("topleft", horiz = TRUE, legend = name, fill = mycols, border = mycols, 
                                cex = 1.5, pt.cex = 1.8)
                        } else {
                            legend("topleft", horiz = TRUE, legend = name.B, fill = mycols, border = mycols, 
                                cex = 1.5, pt.cex = 1.8)
                        }
                    }
                }else {barplot(MSE.mat.plot, col = mycols[so.idx], add = TRUE)}
            }
        }
        
        ## Erzeugen der Titel für die Abbildungen
        if(i < 49) {
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], 
                ", Störeffekt: ", Vplan.Stoer[i], sep = ""), outer = TRUE)
        }else{
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], sep = ""), outer = TRUE)
        }
        dev.off()
    }
    
    
    
    ### Zeichnung der Grafik für je zwei bis drei der Szenarien in einer gemeinsamen Datei ############################
    ## Öffnen eines pdfs
    pdf(paste(file.erg2, "/SimEstHR_Balkendiagramme_Alle.pdf", sep = ""), width = 21/2.54, height = 29.7/2.54)
        
    ## Durchlauf der einzelnen Abbildungen
    for(i in sz.reihf) {
        
        ## Allgemeine Einstellungen für das Layout der Grafiken
        par(mfrow = c(3, 1), lwd = 1.1, cex = 0.6, cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.45, 
        mar = c(4.5, 5, 3, 1), oma = c(0.1, 0.1, 3.5, 0.1))   
        
        ## Erstellung der Szenarioauswahl für eine Abbildung sowie
        ## Erzeugen der Überschriften der zwei bis drei Grafiken in einer Abbildung
        if(i < 49) {
            sz.ausw <- i + c(0, 16, 32)
            main.text <- c("Exp(variierend)/Exp(0.0693) (Gruppe 0/Gruppe 1)", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1)", 
                "Gomp(0.3, variierend)/Gomp(0.3, 0.0109) (Gruppe 0/Gruppe 1)")
        }else{ 
            sz.ausw <- c(49, 50)
            main.text <- c("Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 2", 
                "Weib(0.5, variierend)/Weib(0.5, 0.2192) (Gruppe 0/Gruppe 1), Störeffekt: 4")
        }
        
        ## Setzen eines Zählers für die betrachteten Szenarien innerhalb einer Abbildung
        zaehler <- 0    
        ## Durchlauf der zwei bzw. drei Szenarien einer Szenarioauswahl für eine Abbildung
        for(j in sz.ausw) {
            # Laden der Datei des entsprechenden Szenarios
            load(paste(file.erg1, "/SimEstHR_Sz", idx.sz[j], ".RData", sep = ""))
            
            zaehler <- zaehler + 1 
            
            # Extrahierung der Schätzwerte und Berechnung der MSE-Schätzungen
            MSE.mat <- matrix(0, nrow = 6, ncol = 9)   
            for(l in seq(along = haz.rat)) {
                est.temp <- sim.erg[[l]][, 1:6]
                MSE.mat[, l] <- apply(est.temp, 2, function(x) var(x, na.rm = TRUE)) + 
                    (apply(est.temp, 2, function(x) mean(x, na.rm = TRUE)) - log.haz.rat[l])^2        
            }
            
            # Zeichnen des Balkendiagramms für ein Szenario der aktuellen Szenarioauswahl
            max.bar <- max(apply(MSE.mat, 2, sum))
            for(l in seq(along = haz.rat)) {        
                MSE.mat.plot <- matrix(0, nrow = 6, ncol = 9)
                so.idx <- order(MSE.mat[, l])
                MSE.mat.plot[, l] <- MSE.mat[so.idx, l]
                            
                if(l == 1) {
                    xachse <- barplot(MSE.mat.plot, col = mycols[so.idx], 
                        ylim = c(0, c(1.2 * max.bar, 1.05 * max.bar, 1.05 * max.bar)[zaehler]),
                        xlab = expression(Wahres ~~ ln(omega)), ylab = "MSE", 
                        main = main.text[zaehler])
                    axis(1, at = xachse, labels = round(log.haz.rat, 2))
                    if(j <= 16 | j == 49) {
                        if(Vplan$Bdg[i] == "keine") {
                            legend("topleft", horiz = TRUE, legend = name, fill = mycols, border = mycols, 
                                cex = 1.5, pt.cex = 1.8)
                        } else {
                            legend("topleft", horiz = TRUE, legend = name.B, fill = mycols, border = mycols, 
                                cex = 1.5, pt.cex = 1.8)
                        }
                    }
                }else {barplot(MSE.mat.plot, col = mycols[so.idx], add = TRUE)}
            }
        }
        
        ## Erzeugen der Titel für die Abbildungen
        if(i < 49) {
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], 
                ", Störeffekt: ", Vplan.Stoer[i], sep = ""), outer = TRUE)
        }else{
            title(paste("Stichprobenumfang: ", Vplan[i, 4], "/", Vplan[i, 5], " (Gruppe 0/Gruppe 1),", 
                "\n Zensierungen: ", Vplan.Zens[i], ", Bindungen: ", Vplan.Bdg[i], sep = ""), outer = TRUE)
        }
    }
    dev.off()
}