R/plot_pof.R

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

################################################################################
##                                plot_pof                                    ##
################################################################################
# Author : Rodrigo Marinao Rivas                                              ##
################################################################################
# Created: 2024/08/31                                                         ##
################################################################################
# Function: plot_pof
# Description: A function to generate plots of the Pareto-optimal front (POF) 
#              and the best compromise solution (BCS) in multi-objective 
#              optimization for hydrological models. This function visualizes 
#              the trade-offs between different objectives and identifies the 
#              most balanced solution, aiding in model calibration and 
#              evaluation processes.
################################################################################

plot_pof <- function(
     Results,                # (list) Object containing preprocessed 
                             # hydrological results.
     pof = NULL,             # (matrix or NULL) Dataset representing the filled 
                             # Pareto-optimal front (POF) solutions.
     bcs = NULL,             # (matrix or NULL) Parameters or results 
                             # corresponding to the best compromise solution 
                             # (BCS).
     analysis.period = NULL, # (character or NULL) Time period for analysis 
                             # ("calibration" or "verification").
     dimensions = NULL,      # (numeric vector or NULL) Dimensions of the 
                             # modeled problem, typically a vector indicating 
                             # the number of objectives and variables.
     maxmin = NULL,          # (character or NULL) Indicator of whether 
                             # objectives are to be maximized ("max") or 
                             # minimized ("min").
     obj.thr = NULL,         # (numeric vector or NULL) Objective thresholds.
     obj.names = NULL,       # (character vector or NULL) Names of the 
                             # objectives used in the plot.
     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.
     pch.pof = 21,           # (integer) Symbol type for points in the plot 
                             # representing the Pareto-optimal front solutions.
     pch.bcs = 21,           # (integer) Symbol type for points in the plot 
                             # representing the best compromise solution (BCS).
     col.pof = "#f21b1b",    # (character) Color used for the points 
                             # representing the Pareto-optimal front solutions 
                             # in the plots.
     col.bcs = "#004fcf",    # (character) Color used for the points 
                             # representing the best compromise solution in the
                             # plots.
     legend.pof = c("Pareto-optimal front solutions",
                    "Best compromise solution"),
                             # (character vector) Legend text for the 
                             # Pareto-optimal front and best compromise 
                             # solution.
     cex.pt = 1.25,          # (numeric) Size of points in the plots.
     cex.main = 1,           # (numeric) Size of the main title text in plot.
     cex.lab = 1,            # (numeric) Size of the axis label text in plots.
     cex.axis = 1,           # (numeric) Size of the axis values text in plots.
     do.png = FALSE          # (logical) Boolean value indicating whether the 
                             # plots should be saved as PNG files.
){

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

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

  results.names <- c("ParticlesFilledPOF", "ParticleBestCS", "AnalysisPeriod",
                     "Dimensions", "MaxMin", "ObjsThresholds")

  if(!is.null(Results)){
      if(all(results.names %in% names(Results))){
          flag.results <- TRUE
      }else{
          flag.results <- FALSE
      }
  }else{
      flag.results <- FALSE
  }

  #


  if(is.null(pof)){
    if(flag.results){
      pof <- Results[["ParticlesFilledPOF"]]        
    }else{
      stop("The 'pof' argument was not provided or included in the 'Results'
            argument. Please make sure to include all required arguments, and
            if using the 'Results' argument, ensure 'ParticlesFilledPOF' is
            correctly specified.")
    }
  }

  if(is.null(bcs)){
    if(flag.results){
      bcs <- Results[["ParticleBestCS"]]        
    }else{
      stop("The 'bcs' argument was not provided or included in the 'Results'
            argument. Please make sure to include all required arguments, and if
            using the 'Results' argument, ensure 'ParticleBestCS' is correctly
            specified.")
    }
  }

  if(is.null(analysis.period)){
    if(flag.results){
      analysis.period <- Results[["AnalysisPeriod"]]        
    }else{
      stop("The 'analysis.period' argument was not provided or included in the
            'Results' argument. Please make sure to include all required
            arguments, and if using the 'Results' argument, ensure
            'AnalysisPeriod' is correctly specified.")
    }
  }

  if(is.null(dimensions)){
    if(flag.results){
      dimensions <- Results[["Dimensions"]]        
    }else{
      stop("The 'dimensions' argument was not provided or included in the
            'Results' argument. Please make sure to include all required
            arguments, and if using the 'Results' argument, ensure 'Dimensions'
            is correctly specified.")
    }
  }

  if(is.null(maxmin)){
    if(flag.results){
      maxmin <- Results[["MaxMin"]]        
    }else{
      stop("The 'maxmin' argument was not provided or included in the 'Results'
            argument. Please make sure to include all required arguments, and if
            using the 'Results' argument, ensure 'MaxMin' is correctly
            specified.")
    }
  }

  if(is.null(obj.thr)){
    if(flag.results){
      obj.thr <- Results[["ObjsThresholds"]]        
    }else{
      stop("The 'obj.thr' argument was not provided or included in the 'Results'
            argument. Please make sure to include all required arguments, and if
            using the 'Results' argument, ensure 'ObjsThresholds' is correctly
            specified.")
    }
  }




  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 <- dimensions[1,1]
  nvars <- dimensions[1,2]


  objs.bcs <- bcs[1,c(2:(1+nobjs))]


  if(is.null(obj.names)){
    obj.names <- colnames(objs.bcs)
  }
  

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

  objs.bcs <- bcs[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(obj.thr[1,], ylim_obj_min)
  ylimObj_max  <- pmin(obj.thr[2,], ylim_obj_max)

  ylimObj <- rbind(ylimObj_min, ylimObj_max)

  # Pareto Fronts ##############################################################

  if(analysis.period == "calibration"){
    main_2d <- "Pareto Optimal Front in calibration period"
  }else if(analysis.period == "verification"){
    main_2d <- "Degradation of the Pareto Optimal Front in verification period"
  }

  if(do.png){
    png(filename=paste0(drty.out, "/", analysis.period,
                        "/png.out/Solutions_2D_Compromise.png"),
        width = 3000, height = 3000, res = 500)
  }else{
    dev.new()
  }

  egrid <- expand.grid(1:nobjs,1:nobjs)

  flag_1 <- !(apply(egrid,1,diff) == 0)
  flag_2 <- !duplicated(apply(t(apply(egrid,1,sort)),1,paste,collapse = "_"))

  comb <- egrid[flag_1 & flag_2,]

  nplots <- nrow(comb)

  nrow.lay <- floor(sqrt(nplots))
  ncol.lay <- floor(sqrt(nplots))

  ncol.lay <- ifelse(nrow.lay*ncol.lay<nplots, ncol.lay, ncol.lay)

  par(mfrow = c(nrow.lay, ncol.lay), oma = c(0, 0, 2, 0))

  for(i in 1:nplots){

    aa <- comb[i,1]
    bb <- comb[i,2]

    plot(obj.set.pof[,aa], obj.set.pof[,bb], xlab = obj.names[aa],
         ylab = obj.names[bb],
         main = paste0(main, "\n", main_2d), cex.main = cex.main,
         cex.lab = cex.lab, cex.axis = cex.axis,
         xlim = if(any(is.na(ylimObj[,aa]))) NULL else ylimObj[,aa],
         ylim = if(any(is.na(ylimObj[,bb]))) NULL else ylimObj[,bb])
    grid()
    points(obj.set.pof[,aa], obj.set.pof[,bb], xlab = obj.names[aa],
           ylab = obj.names[bb], pch = pch.pof, bg = col.pof, col = "black",
           cex = cex.pt)
    points(x = objs.bcs[aa], y = objs.bcs[bb], bg = col.bcs, col = "black",
           pch = pch.bcs, cex = cex.pt)

    legend("bottomleft",
           legend = legend.pof, 
           pt.bg = c(col.pof, col.bcs), 
           lty = 0,
           pch = c(pch.pof,pch.bcs),
           lwd = 1, 
           pt.cex = cex.pt,
           cex = 0.75, 
           text.col = "black", 
           horiz = FALSE, 
           inset = c(0.1, 0.1))
  }


  if(do.png){
    dev.off()
  }


}