seedHealth: Scenario Analysis for Seed Health

#' Generate an output matrix for one simulation
#'
#' This function simulates one parameter combination across nseasons once.
#'
#' The truncated random normal variables are generated by \code{\link{altruncnorm}}.
#'
#'
#' @param pHSinit the initial proportion of healthy seed, numeric or numeric vector.
#' @param Kx the total number of plants, positive interger, numeric or numeric vector.
#' @param betax the maximum seasonal transmission rate, numeric or numeric vector.
#' @param wxtnormm the environmental effect on transmission rate (mean of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param wxtnormsd the environmental effect on transmission rate (standard deviation of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param hx the host effect on transmission rate, numeric or numeric vector.
#' @param mxtnormm the vector management effect on transmission rate (mean of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param mxtnormsd the vector management effect on transmission rate (standard deviation of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param axtnormm the roguing effect in terms of decreased DP (mean of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param axtnormsd the roguing effect in terms of decreased DP (standard deviation of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param rx the reversion rate, numeric or numeric vector.
#' @param zxtnormm the proportional selection against diseased plants (mean of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param zxtnormsd the proportional selection against diseased plants (standard deviation of underlying normal distribution prior to truncation), numeric or numeric vector.
#' @param gx the seed production rate in healthy plants, numeric or numeric vector.
#' @param cx the proportional seed production rate in diseased plants, numeric or numeric vector.
#' @param phix the proportion clean seed purchased, numeric or numeric vector.
#' @param nseasons the number of seasons, numeric or numeric vector.
#' @param HPcut the proportion healthy plant number cutoff, numeric or numeric vector.
#' @param pHScut the proportion healthy seed cutoff, numeric or numeric vector.
#' @param maY the maximum attainable yield, end of season, in the absence of disease, numeric or numeric vector.
#' @param miY the minimum yield when all plants are diseased (useable yield despite disease), numeric or numeric vector.
#' @param thetax the rate of decline of Yld with increasing disease incidence, numeric or numeric vector.
#' @param Ex the amount of external inoculum around field, numeric or numeric vector.
#' @keywords seed health
#' @importFrom stats median quantile rnorm var
#' @export
#' @examples
#' onesim() # more examples to be added



# to do - GT

# Columns of output matrix
# col 1 - season timestep (initial time step is season 0)
# col 2 - HP healthy plant number
# col 3 - DP diseased plant number (after roguing)
# col 4 - HS healthy seed number
# col 5 - DS diseased seed number
# col 6 - pHS proportion healthy seed
# col 7 - pDS proportion diseased seed
# col 8 - mx vector management effect on transmission rate
# col 9 - zx proportional selection against diseased plants
# col 10 - ax roguing effect in terms of decreased DP
# col 11 - wx environmental effect on transmission rate
# col 12 - Yld end of season yield
# col 13 - YL end of season yield loss
# col 14 - DPbr diseased plants before roguing

# Columns of output of outfin
# col 15 - 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 16 - 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 17 - HPpseas Healthy plants are calculated from season 1 onwards
# col 18 - pHSpseas Proportion seasons with pHS below pHScut

# Step 1B. Create matrix for output from one simulation (stochastic model)

# Weather (wx), vector management (mx), positive selection (zx) and roguing (zx) are stochastic
# Each have a mean and associated standard deviation

