seedHealth: Scenario Analysis for Seed Health

#' Plots of output for multiple parameter combinations
#'
#' This function plots a heatmap and a 3-D graph of output for multiple parameter combinations across nseasons
#'
#' The truncated random normal variables were generated by \code{\link{altruncnorm}}. 
#' Refer to \code{\link{onesim}}, \code{\link{multisim}}, and \code{\link{multipar}} for relevant information.
#'
#' @inheritParams onesim
#' @param nsim vector of number of simulations
#' @import tibble
#' @import ggplot2
#' @import lattice
#' @importFrom stats median quantile rnorm var
#' @keywords seed health
#' @examples
#' multipar_plot() # more examples to be added
#' @export
#'

# to do - GT

# Columns of output matrix
# col 1 - fHP healthy plant number
# col 2 - fDP diseased plant number (after roguing)
# col 3 - fHS healthy seed number
# col 4 - fDS diseased seed number
# col 5 - fpHS proportion healthy seed
# col 6 - fpDS proportion diseased seed
# col 7 - HPtrans Season in which HP first transitions below HPcut*Kx, if HPtrans is NA i.e., HP never less than HPcut, set to max seasons tested
# col 8 - pHStrans Season in which pHS first transitions below pHScut, if pHStrans is NA i.e., pHS never less than pHScut, set to max seasons tested
# col 9 - HPpseas Healthy plants are calculated from season 1 onwards
# col 10 - pHSpseas Proportion seasons with pHS below pHScut
# col 11 - fYld end of season yield
# col 12 - fYL end of season yield loss

# mean (mean), median (median),  variance (var), quantile 0.05 (0.05), quantile 0.95 (0.95)

multipar_plot <- function(pHSinit=0.2, Kx = 100, betax=0.02, wxtnormm=seq(0,1,0.2), wxtnormsd= 0.3, hx=1, mxtnormm=1,
                          mxtnormsd=0, axtnormm=1, axtnormsd=0, rx=0.1, zxtnormm=seq(0,1,0.2), zxtnormsd= 0.3, gx=4,
                          cx=0.9, phix=0, nseasons=10, nsim=100, HPcut=0.5, pHScut=0.5,maY=100,miY=0, thetax=0.2, Ex=0){

out.multipar <- multipar(pHSinit=pHSinit, Kx = Kx, betax=betax, wxtnormm=wxtnormm, wxtnormsd= wxtnormsd, hx=hx, mxtnormm=mxtnormm,
                         mxtnormsd=mxtnormsd, axtnormm=axtnormm, axtnormsd=axtnormsd, rx=rx, zxtnormm=zxtnormm, zxtnormsd= zxtnormsd, gx=gx,
                         cx=cx, phix=phix, nseasons=nseasons, nsim=nsim, HPcut=HPcut, pHScut=pHScut,maY=maY,miY=miY, thetax=thetax, Ex=Ex)

xvar <- 1-out.multipar$zxtnormm
yvar <- out.multipar$wxtnormm
zvar <- out.multipar$fYLmean

out.multipar.df <- data.frame(xvar, yvar, zvar)

whivarPSlobe80 <- wireframe(zvar ~ xvar * yvar, scales = list(arrows=FALSE, cex= 1.25, col = "black", font = 1,
                                                             distance=c(1,1,1)), screen = list(z =-50, x = -75),
                            xlab = list('Healthy seeds \n selected', rot=-27, cex=1.25, font=2),
                            ylab = list('Disease-conducive \n weather', rot=20, cex=1.25, font=2),
                            zlab = list('Yield loss after 5 seasons (%)', rot=90, cex=1.25, font=2),
                            zlim = range(seq(0, 100,20)), zoom=0.8)

 heatmapplot <- ggplot(out.multipar.df, aes(xvar, yvar, fill=zvar))+
  geom_tile()+
  viridis::scale_fill_viridis(option = "C",
                              guide = guide_colorbar(title = "Yield Loss  After 5 Seasons (%)",
                                                     title.position = "top",
                                                     direction="horizontal",
                                                     barwidth = 20,
                                                     barheight = 2,
                                                     frame.colour = "black")) +
  theme_classic()+
  xlab('Healthy Seeds Selected')+
  ylab('Disease-Conducive Weather')+
  theme(legend.position = "bottom",
        legend.title.align=0.5,
        axis.title = element_text(face = "bold",
                                  size = 20),
        axis.text = element_text(size = 16),
        legend.title = element_text(size = 16,
                                    face="bold"),
        legend.text = element_text(size = 16)
        #legend.background = element_blank(),
        #legend.box.background = element_blank(),
        #panel.grid.major = element_blank(),
        #panel.grid.minor = element_blank(),
        #panel.background = element_rect(fill = "transparent",colour = NA),
        #plot.background = element_rect(fill = "transparent",colour = NA)
  )


#ggsave(filename = "figures/heatmap.png", width = 8, height = 8, bg="transparent", units = "in", dpi = 300, type="cairo-png")
 par(mfrow=c(1,2), mar=c(4,4.4,0.5,0.5), oma=c(0.5,0.5,0.25,0.25))
 print(whivarPSlobe80)
 print(heatmapplot)

}