#### For heatmap
# install.packages("BiocManager")
# BiocManager::install("ComplexHeatmap")

#library(ComplexHeatmap)
library(lattice)
source("params.R")
gcmStart <- 1
gcmLast  <- length(modelFileNames)
forecastNames <- modelFileNames

nmodels <- length(gcmStart:gcmLast)

dirfreq <- paste0(workingDir,"/epsfreq/")

domain_hss <- matrix(NA_real_, nrow = nmodels, ncol = ntest)
dimnames(domain_hss) <- list(modelNamesPlt, as.character(testYrsObs))

yy <- rep(1/3, 3)

for (imodel in 1:nmodels) {
  for (it in 1:ntest) {
    obs_cat <- tblw[, , it] + 2 * tnrm[, , it] + 3 * tabv[, , it]

    pr <- best_probs_by_model[, , it, , imodel]
    fc_cat <- apply(pr, c(1, 2), prob_to_cat, xpb = 0.33)

    ok <- is.finite(obs_cat) & is.finite(fc_cat)
    if (sum(ok) > 10) {
      xx <- table(
         factor(as.numeric( fc_cat[ok]), levels=1:3),
         factor(as.numeric(obs_cat[ok]), levels=1:3)
      )
      Acat <- multi.cont(xx, baseline=yy)
      domain_hss[imodel, it]  <- Acat$hss
    }
  }
}

hssfile <- paste0(workingDir,"/epsfreq/realTimeHSS_L", iLead, FTYPEF,".txt")
write.table(format(domain_hss), file=hssfile,quote=FALSE)

plt<- FALSE
if (plt) {
png(filename = paste0(workingDir,"/pheatrealTimeHSS_L",iLead, FTYPEF,".png"))
 pheatmap(domain_hss, cluster_rows = FALSE,
                 cluster_cols =FALSE,
                 cellwidth = 70, cellheight = 70,
                 color=hcl.colors(12, "YlOrRd", rev = TRUE),
                 breaks = c(0.0, 0.025, 0.05, 0.075, 0.10, 0.125,
                            0.150, 0.175, 0.20, 0.225, 0.25, 0.275,0.3))
 dev.off()
}

atLabels <- c(-Inf, 0.01, 0.025, 0.05, 0.075, 0.10, 0.125,
                            0.150, 0.175, 0.20, 0.225, 0.25, 0.275,0.3,Inf) 
numcols <- length(atLabels) + 3
png(file = paste0(dirfreq,"realTimeHSS_L",iLead, FTYPEF, ".png"))

 pp <-   levelplot(t(domain_hss),
            col.regions = hcl.colors(numcols, "YlOrRd",rev = TRUE),
            xlab = " Test Years ", ylab = " Models ",
	    main = paste0("Domain-wide independent hindcast HSS for ",SEASON),
            at = atLabels,
            scales = list(x = list(rot = 0),y = list(rot = 0))
        )
   print(pp)
  dev.off()

hssmax <- apply(best_hhss_by_model, c(1,2),which_max_na)
myprob <-  array(NA_real_, c(nlon, nlat,ntest, 3))


for (ilon in 1:nlon) {
   for (jlat in 1:nlat) { 
      itmp <- hssmax[ilon,jlat]
      if (is.finite(itmp) ) {
	 for (it in 1:ntest) {
           myprob[ilon, jlat,it,] <- best_probs_by_model[ilon,jlat,it,,itmp]
         }
      }
   }
}


for (it in 1:ntest) {
   tmp <- myprob[,,it,]
   fc_cat <- apply(tmp, c(1, 2), prob_to_cat, xpb = 0.33)
   obs_cat <- tblw[, , it] + 2 * tnrm[, , it] + 3 * tabv[, , it]
       ok <- is.finite(obs_cat) & is.finite(fc_cat)
   if (sum(ok) > 10) {
     xx <- table(
        factor(as.numeric( fc_cat[ok]), levels=1:3),
        factor(as.numeric(obs_cat[ok]), levels=1:3)
     )
     Acat <- multi.cont(xx, baseline=yy)
      print(Acat$hss)
     x <- as.matrix(Acat, ncol=3)
     freqfile <- paste0(dirfreq,"realTimeHSS_L", iLead, "_", testYrsObs[it], FTYPEF,".txt")
      write.csv(x, freqfile, row.names = TRUE)

    }
}



