R/DottyplotParametersPOF.R

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

DottyplotParametersPOF <- function(Results,
                                   legend = NULL,
                                   col = NULL,
                                   name.param = NULL,
                                   lwd = 1,
                                   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)){

        colors[1:(min(length(colors),length(col)))] <- col[1:(min(length(colors),length(col)))]

      }

      #-------


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


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

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

      }



      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"

      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 = c("Pareto-optimal solutions", "Best compromise solution", paste0("Best optimal for single ", obj.names)), 
             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()
      }
    }
  
  }

}