# This function simulates nseasons for one parameter combination once

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


  outm <- as.data.frame(matrix(data=-999, nrow=(nseasons+1), ncol=14, 
                               dimnames = list(1:(nseasons+1), c('season','HP', 'DP','HS','DS','pHS','pDS','mx','zx','ax','wx', 'Yld','YL','DPbr'))))

  outm[1,] <- NA # row one gives initial conditions
  outm$season <- 0:nseasons

  # Initial value for state variables
  outm$pHS[1] <- pHSinit # initial proportion healthy seed (for nseasons=0)
  outm$pDS[1] <- 1 - pHSinit

  # seasons 2 and higher

  for(si in 2:(nseasons+1)) {
    #Generating stochastic  variables
    outm$wx[si] <-ifelse(wxtnormsd==0, wxtnormm, (altruncnorm(1,a=0,b=1,mean=wxtnormm,sd=wxtnormsd)))
    outm$mx[si] <- ifelse(mxtnormsd==0, mxtnormm, (altruncnorm(1,a=0,b=1,mean=mxtnormm,sd=mxtnormsd)))
    outm$ax[si] <- ifelse(axtnormsd==0, axtnormm, (altruncnorm(1,a=0,b=1,mean=axtnormm,sd=axtnormsd)))
    outm$zx[si] <-ifelse(zxtnormsd==0, zxtnormm, (altruncnorm(1,a=0,b=1,mean=zxtnormm,sd=zxtnormsd)))
    #Calculating rate of disease transmission
    tempnewinf<-  betax * outm$wx[si] * hx * outm$mx[si] * ((Kx * outm$pDS[si-1]) * (Kx * outm$pHS[si-1])+ Ex*(Kx * outm$pHS[si-1]))
    #Equation A1
    outm$HP[si] <- min(max(0,Kx * outm$pHS[si-1] - tempnewinf),100)
    # Equation A2
    outm$DPbr[si] <- min(max(0, Kx * outm$pDS[si-1] + tempnewinf),100)
    # Equation A3
    outm$DP[si] <- outm$ax[si]*outm$DPbr[si]
    # Equation B1
    outm$HS[si] <- max(0,gx * (outm$HP[si] + rx * outm$DP[si]))
    # Equation B2
    outm$DS[si] <- max(0,outm$zx[si] * cx * gx * (1 - rx) * outm$DP[si])
    #Equation C1
    outm$Yld[si]<- ((outm$HP[si]+outm$DP[si])/(outm$HP[si]+outm$DPbr[si]))*(miY+(maY-miY)*((1-(outm$DP[si]/(outm$DP[si]+outm$HP[si])))/((1-thetax)+ thetax*(1-(outm$DP[si]/(outm$DP[si]+outm$HP[si]))))^2))
    # Equation C2
    outm$YL[si]<- maY-outm$Yld[si]
    # Equation D1
    outm$pHS[si] <- phix + (1-phix) * outm$HS[si]/(outm$HS[si] + outm$DS[si])
    # Equation D2
    outm$pDS[si] <- 1 - outm$pHS[si]
  }

  outfin <- outm[(nseasons+1),]

  # Season in which HP first transitions below HPcut*Kx
  HPtrans <- outm$season[which(outm$HP < HPcut*Kx)][1]
  # if HPtrans is NA i.e., HP never less than HPcut, set to max seasons tested
  HPtrans[is.na(HPtrans)]<-nseasons
  # Season in which pHS first transitions below pHScut
  pHStrans <- outm$season[which(outm$pHS < pHScut)][1]
  # if pHStrans is NA i.e., pHS never less than pHScut, set to max seasons tested
  pHStrans[is.na(pHStrans)]<-nseasons
  # Proportion seasons with HP below HPcut*Kx
  # Healthy plants are calculated from season 1 onwards so +1 not needed in denominator
  HPpseas <- length(outm$season[which(outm$HP < HPcut*Kx)])/nseasons
  # Proportion seasons with pHS below pHScut
  # if +1 is not included in denominator, if in season 0 pHS<pHScut, the proportion is >1
  pHSpseas <- length(outm$season[2:nseasons][which(outm$pHS < pHScut)])/(nseasons+1)

  # ptempa <- matrix(c(HPtrans,pHStrans,HPpseas,pHSpseas), dimnames=list(c('HPtrans','pHStrans','HPpseas','pHSpseas'),'1'))
  # ptempb <- as.data.frame(t(ptempa))

  ptemp <- as.data.frame(t(matrix(c(HPtrans,pHStrans,HPpseas,pHSpseas), dimnames=list(c('HPtrans','pHStrans','HPpseas','pHSpseas'),'1'))))

  outfin <- cbind(outfin,ptemp)

  #-----------------------------------------------
  # warning message

  if ( pHSinit < 0 | pHSinit > 1){
    warning(paste('pHSinit: your input value is', pHSinit,', it must be between 0 and 1'))
  } else if (betax < 0.001) {
    warning(paste('betax: your input value is', betax,', it must be greater than 0.001'))
  } else if (wxtnormm < 0 | wxtnormm > 1) {
    warning(paste('wxtnormm: your input value is', wxtnormm,', it must be between 0 and 1'))
  } else if (hx < 0 | hx > 1) {
    warning(paste('hx: your input value is', hx,', it must be between 0 and 1'))
  } else if (mxtnormm < 0 | mxtnormm > 1) {
    warning(paste('mxtnormm: your input value is', mxtnormm,', it must be between 0 and 1'))
  } else if (axtnormm < 0 | axtnormm > 1) {
    warning(paste('axtnormm: your input value is', axtnormm,', it must be between 0 and 1'))
  } else if (rx < 0 | rx > 1) {
    warning(paste('rx: your input value is', rx,', it must be between 0 and 1'))
  } else if (zxtnormm < 0 | zxtnormm > 1) {
    warning(paste('zxtnormm: your input value is', zxtnormm,', it must be between 0 and 1'))
  } else if (gx < 0) {
    warning(paste('gx: your input value is', gx,', it must be 0 or a positive integer'))
  } else if (cx < 0 | cx > 1) {
    warning(paste('cx: your input value is', cx,', it must be between 0 and 1'))
  } else if (phix < 0 | phix > 1) {
    warning(paste('phix: your input value is', phix,', it must be between 0 and 1'))
  } else if (thetax < -1 | thetax > 1) {
    warning(paste('thetax: your input value is', thetax,', it must be between -1 and 1'))
  } else if (Ex < 0) {
    warning(paste('Ex: your input value is', Ex,', it must be greater than 0'))
  }else {
    list(outm=outm,outfin=outfin)
  }
  #-----------------------------------------------

}