R/plot_param.R

In hydroMOPSO: Multi-Objective Optimisation with Focus on Environmental Models

################################################################################
##                                plot_param                                  ##
################################################################################
# Author : Rodrigo Marinao Rivas                                              ##
################################################################################
# Created: 2024/08/31                                                         ##
################################################################################
# Function: plot_param
# Description: A function to generate parameter sensitivity plots, including 
#              boxplots and dotty plots, for the analysis of hydrological 
#              model parameters on the Pareto-optimal front. This function 
#              visualizes parameter distribution and variation to aid in 
#              understanding parameter influence on model performance.
################################################################################

plot_param <- function(
     Results,                # (list) Object containing preprocessed 
                             # hydrological results.
     legend.param = NULL,    # (character or NULL) Legend text for the 
                             # parameters in the plots.
     col = NULL,             # (character or NULL) Colors used for points in the
                             # parameter dotty plots.
     col.param = NULL,       # (character or NULL) Specific colors for lines or
                             # points representing parameters in the parameter 
                             # boxplots.
     col.lines = NULL,       # (character or NULL) Specific colors for lines in
                             # the parameter boxplots.
     name.param = NULL,      # (character or NULL) Custom names for the 
                             # parameters to be plotted.
     lwd = 2,                # (numeric) Line width for the parameter lines in
                             # the plots.
     main = "study case #1", # (character) Main title for the plot.
     drty.out = "MOPSO.out", # (character) Output directory where plots will be
                             # saved if specified to save as PNG files.
     cex.pt = 1,             # (numeric) Size of points in the plots.
     cex.main = 1,           # (numeric) Size of the main title text 
                             # in the plot.
     cex.lab = 1,            # (numeric) Size of the axis label text 
                             # in the plots.
     cex.axis = 1,           # (numeric) Size of the axis values text 
                             # in the plots.
     cex.leg = 1,            # (numeric) Size of the legend text in 
                             # the plots.
     do.png = FALSE          # (logical) Boolean value indicating whether the 
                             # plots should be saved as PNG files.
){


      BoxplotParametersOnPOF(Results,
                             legend.param = legend.param,
                             col.param = col.param,
                             col.lines = col.lines,
                             name.param = name.param,
                             lwd = lwd,
                             main = main,
                             drty.out = drty.out,
                             cex.main = cex.main,
                             cex.lab = cex.lab,
                             cex.axis = cex.axis,
                             cex.leg = cex.leg,
                             do.png = do.png)


      DottyplotParameters(Results,
                          legend.param = legend.param,
                          col = col,
                          name.param = name.param,
                          main = main,
                          drty.out = drty.out,
                          cex.pt = cex.pt,
                          cex.main = cex.main,
                          cex.lab = cex.lab,
                          cex.axis = cex.axis,
                          cex.leg = cex.leg,
                          do.png = do.png)


      DottyplotParametersPOF(Results,
                             legend.param = legend.param,
                             col = col,
                             name.param = name.param,
                             main = main,
                             drty.out = drty.out,
                             cex.pt = cex.pt,
                             cex.main = cex.main,
                             cex.lab = cex.lab,
                             cex.axis = cex.axis,
                             cex.leg = cex.leg,
                             do.png = do.png)


}




