rm(list=ls())
library(mgcv)


# Load simulation estimates and aggregate resulting MwSE and PMWSE
n_text <- c("5000")

list_range <- list()
range_m <- matrix(NA,nrow=11,ncol=5)
colnames(range_m) <- c("min","max","diff","ylim_min","ylim_max")
asy_vec <- c("01","02","05","10","20","50","80","90","95","98","99")
rownames(range_m) <- asy_vec

list_range[[2]] <- list_range[[1]] <- range_m
names(list_range) <- n_text


for(ntext_in in n_text) {
    MWSE_matrix_new <- NULL
    PMWSE_matrix_new <- NULL
    for(Namein in c("heterosc_1","heterosc_2","homosc_1","homosc_2")){
        
        # Build name of files
        Namein_1  <- paste(strsplit(Namein,split="_",fixed=T)[[1]][1],"_T_",strsplit(Namein,split="_",fixed=T)[[1]][2],sep="")
        Name_load <- paste("Sim4_",ntext_in,"_100_ti_T_",Namein_1,".RData",sep="")
        
        # Load file
        load(Name_load)
        
        MWSE_temp  <-  NULL
        PMWSE_temp <- NULL
        
        # Choose columns of result matrices
        MWSE_temp  <-  MWSE_matrix[,which(grepl(x=colnames( MWSE_matrix),pattern="bam_gcv") | grepl(x=colnames( MWSE_matrix),pattern="bam_sch_d"))]
        PMWSE_temp <- PMWSE_matrix[,which(grepl(x=colnames(PMWSE_matrix),pattern="bam_gcv") | grepl(x=colnames(PMWSE_matrix),pattern="bam_sch_d"))]
        
        # Rename columns
        colnames( MWSE_temp) <- colnames( MWSE_matrix)[which(grepl(x=colnames( MWSE_matrix),pattern="bam_gcv") | grepl(x=colnames( MWSE_matrix),pattern="bam_sch_d"))]
        colnames(PMWSE_temp) <- colnames(PMWSE_matrix)[which(grepl(x=colnames(PMWSE_matrix),pattern="bam_gcv") | grepl(x=colnames(PMWSE_matrix),pattern="bam_sch_d"))]
        colnames( MWSE_temp) <- paste(Namein,colnames( MWSE_temp),sep="_")
        colnames(PMWSE_temp) <- paste(Namein,colnames(PMWSE_temp),sep="_")
        MWSE_matrix_new <- cbind( MWSE_matrix_new, MWSE_temp)
        PMWSE_matrix_new <- cbind(PMWSE_matrix_new,PMWSE_temp)

    }
    # Assign matrix based on number of observations
    if(ntext_in == "2000") {
        MWSE_matrix_2000 <-  MWSE_matrix_new
        PMWSE_matrix_2000 <- PMWSE_matrix_new
    }
    if(ntext_in == "5000") {
        MWSE_matrix_5000 <-  MWSE_matrix_new
        PMWSE_matrix_5000 <- PMWSE_matrix_new
    }
}

# To save space and to avoid cluttering MWSE results are not published


for(ntext_in in n_text) {
    range_m <- list_range[[ntext_in]]
    for(MWS in c("PMWSE")) {
        for(asy_in in asy_vec){
            
            # Get the current result matrix
            MWSE_matrix_new <- get(paste( MWS,"_matrix_",ntext_in,sep=""))
            
            
            # Shrink colnames to better fit in plots
            colnames( MWSE_matrix_new) <- gsub(x=colnames( MWSE_matrix_new),pattern="_bam"    ,replacement="",fixed=T)
            colnames( MWSE_matrix_new) <- gsub(x=colnames( MWSE_matrix_new),pattern="Sim_100_",replacement="",fixed=T)
            
            # Extract estimates of different asymmetries to rearange for plots
            MWSE_matrix_new_a <-  MWSE_matrix_new[,grepl(x=colnames( MWSE_matrix_new),pattern=paste("_",asy_in,"_",sep=""),fixed=T)]
            
            range_m[asy_in,c("min","max")] <- range(MWSE_matrix_new_a)
             
            
        }
    }
    range_m[,"diff"]     <- range_m[,"max"] - range_m[,"min"]
    range_m[,"ylim_min"] <- floor(range_m[,"min"]*100)/100
    if(ntext_in == "5000") {
        range_m[,"ylim_max"] <- floor(range_m[,"min"]*100)/100 + 0.1
    } else {
        range_m[,"ylim_max"] <- floor(range_m[,"min"]*100)/100 + 0.3
    }   
    
    list_range[[ntext_in]] <- range_m 
    
    
}



# Plot results as boxplots, for each number of observations (ntext) each 
# goodness-of-fit criterion (MWS) and each asymmetry level (asy) a seperate boxplot

pdf("Aggregate_Simulation_Paper.pdf",width=4,height=4.5)

for(ntext_in in n_text) {
    for(MWS in c("PMWSE")) {
        for(asy_in in asy_vec){
            
            # Get the current result matrix
            MWSE_matrix_new <- get(paste( MWS,"_matrix_",ntext_in,sep=""))
            
            
            # Shrink colnames to better fit in plots
            colnames( MWSE_matrix_new) <- gsub(x=colnames( MWSE_matrix_new),pattern="_bam"    ,replacement="",fixed=T)
            colnames( MWSE_matrix_new) <- gsub(x=colnames( MWSE_matrix_new),pattern="Sim_100_",replacement="",fixed=T)
            
            # Extract estimates of different asymmetries to rearange for plots
            MWSE_matrix_new_a <-  MWSE_matrix_new[,grepl(x=colnames( MWSE_matrix_new),pattern=paste("_",asy_in,"_",sep=""),fixed=T)]
            
            # Shrink colnames to better fit in plots
            colnames( MWSE_matrix_new_a) <- gsub(x=colnames( MWSE_matrix_new_a),pattern=paste("_",asy_in,sep=""),replacement="",fixed=T)
            
            colnames( MWSE_matrix_new_a) <- gsub(x=colnames( MWSE_matrix_new_a),pattern="sc_",replacement=" ",fixed=T)
            
            colnames( MWSE_matrix_new_a) <- gsub(x=colnames( MWSE_matrix_new_a),pattern="_",replacement=" ",fixed=T)
            
            colnames( MWSE_matrix_new_a) <- gsub(x=colnames( MWSE_matrix_new_a),pattern="sch",replacement="schall",fixed=T)
            
            # Plot results in one picture
            # PMWSE has range of 0.1 
            ylim2 <- NULL
            
            ylim2 <- list_range[[ntext_in]][asy_in, c("ylim_min","ylim_max")]
            par(mar=c(6,4,2,1))    
            ylab1 <- MWS
            boxplot(MWSE_matrix_new_a,las=2,main=paste(asy_in,"% Expectile", sep=" "),names=rep(c("gcv","schall"),times=4),ylim=ylim2,ylab=ylab1)#
            mtext(at=c(1.5),text="1",side=1,line=3.5)
            mtext(at=c(3.5),text="2",side=1,line=3.5)
            mtext(at=c(5.5),text="1",side=1,line=3.5)
            mtext(at=c(7.5),text="2",side=1,line=3.5)
            mtext(at=c(2.5),text="heterosc",side=1,line=4.5)
            mtext(at=c(6.5),text="homosc",side=1,line=4.5)
            abline(v=2.5,lty=3)
            abline(v=4.5,lty=2)
            abline(v=6.5,lty=3)
            
            
        }
    }
}
dev.off()