BoxplotParametersOnPOF <- function(Results,
                                   legend.param = NULL,
                                   col.param = NULL,
                                   col.lines = NULL,
                                   name.param = NULL,
                                   lwd = 2,
                                   main = "study case #1",
                                   drty.out = "MOPSO.out",
                                   cex.main = 1,
                                   cex.lab = 1,
                                   cex.axis = 1,
                                   cex.leg = 1,
                                   do.png = FALSE
                              ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(!is.null(Results[["hydroResults"]])){
    Results <- Results[["hydroResults"]]
  }


  analysis.period <- Results[["AnalysisPeriod"]]

  if(do.png){

    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }

 
  nobjs <- Results[["Dimensions"]][1,1]
  
  obj.names <- as.character(Results[["ObjsNames"]])

  parameter.set.pof <- Results[["ParticlesFilledPOF"]][,-c(1:(nobjs+1))]

  
  raw.some.particles <- matrix(NA, ncol = ncol(Results[["ParticleBestCS"]]),
                               nrow = 1 + length(Results[["ParticleBestObjs"]]))
  colnames(raw.some.particles) <- colnames(Results[["ParticleBestCS"]])
  raw.some.particles[1,] <- as.numeric(Results[["ParticleBestCS"]][1,])
  
  for(i in 1:length(Results[["ParticleBestObjs"]])){
    raw.some.particles[i+1,] <- as.numeric(Results[["ParticleBestObjs"]][[i]][1,])
  }

  some.particles <- raw.some.particles[,-c(1:(nobjs+1))]
  


  particle.bcs <- Results[["ParticleBestCS"]][-c(1:(1+nobjs))]
  objs.bcs <- Results[["ParticleBestCS"]][c(2:(1+nobjs))]


  particles.best.objs <- lapply(Results[["ParticleBestObjs"]], "[",
                                -c(1:(1+nobjs)))
  objs.best.objs <- lapply(Results[["ParticleBestObjs"]], "[", c(2:(1+nobjs)))




  nparam <- ncol(parameter.set.pof)



  samp.colors1 <- c("#A6CEE3","#1F78B4","#B2DF8A","#33A02C","#FB9A99","#E31A1C",
                   "#FDBF6F","#FF7F00","#CAB2D6","#6A3D9A","#FFFF99","#B15928",
                   "#8DD3C7","#FFFFB3","#BEBADA","#FB8072","#80B1D3","#FDB462",
                   "#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")

  samp.colors2 <- c("#E41A1C","#377EB8","#4DAF4A","#984EA3","#FF7F00",
                    "#FFFF33","#A65628","#F781BF","#999999")

  
  # Parameters on Pareto Front #################################################
  
  if(do.png){
    png(filename=paste0(drty.out, "/", analysis.period,
                    "/png.out/Parameters_on_Pareto_Optimal_Front_Boxplots.png"),
        width = 3840, height = 2160, res = 280)
  }
  
  nplots <- nparam + 1
  
  nrow.lay <- floor(sqrt(nplots))
  ncol.lay <- floor(sqrt(nplots))+1
  
  ncol.lay <- ifelse(nrow.lay*ncol.lay<nplots, ncol.lay + 1, ncol.lay)
  
  par(mfrow = c(nrow.lay, ncol.lay), oma = c(0, 0, 2, 0), mar = c(2,4,3,2))
  
  ###



  colors.sol <- rep("", length(nobjs))

  if(nparam<=24){
    colors.sol1 <- samp.colors1[1:nparam]
  }else{
    colors.sol1 <- c(samp.colors1, sample(samp.colors1, size = nparam - 24,
                                          replace = TRUE))
  }
  

  if(!is.null(col.param)){

    min.col <- min(length(col.param),length(colors.sol1))

    colors.sol1[1:min.col] <- col.param[1:min.col]

  }


  #--

  if(nrow(some.particles)<=9){
    colors.sol2 <- samp.colors2[1:nrow(some.particles)]
  }else{
    colors.sol2 <- c(samp.colors2, sample(samp.colors2,
                                          size = nrow(some.particles) - 9,
                                          replace = TRUE))
  }


  if(!is.null(col.lines)){

    min.col <- min(length(col.lines),length(colors.sol2))

    colors.sol2[1:min.col] <- col.lines[1:min.col]

  }

  #----


  param.names <- colnames(parameter.set.pof)



  if(!is.null(name.param)){

    min.col <- min(length(name.param),length(param.names))

    param.names[1:min.col] <- name.param[1:min.col]

  }


  if(file.exists(paste0("MOPSO.in/ParamRanges.txt"))){

    file.ranges <- read.table(paste0("MOPSO.in/ParamRanges.txt"), header = TRUE,
                              row.names = 2)

    type.change <- file.ranges[param.names, "TypeChange"]

    param.names <- apply(cbind(param.names, paste0(type.change, " change")), 1,
                         paste, collapse = "\n")

  }



  #---------


  
  ltype <- sample(c("dashed", "dotted", "dotdash", "longdash", "twodash"),
                  size = nrow(some.particles), replace = TRUE)
  
  for(i in 1:nparam){
    boxplot(parameter.set.pof[,i], main = param.names[i],
            col = colors.sol1[i])
    
    for(j in 1:nrow(some.particles)){
      abline(h = some.particles[j,i], col = colors.sol2[j], lty = ltype[j],
             lwd = lwd)
    }
    
    
  }
  

  names.sol <- c("Best Compromise Solution", paste0("Best ",obj.names))

  if(!is.null(legend.param)){

    min.leg1 <- min(length(legend.param),length(names.sol))
    min.leg2 <- min(length(legend),length(names.sol))

    names.sol[1:min.leg1] <- legend[1:min.leg2]

  }

  plot(c(0,1),c(0,1),type = "n", axes = FALSE, xlab = "", ylab = "",
      cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis)

  legend(x = 0.5, y = 0.5,xjust = 0.5, yjust = 0.5,
         legend = names.sol, 
         col = colors.sol2, 
         lty = ltype,
         lwd = lwd,
         bty = "n", 
         cex = cex.leg, 
         text.col = "black", 
         horiz = FALSE, 
         inset = c(0.1, 0.1))
  
  mtext("Parameters on Pareto Front", outer = TRUE, cex = 1.5, font = 2)
  
  if(do.png){
    dev.off()
  }
  
}







DottyplotParameters <- function(Results,
                                   legend.param = legend.param,
                                   col = NULL,
                                   name.param = NULL,
                                   main = "study case #1",
                                   drty.out = "MOPSO.out",
                                   cex.pt = 1,
                                   cex.main = 1,
                                   cex.lab = 1,
                                   cex.axis = 1,
                                   cex.leg = 1,
                                   do.png = FALSE
                              ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(!is.null(Results[["hydroResults"]])){
    Results <- Results[["hydroResults"]]
  }


  analysis.period <- Results[["AnalysisPeriod"]]

  if(do.png){

    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }

 
  nobjs <- Results[["Dimensions"]][1,1]
  ObjThr <- as.matrix(Results[["ObjsThresholds"]])
  
  obj.names <- as.character(Results[["ObjsNames"]])

  parameter.set.pof <- Results[["ParticlesFilledPOF"]][,-c(1:(nobjs+1))]
  obj.set.pof <- Results[["ParticlesFilledPOF"]][,c(2:(nobjs+1))]

  parameter.set.sub <- Results[["ParticlesFull"]][,-c(1:(nobjs+1))]
  obj.set.sub <- Results[["ParticlesFull"]][,c(2:(nobjs+1))]
  
  objs <- Results[["FilledPOF"]][,-1]

  
  raw.some.particles <- matrix(NA, ncol = ncol(Results[["ParticleBestCS"]]),
                               nrow = 1 + length(Results[["ParticleBestObjs"]]))
  colnames(raw.some.particles) <- colnames(Results[["ParticleBestCS"]])
  raw.some.particles[1,] <- as.numeric(Results[["ParticleBestCS"]][1,])
  
  for(i in 1:length(Results[["ParticleBestObjs"]])){

  raw.some.particles[i+1,] <- as.numeric(Results[["ParticleBestObjs"]][[i]][1,])

  }

  some.particles <- raw.some.particles[,-c(1:(nobjs+1))]
  


  particle.bcs <- Results[["ParticleBestCS"]][-c(1:(1+nobjs))]
  objs.bcs <- Results[["ParticleBestCS"]][c(2:(1+nobjs))]


  particles.best.objs <- lapply(Results[["ParticleBestObjs"]], "[",
                                -c(1:(1+nobjs)))
  objs.best.objs <- lapply(Results[["ParticleBestObjs"]], "[", c(2:(1+nobjs)))


  ylim_obj_min <- apply(obj.set.pof,2,min)
  ylim_obj_max <- apply(obj.set.pof,2,max)

  ylimObj_min  <- pmax(ObjThr[1,], ylim_obj_min)
  ylimObj_max  <- pmin(ObjThr[2,], ylim_obj_max)

  ylimObj <- rbind(ylimObj_min, ylimObj_max)


  nparam <- ncol(parameter.set.pof)



  samp.colors3 <- c("#a000b8","#0aab05","#d16800", "#FFFF33","#A65628",
                    "#F781BF","#999999")
  
  
  # Parameter Dottyplots #######################################################
  
  if(analysis.period == "calibration"){
  
    for(j in 1:nobjs){

      if(do.png){
        png(filename=paste0(drty.out, "/", analysis.period,
                            "/png.out/Parameter_Dottyplots_Obj",j,".png"),
            width = 3840, height = 2160, res = 280)
      }else{
        dev.new()
      }

      nplots <- nparam + 1
      
      nrow.lay <- floor(sqrt(nplots))
      ncol.lay <- floor(sqrt(nplots))+1
      
      ncol.lay <- ifelse(nrow.lay*ncol.lay<nplots, ncol.lay + 1, ncol.lay)
      
      par(mfrow = c(nrow.lay, ncol.lay), oma = c(0, 0, 2, 0), mar = c(3,4,3,2))
      
      ###
      

      colors <- c("orangered", "gray20", "#004fcf", samp.colors3)
      

      if(!is.null(col)){

        min.col <- min(length(colors),length(col))

        colors[1:min.col] <- col[1:min.col]

      }

      #-------


      param.names <- colnames(parameter.set.sub)

min(length(colors),length(col))
      if(!is.null(name.param)){

        min.col <- min(length(name.param),length(param.names))

        param.names[1:(min.col)] <- name.param[1:(min.col)]

      }


      if(file.exists(paste0("MOPSO.in/ParamRanges.txt"))){

        file.ranges <- read.table(paste0("MOPSO.in/ParamRanges.txt"),
                                  header = TRUE, row.names = 2)

        type.change <- file.ranges[param.names, "TypeChange"]

        param.names <- apply(cbind(param.names, paste0(type.change, " change")),
                             1, paste, collapse = "\n")

      }

      


      for(i in 1:nparam){
        plot(x = parameter.set.pof[,i], y = obj.set.pof[,j],
             cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis,
             xlab = colnames(parameter.set.sub)[i], ylab = obj.names[j], 
             main = param.names[i], type = "n")
        points(x = parameter.set.sub[,i], y = obj.set.sub[,j], col = colors[2],
               cex = cex.pt)
        points(x = parameter.set.pof[,i], y = obj.set.pof[,j], col = colors[1],
               cex = cex.pt)
        points(x = particle.bcs[i], y = objs.bcs[j], bg = colors[3], pch = 22,
               cex = 1.75*cex.pt)
        for(h in 1:length(objs)){
          points(x = particles.best.objs[[h]][i], y = objs.best.objs[[h]][j],
                 bg = colors[h+3], pch = 24, cex = 1.75*cex.pt)
        }

      }

      colors_black <- colors
      colors_black[-c(1:2)] <- "black"


      names.sol <- c("Pareto-optimal solutions","Sub-optimal solutions", 
                     "Best compromise solution",
                     paste0("Best optimal for single ", obj.names))

      if(!is.null(legend.param)){

        min.leg1 <- min(length(legend.param),length(names.sol))
        min.leg2 <- min(length(legend),length(names.sol))

        names.sol[1:min.leg1] <- legend.param[1:min.leg2]

      }
  

      plot(c(0,1),c(0,1),type = "n", axes = FALSE, xlab = "", ylab = "")
      legend(x = 0.5, y = 0.5,xjust = 0.5, yjust = 0.5,
             legend = names.sol, 
             col = colors_black,
             pt.bg = colors,
             lty = 0,
             pch = c(1,1, 22, rep(24, nobjs)),
             lwd = c(1,1, NA, rep(NA, nobjs)),
             bty = "n", 
             pt.cex = c(cex.pt, cex.pt, 1.5*cex.pt, rep(1.75*cex.pt, nobjs)), 
             cex = cex.leg, 
             text.col = "black", 
             horiz = FALSE, 
             inset = c(0.1, 0.1),
             y.intersp = 1.4)
        
      
      mtext(paste0("Parameters (Obj ",j,": ", obj.names[j],")"), outer = TRUE,
            cex = 1.5, font = 2)
      
      if(do.png){
        dev.off()
      }
    }
  
  }

}




DottyplotParametersPOF <- function(Results,
                                   legend.param = legend.param,
                                   col = NULL,
                                   name.param = NULL,
                                   main = "study case #1",
                                   drty.out = "MOPSO.out",
                                   cex.pt = 1,
                                   cex.main = 1,
                                   cex.lab = 1,
                                   cex.axis = 1,
                                   cex.leg = 1,
                                   do.png = FALSE
                              ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(!is.null(Results[["hydroResults"]])){
    Results <- Results[["hydroResults"]]
  }


  analysis.period <- Results[["AnalysisPeriod"]]

  if(do.png){

    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }

 
  nobjs <- Results[["Dimensions"]][1,1]
  ObjThr <- as.matrix(Results[["ObjsThresholds"]])
  
  obj.names <- as.character(Results[["ObjsNames"]])

  parameter.set.pof <- Results[["ParticlesFilledPOF"]][,-c(1:(nobjs+1))]
  obj.set.pof <- Results[["ParticlesFilledPOF"]][,c(2:(nobjs+1))]

  parameter.set.sub <- Results[["ParticlesFull"]][,-c(1:(nobjs+1))]
  obj.set.sub <- Results[["ParticlesFull"]][,c(2:(nobjs+1))]
  
  objs <- Results[["FilledPOF"]][,-1]

  
  raw.some.particles <- matrix(NA, ncol = ncol(Results[["ParticleBestCS"]]),
                               nrow = 1 + length(Results[["ParticleBestObjs"]]))
  colnames(raw.some.particles) <- colnames(Results[["ParticleBestCS"]])
  raw.some.particles[1,] <- as.numeric(Results[["ParticleBestCS"]][1,])
  
  for(i in 1:length(Results[["ParticleBestObjs"]])){
  raw.some.particles[i+1,] <- as.numeric(Results[["ParticleBestObjs"]][[i]][1,])
  }

  some.particles <- raw.some.particles[,-c(1:(nobjs+1))]
  


  particle.bcs <- Results[["ParticleBestCS"]][-c(1:(1+nobjs))]
  objs.bcs <- Results[["ParticleBestCS"]][c(2:(1+nobjs))]


  particles.best.objs <- lapply(Results[["ParticleBestObjs"]], "[",
                                -c(1:(1+nobjs)))
  objs.best.objs <- lapply(Results[["ParticleBestObjs"]], "[",
                           c(2:(1+nobjs)))


  ylim_obj_min <- apply(obj.set.pof,2,min)
  ylim_obj_max <- apply(obj.set.pof,2,max)

  ylimObj_min  <- pmax(ObjThr[1,], ylim_obj_min)
  ylimObj_max  <- pmin(ObjThr[2,], ylim_obj_max)

  ylimObj <- rbind(ylimObj_min, ylimObj_max)


  nparam <- ncol(parameter.set.pof)



  samp.colors3 <- c("#a000b8","#0aab05","#d16800", "#FFFF33","#A65628",
                    "#F781BF","#999999")
  
  
  # Parameter Dottyplots #######################################################
  
  if(analysis.period == "calibration"){
  
    for(j in 1:nobjs){

      if(do.png){
        png(filename=paste0(drty.out, "/", analysis.period,
                           "/png.out/Parameter_Dottyplots_in_POF_Obj",j,".png"),
                            width = 3840, height = 2160, res = 280)
      }else{
        dev.new()
      }

      nplots <- nparam + 1
      
      nrow.lay <- floor(sqrt(nplots))
      ncol.lay <- floor(sqrt(nplots))+1
      
      ncol.lay <- ifelse(nrow.lay*ncol.lay<nplots, ncol.lay + 1, ncol.lay)
      
      par(mfrow = c(nrow.lay, ncol.lay), oma = c(0, 0, 2, 0), mar = c(3,4,3,2))
      
      ###
      

      colors <- c("orangered", "#004fcf", samp.colors3)
      

      if(!is.null(col)){

        min.col <- min(length(colors),length(col))

        colors[1:(min.col)] <- col[1:(min.col)]

      }

      #-------


      param.names <- colnames(parameter.set.sub)


      if(!is.null(name.param)){

        min.col <- min(length(name.param),length(param.names))

        param.names[1:(min.col)] <- name.param[1:(min.col)]

      }



      if(file.exists(paste0("MOPSO.in/ParamRanges.txt"))){

        file.ranges <- read.table(paste0("MOPSO.in/ParamRanges.txt"),
                                  header = TRUE, row.names = 2)

        type.change <- file.ranges[param.names, "TypeChange"]

        param.names <- apply(cbind(param.names, paste0(type.change, " change")),
                             1, paste, collapse = "\n")

      }





      for(i in 1:nparam){
        plot(x = parameter.set.pof[,i], y = obj.set.pof[,j],
             ylim = if(any(is.na(ylimObj[,j]))) NULL else ylimObj[,j],
             cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis,
             xlab = colnames(parameter.set.sub)[i], ylab = obj.names[j], 
             main = param.names[i], cex = cex.pt, type = "n")
        points(x = parameter.set.pof[,i], y = obj.set.pof[,j], col = colors[1],
               cex = cex.pt)
        points(x = particle.bcs[i], y = objs.bcs[j], bg = colors[2], pch = 22,
               cex = 1.75*cex.pt)
        for(h in 1:length(objs)){
          points(x = particles.best.objs[[h]][i], y = objs.best.objs[[h]][j],
                 bg = colors[h+2], pch = 24, cex = 1.75*cex.pt)
        }

      }

      colors_black <- colors
      colors_black[-c(1:2)] <- "black"



      names.sol <- c("Pareto-optimal solutions", "Best compromise solution",
                     paste0("Best optimal for single ", obj.names))

      if(!is.null(legend.param)){

        min.leg1 <- min(length(legend.param),length(names.sol))
        min.leg2 <- min(length(legend),length(names.sol))

        names.sol[1:min.leg1] <- legend.param[1:min.leg2]

      }

      plot(c(0,1),c(0,1),type = "n", axes = FALSE, xlab = "", ylab = "")
      legend(x = 0.5, y = 0.5,xjust = 0.5, yjust = 0.5,
             legend = names.sol, 
             col = colors_black,
             pt.bg = colors,
             lty = 0,
             pch = c(1, 22, rep(24, nobjs)),
             lwd = c(1, NA, rep(NA, nobjs)),
             bty = "n", 
             pt.cex = c(cex.pt, 1.5*cex.pt, rep(1.75*cex.pt, nobjs)), 
             cex = cex.leg, 
             text.col = "black", 
             horiz = FALSE, 
             inset = c(0.1, 0.1),
             y.intersp = 1.4)
        
      
      mtext(paste0("Parameters on POF (Obj ",j,": ", obj.names[j],")"),
            outer = TRUE, cex = 1.5, font = 2)
      
      if(do.png){
        dev.off()
      }
    }
  
  }

}