R/plotting.R

In luke-a-rogers/pbsedm: An R package to implement some of the methods of empirical dynamic modelling

##' Obtain axes ranges (adapted from function by Uwe Ligges)
##'
##' Obtain axes ranges for scatterplot3d plot (and maybe does others), adapted
##' from http://r.789695.n4.nabble.com/Axis-Limits-in-Scatterplot3d-td892026.html<description>
##'
##' @param s3dobject scatterplot3d plot presumably
##' @return Ranges of x, y and z axes
##' @export
##' @author Andrew Edwards (adapted from Uwe Ligges' function)
##' @examples
##' \donttest{
##'   s3d <- scatterplot3d::scatterplot3d(rnorm(5), rnorm(5), rnorm(5))
##'   gets3dusr(s3d)
##' }
gets3dusr = function(s3dobject)
  {
  es3d = environment(s3dobject$points3d)
  g3d  = function(name) get(name, envir=es3d)
  c(c(g3d("x.min"),
      g3d("x.max")) * g3d("x.scal"),
    c(0,
      g3d("y.max")) * g3d("y.scal") + g3d("y.add"),
    c(g3d("z.min"),
      g3d("z.max")) * g3d("z.scal"))
}


##' Plot data and results of an object of class `pbsEDM` as time series and
##' phase plots
##'
##' Plots time series of `N_t` and `Y_t`, 2d phase plots of `N_t` vs `N_{t-1}`
##' and `Y_t` vs `Y_{t-1}`, and 3d phase plot of `Y_t` vs `Y_{t-1}` vs `Y_{t-2}`.
##' See vignette "Analyse a simple time series".
##' Only working if `N$N_t` values are in there (needs another switch if
##' not). And very likely only works for univariate time series for now, not
##' pred-prey for example.
##'
##' @param x list object of class `pbsEDM`, an output from `pbsEDM()`
##' @param portrait if TRUE then plots panels in portrait mode for manuscripts/Rmarkdown (3x2), false is
##'   landscape (2x3, filled columnwise) for presentations.
##' @param ... additional arguments to be passed onto other plotting functions,
##'   in particular `last.time.to.plot` to plot only up to that time step (to
##'   loop through in a movie), and late.num to plot the final number of time
##'   steps in a different colour, or label.cex for (a) etc. label sizes.
##' @return Five panels in 2x3 or 3x2 format, in current plot environment
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'  aa <- pbsEDM(NY_lags_example,
##'               lags = list(N_t = 0:1),
##'               first_difference = TRUE)
##'  plot(aa)
##' }
plot.pbsEDM = function(x,
                       portrait = TRUE,
                       ...){
  stopifnot(attr(x, "class") == "pbsEDM")

  ifelse(portrait,
         par(mfrow = c(3,2)),
         par(mfcol = c(2,3)))

  # Commenting as isn't passed on anyway
  #if(!exists("last.time.to.plot")) last.time.to.plot <- length(x$X_observed)
  #                                          # though last will have NA # not
  #                                          # incorporated fully yet

  plot_time_series(values = x$N,  # $N_t
                   X.or.N = "N",
                   label = "(a)",
                   ...)

  plot_time_series(values = x$X_observed,
                   X.or.N = "X",
                   label = "(c)",
                   ...)

  plot_phase_2d(values = x$N,    # $N_t,
                X.or.N = "N",
                label = "(b)",
                ...)

  plot_phase_2d(values = x$X_observed,
                X.or.N = "X",
                label = "(d)",
                ...)

  plot_phase_3d(x,
                label = "(e)",
                ...)
  invisible()
}

##' Plot output from `pbsEDM_Evec()` as six panel plot, with data shown the
##'  same as for `plot.pbsEDM()` and sixth panel with predictions vs observation
##'  for each `E`
##'
##' Plots time series of `N_t` and `Y_t`, 2d phase plots of `N_t` vs `N_{t-1}`
##' and `Y_t` vs `Y_{t-1}`, and 3d phase plot of `Y_t` vs `Y_{t-1}` vs
##' `Y_{t-2}`. Sixth panel shows predictions vs observations for different values
##' of `E`. See vignette "Analyse a simple time series" for full details.
##'
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param ... extra arguments to `plot.pbsEDM()`, such as label.cex
##' @return Six-panel figure in current plot environment
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   aa_Evec <- pbsEDM_Evec(NY_lags_example$N_t)
##'   plot_pbsEDM_Evec(aa_Evec, portrait = FALSE)
##'
##' # For manuscript figures:
##' E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##'
##' postscript("six_panels_all.eps",
##'            height = 5.36,
##'            width = 9,
##'            horizontal=FALSE,
##'            paper="special")
##' plot_pbsEDM_Evec(E_results,
##'     last.time.to.plot = NULL, # to automatically plot all
##'     portrait = FALSE)
##' dev.off()
##' }
plot_pbsEDM_Evec <- function(E_res,
                             ...){
  # The data are the same for each value of E, and so use one for the first five panels:
  plot.pbsEDM(E_res[[1]],
              ...)

  plot_pred_obs(E_res,
                label = "(f)",
                ...)
  invisible()
}


##' Plot output from `plot_pbsEDM_Evec()` as six panel plot saved as
##'  encapsulated postscript file for manuscript.
##'
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param eps.filename filename to save to
##' @param ... extra arguments to `plot_pbsEDM_Evec()`.
##' @return Saved .eps file to use in manuscript
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##' E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##' plot_pbsEDM_Evec_save(E_results)
##' }
plot_pbsEDM_Evec_save <- function(E_res,
                                  eps.filename = "six_panels_all.eps",
                                  ...){

  postscript(eps.filename,
             height = 5.36,
             width = 9,
             horizontal=FALSE,
             paper="special")

  plot_pbsEDM_Evec(E_res,
                   last.time.to.plot = NULL, # to automatically plot all
                   portrait = FALSE,
                   ...)

  dev.off()
}

##' Plot rho versus `E` for output from `pbsEDM_Evec()`
##'
##' Plot correlation coefficient against `E`.
##'
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param ... Extra arguments to pass to `plot()`
##' @return plot of rho versus E
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##' E_results_17 <- pbsEDM_Evec(NY_lags_example$N_t,
##'                             E_vec = 1:16)  # go up to E=17
##'   postscript("rho_Evec.eps",
##'             height = 4.5,
##'             width = 4.5,
##'             horizontal=FALSE,
##'             paper="special")
##'   plot_rho_Evec(E_results_17)
##'   dev.off()
##' }
plot_rho_Evec <- function(E_res,
                          ...){
  rho <- vector()
  E <- vector()
  for(i in 1:length(E_res)){
    rho[i] <- E_res[[i]]$results$X_rho
    E[i] <- E_res[[i]]$results$E
  }

  plot(E,
       rho,
       xlim = c(0, max(E)),
       ylim = c(0,1),
       xlab = expression("Embedding dimension," * italic(E)),
       ylab = expression("Forecast skill, " * rho),
       ...
       )
  box()
  axis(1, 1:max(E),
       tcl = -0.2,
       labels = rep("", max(E)))
  axis(2,
       seq(0, 1, 0.1),
       tcl = -0.2,
       labels = rep("", 11))

  # see EDMsimulate/report/sockeye-sim-edm.rnw for adding other plots
  invisible()
}

##' Save `rho` versus `E` .eps figure for manuscript.
##'
##' @param eps.filename filename to save to
##' @return plot of rho versus E saved as .eps
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##' }
plot_rho_Evec_save <- function(eps.filename = "rho_Evec.eps"){

  # Use high E values to show the decrease
  E_results_17 <- pbsEDM_Evec(NY_lags_example$N_t,
                              E_vec = 1:16)  # go up to E=17

  postscript(eps.filename,
             height = 4.5,
             width = 4.5,
             horizontal=FALSE,
             paper="special")

  plot_rho_Evec(E_results_17)

  dev.off()
}


##' Plot predictions versus observed values of `Y_t` for list of `pbsEDM`
##' objects, for various values of `E`
##'
##' Plot the predicted versus observed values of `Y_t` for various values
##' of `E`, using the output (a list of `pbsEDM` objects) from
##' `pbsEDM_Evec()`. Each value of `t` corresponds to focal point `t* = t - 1`.
##'
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param E_components How many of the first E_res components to show
##' @param E_cols Vector of colours, one for each value of E
##' @param last.time.to.plot Last time value of `t` to use when plotting.
##' @param E_cex Vector of sizes, one for each value of E
##' @param portrait dummy argument that allows `...` to be passed from
##'   `plot_pbsEDM_Evec()`, from which we want `last.time.to.plot`.
##' @param label label to annotate plot, such as `(a)` etc. for six-panel figure
##' @param label.cex size of label annotation
##' @return Figure in the current plot enivironment
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   aa <- pbsEDM_Evec(NY_lags_example$N_t)
##'   plot_pred_obs(aa)
##' }
plot_pred_obs <- function(E_res,
                          E_components = 5,
                          E_cols = c("orange",
                                     "blue",
                                     "green",
                                     "red",
                                     "black"),
                          last.time.to.plot = NULL,
                          E_cex = seq(1.7, 0.6, length = 5),
                          portrait = NULL,
                          label = NULL,
                          label.cex = 0.7
                          ){

  stopifnot("Need more distinct colours in E_cols"=
              length(E_cols) <= E_components)

  stopifnot("Need more values in E_cex"=
              length(E_cex) <= E_components)

  if(is.null(last.time.to.plot)) last.time.to.plot <- length(E_res[[1]]$X_observed)
  if(is.null(portrait)){portrait = NULL} # to use the dummy argument portrait
  # Determine range for axes
  obs.max.abs = max( abs( range(E_res[[1]]$X_observed,
                                na.rm=TRUE) ) ) # max abs observed value
  forecast.max.abs = 0         # max abs forecast value
  for(j in 1:E_components){
    forecast.max.abs = max(c(forecast.max.abs,
                             abs( range( range(E_res[[j]]$X_forecast,
                                               na.rm=TRUE) ) ) ) )
  }

  max.abs = max(obs.max.abs, forecast.max.abs)
  axes.range = c(- max.abs, max.abs)

  plot(0, 0,
       xlab = expression("Observation of " * italic(Y) [t]),
       ylab = expression("Prediction of " * italic(Y) [t]),
       xlim = axes.range,
       ylim = axes.range,
       asp = 1,
       type = "n")

  if(!is.null(label)){
    mtext(text = label,
          at = axes.range[1],
          adj = 1,
          cex = label.cex)
  }

  abline(0, 1, col="grey")
  leg = vector()

  for(j in 1:E_components){
    points(E_res[[j]]$X_observed[1:last.time.to.plot],
           E_res[[j]]$X_forecast[1:last.time.to.plot],
           pch = 16, # could note the final one differently (but not a star,
                     # since that's last N[t] not Y[t])
           col = E_cols[j],
           cex = E_cex[j])
    leg = c(leg,         # Could try and make E italic but likely fiddly
            as.expression(bquote(italic(E) * "=" *.(E_res[[j]]$results$E) *
    ", " * rho * "=" *.(format(round(E_res[[j]]$results$X_rho,
                                2),
                          nsmall = 2)))))
  }

  legend("topleft",
         pch = 16,
         leg,
         col = E_cols,
         pt.cex = E_cex)
  invisible()
}

##' Plot the observed time series as either `N_t` or `Y_t`
##'
##' First value must be `t=1`. For non-differenced values `N_t`, shows the time
##' series with the final values in a different colour, and a title showing the
##' final time step. For first-differenced values `Y_t`, shows the time series
##' of those, plus a one-dimensional phase plot.
##'
##' @param values vector of values to be plotted
##' @param X.or.N "N" if raw non-differenced data, "X" for differenced data
##' @param par.mar.ts `par(mar)` values
##' @param max_time maximum time value for the time axis
##' @param t.axis.range range of time axis
##' @param last.time.to.plot last time value of N[t] to use when plotting, so
##'   final Y[t] used will be Y[t-1] (since Y[t] uses N[t+1])
##' @param late.num final number of `N[t]` time steps to plot in a different colour
##' @param late.col colour in which to plot final `late.num` time steps
##' @param early.col colour in which to plot earlier time step points
##' @param early.col.lines colour in which to plot earlier time step points
##' @param start first time step (must be 1)
##' @param pt.type `type` value for `points()`
##' @param par.mgp `par("mgp")` values
##' @param add.legend logical, whether to add legend to the `N[t]` time series
##' @param label label to annotate plot, such as `(a)` etc. for six-panel figure
##' @param label.cex size of label annotation
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   plot_time_series(NY_lags_example$N_t, X.or.N = "N")
##'   plot_time_series(NY_lags_example$Y_t, X.or.N = "X")
##' }
plot_time_series <- function(values,
                             X.or.N,
                             par.mar.ts = c(3, 3, 1, 1),
                             max_time = NULL,
                             t.axis.range = NULL,
                             last.time.to.plot = NULL,
                             late.num = 3,
                             late.col = "red",
                             early.col = "black",
                             early.col.lines = "lightgrey",
                             start = 1, # may not work for others
                             pt.type = "p",
                             par.mgp = c(1.5, 0.5, 0),
                             add.legend = TRUE,
                             label = NULL,
                             label.cex = 0.7
                             ){
  stopifnot(start == 1)

  if(is.null(max_time)) max_time <- length(values)

  if(is.null(last.time.to.plot)){
    if(X.or.N == "N"){
       last.time.to.plot <- max(which(!is.na(values)))} else {
                                                      last.time.to.plot <-
                                                        length(values)
                                                      # Want the NA included for
                                                      #  X, but N has an extra
                                                      #  time step added
                                                      }
  }


  if(is.null(t.axis.range)) {
      t.axis.range <- c(start, max_time)
  }

  par("mgp" = par.mgp) # first val sets axis title dist (mex units)
                       #  to axes, second sets labels
  par(pty = "m")                 # maximal plotting region, not square
                                 #  like for phase plots
  par(mar = par.mar.ts)

  col.plot = c(rep(early.col,
                   max(c(0, last.time.to.plot - late.num))),
               rep(late.col,
                   min(c(last.time.to.plot, late.num))) )   # colours of points
  col.plot.lines = col.plot                            # colours of lines
  col.plot.lines[col.plot.lines == early.col] = early.col.lines
  pch.plot = (col.plot == early.col) * 1 + (col.plot == late.col) * 16
                                        # filled circles for latest
  pch.plot[length(pch.plot)] = 8     # latest one a star
  pch.plot[length(pch.plot) - late.num + 1] = 1 # late.num'th from end an open circle


  if(X.or.N == "N"){
    Nt.max.abs = max( abs( range(values[start:max_time],
                                 na.rm=TRUE) ) )
    Nt.axes.range = c(0, Nt.max.abs*1.04)    # Expand else points can hit edge

    plot(0, 0,
         xlab = expression("Time, " * italic(t)),
         ylab = expression(italic(N) [t]),
         xlim = c(0, max(t.axis.range)),
         ylim = Nt.axes.range,
         type = "n",                           # empty plot
         main = as.expression(bquote("Time " * italic(t) * "=" *.(last.time.to.plot))))

    if(add.legend){
      # for the red late.num points
      late.num.ind = (length(pch.plot) - late.num + 1):length(pch.plot)
      late.num.ind = late.num.ind[late.num.ind > 0]   # else t=0 shown for t=2 plot
      legend("topright",
             pch = pch.plot[late.num.ind],
             col = col.plot[late.num.ind],
             leg = paste0("t=",
                          late.num.ind),
             cex = 1,
             bty = "n")
    }

    if(!is.null(label)){
      mtext(text = label,
            at = 0,
            adj = 1,
            cex = label.cex)
    }

  } else {
    XtLoc = -0.05 * max(t.axis.range)  # location to plot Xt on a vertical line,
    Xt.max.abs = max(abs( range(values[start:max_time],
                                na.rm=TRUE) ),
                     abs( range(values[start:max_time],
                                na.rm=TRUE)) )
    Xt.axes.range = c(-Xt.max.abs, Xt.max.abs)

    plot(0, 0,
         xlab = expression("Time, " * italic(t)),
         ylab = expression(italic(Y) [t] ~ "=" ~ italic(N) [t+1] ~ "-" ~ italic(N) [t]),
         xlim = c(XtLoc, max(t.axis.range)),
         ylim = Xt.axes.range,
         type = "n")                           # empty plot
    abline(v = 0.5*XtLoc, col="black")

    if(!is.null(label)){
      mtext(text = label,
            at = XtLoc,
            adj = 1,
            cex = label.cex)
    }
  }

  iii = last.time.to.plot             # use iii since simpler
  if(iii > 1.5){
    segments(start:(iii-1),
             values[start:(iii-1)],    # dplyr::pull(Nx.lags.use[start:(iii-1), "Y_t"]),
             (start+1):iii,
             values[(start+1):iii],
             col = col.plot.lines)      # lines() will not use vector col
  }

  points(start:iii,
         values[start:iii],
         type = pt.type,
         pch = pch.plot,
         col = col.plot)

  # '1d phase plot':
  if(X.or.N == "X"){
    points(rep(XtLoc, iii-start+1),
           values[start:iii],
           type = pt.type,
           pch = pch.plot,
           col = col.plot)
  }
  invisible()
}

##' Plot 2d phase plot of `N_t` vs `N_{t-1}` or `Y_t` vs `Y_{t-1}`
##'
##' Shows the values in two-dimensional phase space (first differenced or not)
##'   with lag of 1. Cobwebbing shows how one point iterates to the next point
##'   in the time series.
##'
##' @param values vector of `N_t` or `Y_t` values
##' @param X.or.N "N" if raw non-differenced data, "X" for differenced data
##' @param par.mar.phase `par(mar)` values
##' @param axis.range range of axes, if NA then calculated from values
##' @param start first time step to plot (currently must be 1)
##' @param last.time.to.plot last time value of N[t] to use when plotting, so
##'   final Y[t] used will be Y[t-1] (since Y[t] uses N[t+1])
##' @param pt.type `type` value for `points()`
##' @param cobwebbing if TRUE then add cobwebbing lines to phase plot
##' @param late.col colour in which to plot final `late.num` time steps
##' @param early.col colour in which to plot earlier time step points
##' @param early.col.lines colour in which to plot earlier time step points
##' @param late.num final number of `N[t]` time steps to plot in a different colour
##' @param label label to annotate plot, such as `(a)` etc. for six-panel figure
##' @param label.cex size of label annotation
##' @return Plots figure to current device
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   plot_phase_2d(NY_lags_example$N_t, X.or.N = "N")
##'   plot_phase_2d(NY_lags_example$Y_t, X.or.N = "X")
##' }
plot_phase_2d <- function(values,
                          X.or.N = "X",
                          par.mar.phase = c(3, 0, 1, 0),
                          axis.range = NA,
                          start = 1,
                          last.time.to.plot = NULL,
                          pt.type = "p",
                          cobwebbing = TRUE,
                          late.col = "red",
                          early.col = "black",
                          early.col.lines = "lightgrey",
                          late.num = 3,
                          label = NULL,
                          label.cex = 0.7
                          ){

  if(is.na(axis.range)) {
    expand <- ifelse(X.or.N == "N", 1.04, 1)   # expand N axis as in plot_time_series()
    axis.range <- c(min(0, min(values, na.rm = TRUE)),
                    max(values, na.rm = TRUE) * expand)   # for N_t or Y_t
  }
  if(is.null(last.time.to.plot)){
    if(X.or.N == "N"){
      last.time.to.plot <- max(which(!is.na(values)))} else {
                                                       last.time.to.plot <-
                                                        length(values)
                                                       # Want the NA included for
                                                       #  X, but N has an extra
                                                       #  time step
                                                       #  added. i.e. want 100 for both
                                                     }
  }

  stopifnot(start == 1)
  # Now copying from plot_observed:
  col.plot = c(rep(early.col,
                   max(c(0, last.time.to.plot - late.num))),
               rep(late.col,
                   min(c(last.time.to.plot, late.num))) )   # colours of points
  col.plot.lines = col.plot                            # colours of lines
  col.plot.lines[col.plot.lines == early.col] = early.col.lines
  pch.plot = (col.plot == early.col) * 1 + (col.plot == late.col) * 16
                                        # filled circles for latest
  pch.plot[length(pch.plot)] = 8     # latest one a star
  pch.plot[length(pch.plot) - late.num + 1] = 1 # late.num'th from end an open circle

  # Copying from plotPanelMovie.df2

  par(pty="s")             # set following plot types to be square
                           #  (without this the axes don't seem to
                           #  be the same, even with the settings below)
  par(mar = par.mar.phase) # margins

  # N_t vs N{t-1}:
  # Empty plot to get started
  if(X.or.N == "N"){
    plot(0, 0,
         xlab = expression(italic(N) [t-1]),
         ylab = expression(italic(N) [t]),
         xlim = axis.range,
         ylim = axis.range,
         type = "n")
    values.to.plot <- values[start:last.time.to.plot]
  } else {
    plot(0, 0,
         xlab = expression(italic(Y) [t-1]),
         ylab = expression(italic(Y) [t]),
         xlim = axis.range,
         ylim = axis.range,
         type = "n")
    values.to.plot <- values[start:last.time.to.plot]
  }

  if(!is.null(label)){
    mtext(text = label,
          at = axis.range[1],
          adj = 1,
          cex = label.cex)
  }

  if(cobwebbing) abline(0, 1, col="darkgrey")

  # Draw lines first so they get overdrawn by points
  if(last.time.to.plot > 2.5){
    if(cobwebbing){
      # Do lines for cobwebbing
      vals = rep(values.to.plot,
                 each = 2)
      vals = vals[-1]
      vals = vals[-length(vals)]
      len = length(vals)

      col.cobweb.lines = rep(early.col.lines, len)
      col.cobweb.lines[(max(len - 2*late.num + 1,
                            1)):len] = late.col
      # Take off two more for "X" but only if final one is NA
      segments(vals[1:(len - 2 - 2*(X.or.N == "X" & is.na(vals[len])))],
               vals[2:(len - 1 - 2*(X.or.N == "X" & is.na(vals[len])))],
               vals[2:(len - 1 - 2*(X.or.N == "X" & is.na(vals[len])))],
               vals[3:(len - 2*(X.or.N == "X" & is.na(vals[len])))],
               col = col.cobweb.lines)
    } else {
      # Join each point to the next   (N(5), N(6)) to (N(6), N(7))
      #  Not checked.
      segments(
        # dplyr::pull(Nx.lags.use[start:(iii-1), "N_tmin1"]),
        # dplyr::pull(Nx.lags.use[start:(iii-1), "N_t"]),
        # dplyr::pull(Nx.lags.use[(start+1):iii, "N_tmin1"]),
        # dplyr::pull(Nx.lags.use[(start+1):iii, "N_t"]),
        # not fully tested:
        pbsLag(values.to.plot)[-length(values.to.plot)],   # N_{t-1}
        values.to.plot[-length(values.to.plot)],
        pbsLag(values.to.plot)[-1],
        values.to.plot[-1],
        col = col.plot.lines) # lines() will not use vector of col
    }
  }
  if(last.time.to.plot > 1.5){
    points(pbsLag(values.to.plot),
           values.to.plot,
           type = pt.type,
           pch = pch.plot,
           col = col.plot)          # start row has NA's, gets ignored
  }
      # legend("topright", legend=paste0("Time t=", iii), box.col = "white",
      #        inset = 0.01)  # inset to stop white overwriting outer box


      # x_t vs x{t-1}:

#      if(iii > 2.5)
#        {
#          if(cobwebbing)
#            {
#               xvals = rep( dplyr::pull(Nx.lags.use[start:iii, "Y_t"]), each = 2)
#               xvals = xvals[-1]
#               xvals = xvals[-length(xvals)]
#               lenx = length(xvals)
#              segments(xvals[1:(lenx-2)],
#                        xvals[2:(lenx-1)],
#                       xvals[2:(lenx-1)],
#                        xvals[3:lenx],
#                        col = col.cobweb.lines)
#           } else
#           {  # Just join consecutive points with lines
#               segments(dplyr::pull(Nx.lags.use[start:(iii-1), "Y_tmin1"]),
#                        dplyr::pull(Nx.lags.use[start:(iii-1), "Y_t"]),
#                        dplyr::pull(Nx.lags.use[(start+1):iii, "Y_tmin1"]),
#                        dplyr::pull(Nx.lags.use[(start+1):iii, "Y_t"]),
#                        col = col.plot.lines)
#            }
#        }
#      if(iii > 1.5)
#        {
#          points(dplyr::pull(Nx.lags.use[start:iii, "Y_tmin1"]),
#                 dplyr::pull(Nx.lags.use[start:iii, "Y_t"]),
#                 type = pt.type, pch = pch.plot,
#                 col = col.plot)           # start row has NA's, get ignored
#        }
      # legend("topleft", legend=paste("Time", iii), border = NULL)
  invisible()
}

##' Plot 3d phase plot of `Y_t` vs `Y_{t-1}` vs `Y_{t-2}`
##'
##' Shows the first-differenced values in three-dimensional phase space, as
##'   points showing lags of 0, 1 and 2. Cobwebbing shows how one point iterates
##'   to the next point in the time series. Can also highlight a $t^*$ value and
##'   the points that should be excluded from the library of nearest neighbours.
##'
##' @param obj `pbsEDM` object
##' @param par.mgp.3d `par(mgp)` values
##' @param par.mai.3d `par(mai)` values
##' @param par.mar.3d `par(mar)` values
##' @param x.lab x axis label
##' @param y.lab y axis label
##' @param z.lab z axis label
##' @param last.time.to.plot last time value of N[t] to use when plotting, so
##'   final Y[t] used will be Y[t-1] (since Y[t] uses N[t+1])
##' @param axis.range range of axes, if NA then calculated from values; all
##'   three axes have the same range
##' @param late.num final number of `N_t` time steps to plot in a different
##'   colour; not showing `N_t` here but keeping consistency with other plots
##' @param pt.type `type` value for `points()`
##' @param late.col colour in which to plot final `late.num` time steps
##' @param early.col colour in which to plot earlier time step points
##' @param early.col.lines colour in which to plot earlier time step points
##' @param axes.col colour of orthogonal origin axes that go through (0, 0, 0)
##' @param par.mar.phase `par(mar)` to reset to after plotting 3d figure
##' @param par.mgp `par(mgp)` to reset to after plotting 3d figure
##' @param tstar focal point to highlight (for `pbsSmap` vignette)
##' @param angle_view manually specify the angle for viewing (to rotate plot)
##' @param label label to annotate plot, such as `(a)` etc. for six-panel figure
##' @param label.cex size of label annotation
##' @return Plots figure to current device
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'  aa <- pbsEDM(NY_lags_example,
##'               lags = list(N_t = 0:2),
##'               first_difference = TRUE)
##'  plot_phase_3d(aa)
##' }
plot_phase_3d <- function(obj,
                          par.mgp.3d = c(3, 10, 0),
                          par.mai.3d = c(0.1, 0.1, 0.1, 0.1),
                          par.mar.3d = c(3, 0, 0, 0),
                          x.lab = expression(italic(Y) [t-2]),
                          y.lab = expression(italic(Y) [t-1]),
                          z.lab = expression(italic(Y) [t]),
                          last.time.to.plot = NULL,
                          axis.range = NA,
                          late.num = 3,
                          pt.type = "p",
                          late.col = "red",
                          early.col = "black",
                          early.col.lines = "lightgrey",
                          axes.col = "darkblue",
                          par.mar.phase = c(3, 0, 1, 0),  # to reset for normal figs
                          par.mgp = c(1.5, 0.5, 0),
                          tstar = NA,
                          angle_view = NA,
                          label = NULL,
                          label.cex = 0.7
                          ){
  if(is.na(axis.range)) {
    axis.range <- c(min(0, min(obj$X_observed, na.rm = TRUE)),
                    max(obj$X_observed, na.rm = TRUE))
  }

  if(is.null(last.time.to.plot)) last.time.to.plot <- length(obj$X_observed)

  if(is.na(angle_view)){
    if(is.na(tstar)){
      angle_view <- 40 + last.time.to.plot} else {
                                            angle_view <- 40 + tstar
                                          } # to match movie
  }


  if(!is.na(tstar)) late.num = 1

  start = 1     # currently only works for 1

  if(last.time.to.plot > 2){
    values.to.plot <- obj$X_observed[start:last.time.to.plot]
  } else {
    values.to.plot <- NA          # nothing to plot
  }

  # Copied from plot_observed:
  col.plot = c(rep(early.col,
                   max(c(0, last.time.to.plot - late.num))),
               rep(late.col,
                   min(c(last.time.to.plot, late.num))) )   # colours of points
  col.plot.lines = col.plot                            # colours of lines
  col.plot.lines[col.plot.lines == early.col] = early.col.lines
  pch.plot = (col.plot == early.col) * 1 + (col.plot == late.col) * 16
                                        # filled circles for latest
  pch.plot[length(pch.plot)] = 8     # latest one a star
  pch.plot[length(pch.plot) - late.num + 1] = 1 # late.num'th from end an open circle

  # Empty plot to get started
  par(mgp = par.mgp.3d)
  par(mai = par.mai.3d)  # scat..3d resets mar, think mai still has an effect
  par.mar.3d = c(3, 0, 0, 0)
  scat = scatterplot3d::scatterplot3d(0,
                                      0,
                                      0,
                                      xlab = x.lab,
                                      ylab = y.lab,
                                      zlab = z.lab,
                                      xlim = axis.range,
                                      ylim = axis.range,
                                      zlim = axis.range,
                                      type = "n",
                                      box = FALSE,
                                      angle = angle_view,
                                      mar = par.mar.3d)
                                              # mar=c(5,3,4,3)+0.1 is default,set
                                              #  within scatterplot3d
                                              # Add axes so can see origin:
      actual.axes.ranges = gets3dusr(scat)    # Get the actual values used
      scat$points3d(actual.axes.ranges[1:2], c(0,0), c(0,0), type = "l",
                    col = axes.col)
      scat$points3d(c(0,0), actual.axes.ranges[3:4], c(0,0), type = "l",
                    col = axes.col)
      scat$points3d(c(0,0), c(0,0), actual.axes.ranges[5:6], type = "l",
                    col = axes.col)
      # Obtain x-y co-ords of points for segments:
#**      proj.pts = scat$xyz.convert(dplyr::pull(Nx.lags.use[start:iii,
#                                                            "Y_tmin2"]),
#                                    dplyr::pull(Nx.lags.use[start:iii,
#                                                            "Y_tmin1"]),
  #                                    dplyr::pull(Nx.lags.use[start:iii, "Y_t"]) )
  # maybe this is okay with NA's:
#  proj.pts = scat$xyz.convert(pbsLag(obj$X_observed,
#                                     2)[start:iii],  # "Y_tmin2", index wrong?
#                              pbsLag(obj$X_observed,
#                                     1)[start:iii],  # "Y_tmin1"
  #                              pbsLag(obj$X_observed)[start:iii])  # "Y_t"

  if(!is.null(label)){
    legend("topleft",    # Fake legend ensures it's fixed in the animation
           bty = "n",
           legend = label,
           cex = label.cex * 1.5) # Needs to be larger for 3-d
  }

  if(!all(is.na(values.to.plot))){
    proj.pts = scat$xyz.convert(pbsLag(values.to.plot,
                                       2),  # "Y_tmin2"
                                pbsLag(values.to.plot,
                                       1),  # "Y_tmin1"
                                values.to.plot)  # "Y_t"
    if(last.time.to.plot > 3.5){
      segments(proj.pts$x[1:(last.time.to.plot - start)],
               proj.pts$y[1:(last.time.to.plot - start)],
               proj.pts$x[2:(last.time.to.plot - start + 1)],
               proj.pts$y[2:(last.time.to.plot - start + 1)],
               col = col.plot.lines) # lines() will not use vector
    }

    # The points
    if(last.time.to.plot > 2.5){
      scat$points3d(pbsLag(values.to.plot,
                           2),  # "Y_tmin2"
                    pbsLag(values.to.plot,
                           1),  # "Y_tmin1"
                    values.to.plot,  # "Y_t"
                    type = pt.type,
                    pch = pch.plot,
                    col = col.plot)
    }
  }

  # Points to highlight for pbsSmap vignette, tstar and subsequent points that
  #  are excluded from library as candidate neighours. Need the time increments.
  if(!is.na(tstar)){
    scat$points3d(obj$X[tstar:(tstar + 3), "Y_t_2"], # "Y_tmin2"
                  obj$X[tstar:(tstar + 3), "Y_t_1"], # "Y_tmin1"
                  obj$X[tstar:(tstar + 3), "Y_t_0"], # "Y_t"
                  type = "h",
                  pch = 16,
                  col = c("blue", rep("red", 3)))
  }

  par(mar = par.mar.phase)   # scatterplot3d changes mar
  par(mgp = par.mgp)         #  so set back to usual for 2d figures
  invisible()                   # else returns the above line
}

##' 2-d phase plot of x_t v x_{t-1} with coloured points to explain EDM for E=2
##'
##' Highlights a point to be projected, its nearest `E+1` neighbours, and then
##' draw arrows to show where they go and so where the projection goes. Very
##' useful for understanding and checking what EDM is doing. Type
##' `plot_explain_edm_movie` for example calls, and use that function to save
##' the movie. Use `plot_explain_edm_save()` for single figure for manuscript.
##'
##' @param obj list object of class `pbsEDM`, an output from `pbsEDM()`
##' @param tstar the time index, `t*` for `x(t*)` of the target point for which to make a
##'   projection from
##' @param x.lab label for x axis
##' @param y.lab label for x axis
##' @param main main title
##' @param tstar.col colour for `x(t*)`
##' @param tstar.pch pch value (point style) for `x(t*)`
##' @param tstar.cex cex value (size) for `x(t*)`
##' @param lib.remove if TRUE then shows which values to exclude from library
##' @param neigh.plot if TRUE then highlight neighbours to `x(t*)`
##' @param neigh.proj if TRUE then highlight projections of neighbours to `x(t*)`
##' @param pred.plot if TRUE then highlight forecasted value `x(t*+1)`
##' @param pred.rEDM if TRUE then show predictions from `rEDM` for
##'   `NY_lags_example` saved values; obj must be `NY_lags_example`
##' @param true.val if TRUE then plot the true value of `x(t*+1)`
##' @param legend.plot if TRUE then do a legend and print value of `t*`
##' @param legend.inc.rEDM if TRUE then include `rEDM` in the legend (not wanted
##'   for first manuscript Figure).
##' @param legend.bg background colour for legend (default white overrides the
##'   grey lines, seems better)
##' @param font.main font for title, can't have bold with subscripts it seems,
##'   so set this for all plots.
##' @return single plot that explains one part of EDM, link together in a movie
##'   using `pbs_explain_edm_movie()` and `pbs_explain_edm_movie_save()`
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   aa <- pbsEDM(NY_lags_example,
##'               lags = list(N_t = 0:1),
##'               first_difference = TRUE)
##'   plot_explain_edm(aa,
##'                    tstar = 15,
##'                    main = paste0(
##'                            "See where neighbours go"),
##'                    tstar.pch = 19,
##'                    tstar.cex = 1.2,
##'                    neigh.plot = TRUE,
##'                    neigh.proj = TRUE)
##' # Type plot_explain_edm_movie to see other examples
##' }
plot_explain_edm <- function(obj,
                             tstar,
                             x.lab = expression(italic(Y) [t-1]),
                             y.lab = expression(italic(Y) [t]),
                             main = "Plot all the points in lagged space.",
                             tstar.col = "blue",
                             tstar.pch = 1,
                             tstar.cex = 1,
                             lib.remove = FALSE,
                             neigh.plot = FALSE,
                             neigh.proj = FALSE,
                             pred.plot = FALSE,
                             pred.rEDM = FALSE,
                             true.val = FALSE,
                             legend.plot = TRUE,
                             legend.inc.rEDM = TRUE,
                             legend.bg = "white",
                             font.main = 1){
  if(pred.rEDM){
    testthat::expect_equal(obj$X_observed, NY_lags_example$Y_t)
    #  ideally want this message:
    #  stop("pred.rEDM can only be TRUE when plotting NY_lags_example")
  }

  par(pty="s")
                             # TODO maybe. Andy NOT CHANGED THESE (yet) in notation update
  Xt <- obj$X[, "N_t_0"]      # Should change Nt in pbsEDM() as confusing. Luke
                              # was still tweaking that. Think here want Xt ->
                              # Y_t, and Xtmin1 -> X_tmin1
  Xtmin1 <- obj$X[, "N_t_1"]

  Xt.max.abs = max(abs( range( c(Xt,
                                 Xtmin1,
                                 obj$X_forecast),
                              na.rm=TRUE)))
  Xt.axes.range = c(-Xt.max.abs, Xt.max.abs)

  # Plot all the points
  plot(Xtmin1,
       Xt,
       xlab = x.lab,
       ylab = y.lab,
       xlim = Xt.axes.range,
       ylim = Xt.axes.range,
       main = main,
       font.main = font.main)

  # Highlight x[tstar]
  points(Xtmin1[tstar],
         Xt[tstar],
         col = tstar.col,
         pch = tstar.pch,
         cex = tstar.cex)

  if(lib.remove){
    lib.remove.ind <- c(1,       # points 1 and T aren't plotted anyway since
                                 # not defined
                        tstar + 1,
                        min(tstar + 2, nrow(obj$X)),
                        nrow(obj$X) - 1,
                        nrow(obj$X))
    lib.remove.ind <- lib.remove.ind[ - which(lib.remove.ind == tstar)]   # don't want to cross out tstar
                                                                          # (happens for tstar = T-1 for forecast)
    points(Xtmin1[lib.remove.ind],
           Xt[lib.remove.ind],
           pch = 4,
           cex = 2,
           col = "blue")
    }

  psi_vec <- obj$neighbour_index[tstar,]     #  vector of indices of points closest
                                        #  to x[tstar]

  # Highlight neighbours
  if(neigh.plot){
    points(Xtmin1[psi_vec],
           Xt[psi_vec],
           pch = 19,
           col = "red")
  }

  # Highlight projections of neighbours
  if(neigh.proj){
    points(Xtmin1[psi_vec + 1],
           Xt[psi_vec + 1],
           pch = 19,
           col = "purple",
           cex = 0.5)

    for(i in 1:length(psi_vec)){
      igraph:::igraph.Arrows(Xtmin1[psi_vec[i]],
                              Xt[psi_vec[i]],
                              Xtmin1[psi_vec[i] + 1],
                              Xt[psi_vec[i] + 1],
                              curve = 1,
                              size = 0.7,
                              h.lwd = 2,
                              sh.lwd = 2,
                              width = 1,
                              sh.col = "purple")
               # Example:
               # iArrows(0, 0, 4, 4, curve = 1, size = 0.7, h.lwd = 2,
               #         sh.lwd=2, width = 1, sh.col="red")
    }
  }

  # plot forecast value
  if(pred.plot){
    points(Xt[tstar],
           obj$X_forecast[tstar + 1],    # CHECK not tstar+1, think it's correct
           col = tstar.col,
           pch = 8)

    igraph:::igraph.Arrows(Xtmin1[tstar],
                           Xt[tstar],
                           Xt[tstar],
                           obj$X_forecast[tstar + 1],
                           curve = 1,
                           size = 0.7,
                           h.lwd = 2,
                           sh.lwd = 2,
                           width = 1,
                           sh.col = "darkgrey")

    abline(h = Xt[psi_vec + 1],
           col = "lightgrey")
  }

  # show predicted value from rEDM
  if(pred.rEDM){
    points(Xt[tstar],
           dplyr::pull(NY_lags_example[tstar+1, "rEDM.pred"]),
           col = "red",
           cex = 1.5,
           lwd = 2)
  }

  # plot the true value of xt(t^*+1)
  if(true.val){
    points(Xtmin1[tstar + 1],
           Xt[tstar + 1],
           cex = 1.5,
           lwd = 2,
           col = "darkgreen")
  }

  if(legend.plot){
    if(legend.inc.rEDM){
      legend("bottomleft",
             pch=c(tstar.pch, 19, 8, 1, 1),
             leg=c(expression(italic(t) * "*"),
                   "neighbours",
                   expression(italic(t) * "*+1 prediction"),
                   "rEDM pred",
                   expression("true " * italic(t) * "*+1")),
             col=c(tstar.col,
                   "red",
                   tstar.col,
                   "red",
                   "darkgreen"),
             cex=0.85,
             bg = legend.bg)} else {
                        legend("bottomleft",
                               pch=c(tstar.pch, 19, 8, 1),
                               leg=c(expression(italic(t) * "*"),
                                     "neighbours",
                                     expression(italic(t) * "*+1 prediction"),
                                     expression("true " * italic(t) * "*+1")),
                               # This was original with bold x:
                               # leg=c(expression(bold(x)* "(" * italic(t) * "*)"),
                               #       "neighbours",
                               #       expression(bold(x)* "(" *
                               #                            italic(t) *
                               #                            "*+1) pred (wt avge)"),
                               #       expression("true " * bold(x)* "(" * italic(t) * "*+1)")),
                               col=c(tstar.col,
                                     "red",
                                     tstar.col,
                                     "darkgreen"),
                               cex=0.85,
                               bg = legend.bg)
                      }

    legend("topleft",
           leg = as.expression(bquote(italic(t)* "*=" *.(tstar))),
           col = "white",
           bty = "n")
  }
}

##' Create and save encapsulated postscript file for manuscript to explain EDM.
##'
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param eps.filename filename to save to
##' @param tstar focal point to make projection from
##' @return saved .eps file to use in manuscript.
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##' # For manuscript:
##' E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##' plot_explain_edm_save(E_results[[1]])
##' }
plot_explain_edm_save <- function(E_res,
                                  eps.filename = "explain_EDM.eps",
                                  tstar = 39){
  postscript(eps.filename,
             height = 6.5,
             width = 6.5,
             horizontal=FALSE,
             paper="special")

  plot_explain_edm(E_res,
                   tstar = tstar,
                   main = "",
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   neigh.plot = TRUE,
                   neigh.proj = TRUE,
                   pred.plot = TRUE,
                   pred.rEDM = FALSE,
                   true.val = TRUE,
                   legend.inc.rEDM = FALSE)
  dev.off()
}


##' Create movie to explain EDM
##'
##' Use with gifski in .Rmd - see `analyse_simple_time_series` vignette. Also
##'  called from `plot_explain_edm_movie_save` to save a pdf that can be
##'  incremented through in the Appendix.
##'
##' @param obj list object of class `pbsEDM`, an output from `pbsEDM()`
##' @param inc.rEDM logical whether to include rEDM calculation figure and in legend.
##' @param ... extra values to use in calls to `plot_explain_edm()`; needs to include tstar
##' @return movie (if used with gifski) to explain EDM; if combined with
##'   `par(mfrow=c(4,2))` (and not gifski) then will give multiple panels, but
##'   borders etc have not been set up properly for this and will need adjusting
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   aa <- pbsEDM(NY_lags_example,
##'               lags = list(N_t = 0:1),
##'               first_difference = TRUE)
##'   plot_explain_edm_movie(aa, tstar = 15) # will only show last panel; see
##'                                          # `analyse_simple_time_series`
##'                                          #  vignette.
##' }
plot_explain_edm_movie <- function(obj,
                                   inc.rEDM = FALSE,
                                   ...){

  # Need tstar here
  # https://stackoverflow.com/questions/64851094/extract-value-of-argument-from-ellipsis-without-evaluating-other-arguments
  if(hasArg(tstar)) {
    tstar <- eval.parent(match.call()[["tstar"]])
  }

  plot_explain_edm(obj,
                   tstar.col = "black",
                   legend.inc.rEDM = inc.rEDM,
                   legend.plot = FALSE,
                   ...)    # includes tstar that  gets carried through

  plot_explain_edm(obj,
                   main = as.expression(bquote("Want to predict" * italic(Y) [ .(tstar+1)] *
                                        " by locating (" *
                                        italic(Y) [ .(tstar -1)] * ", " *
                                        italic(Y) [ .(tstar)] * ").")),
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   legend.inc.rEDM = inc.rEDM,
                   ...)

  # Remove all unusable potential neighbours
  plot_explain_edm(obj,
                   tstar.col = "blue",,
                   main = paste0(
                     "Exclude invalid potential neighbours ..."),
                   legend.inc.rEDM = inc.rEDM,
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   lib.remove = TRUE,
                   ...)

  plot_explain_edm(obj,
                   main = paste0(
                     "... and determine the three nearest valid neighbours."),
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   neigh.plot = TRUE,
                   legend.inc.rEDM = inc.rEDM,
                   lib.remove = TRUE,
                   ...)

  plot_explain_edm(obj,
                   main = paste0(
                     "See where they get projected to in their next time step ..."),
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   neigh.plot = TRUE,
                   neigh.proj = TRUE,
                   legend.inc.rEDM = inc.rEDM,
                   ...)

  plot_explain_edm(obj,
                   main =
                   as.expression(bquote("... and take a weighted average of " *
                                        italic(Y) [t] * " to give " *
                                        italic(Y) [ .(tstar+1)])),
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   neigh.plot = TRUE,
                   neigh.proj = TRUE,
                   pred.plot = TRUE,
                   legend.inc.rEDM = inc.rEDM,
                   ...)

  if(inc.rEDM){
    plot_explain_edm(obj,
                     main = paste0("... should agree with prediction from rEDM and ..."),
                     tstar.pch = 19,
                     tstar.cex = 1.2,
                     neigh.plot = TRUE,
                     neigh.proj = TRUE,
                     pred.plot = TRUE,
                     pred.rEDM = TRUE,
                     legend.inc.rEDM = inc.rEDM,
                     ...)
  }

  plot_explain_edm(obj,
                   main = paste0(
                     "... which is hopefully close to the true known value."),
                   tstar.pch = 19,
                   tstar.cex = 1.2,
                   neigh.plot = TRUE,
                   neigh.proj = TRUE,
                   pred.plot = TRUE,
                   pred.rEDM = TRUE,
                   true.val = TRUE,
                   legend.inc.rEDM = inc.rEDM,
                   ...)
}




##' Make pdf movie of explain EDM figures to go into Supp pdf.
##'
##' Creates a single .pdf with one frame for each figure,
##'  included as a pausable movie of Supp pdf for primer manuscript.
##'  See `analyse_simple_time_series.Rmd` vignette for creating the similar .gif
##'  that is shown in the vignette.
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param pdf.filename filename to save to
##' @param tstar focal point to make projection from
##' @param ... extra arguments to pass to `plot_explain_edm_movie()`
##' @return saved pdf file with one page for each figure, that can become
##'  movie in Supp of primer manuscript.
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##'   plot_explain_edm_movie_save(E_results[[1]])
##' }
plot_explain_edm_movie_save <- function(E_res,
                                        pdf.filename = "explain_EDM_movie.pdf",
                                        tstar = 39,
                                        ...){
  pdf(pdf.filename,
      height = 6,
      width = 5.36)

  plot_explain_edm_movie(E_res,
                         tstar = tstar)

  dev.off()
}

##' Make pdf movie of explain EDM figure for each tstar to go into Supp pdf.
##'
##' Creates a single .pdf with one frame for each figure (value of tstar),
##'  included as a pausable movie of Supp pdf for primer manuscript.
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param pdf.filename filename to save to
##' @param ... extra arguments to pass to `plot_explain_edm_movie()`
##' @return saved pdf file with one page for each figure, that can become
##'  movie in Supp of primer manuscript.
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##'   plot_explain_edm_all_tstar_movie_save(E_results[[1]])
##' }
plot_explain_edm_all_tstar_movie_save <- function(E_res,
                                                  pdf.filename =
                                                    "explain_EDM_all_tstar_movie.pdf",
                                                  ...){
  pdf(pdf.filename,
      height = 6,
      width = 5.36)

  for(tstar_val in 2:(nrow(E_res$X) - 1)){    # For E=2
    plot_explain_edm(E_res,
                     tstar = tstar_val,
                     main = "",
                     tstar.pch = 19,
                     tstar.cex = 1.2,
                     neigh.plot = TRUE,
                     neigh.proj = TRUE,
                     pred.plot = TRUE,
                     pred.rEDM = TRUE,
                     true.val = TRUE,
                     legend.inc.rEDM = TRUE)
  }
  dev.off()
}

##' Make pdf movie of six-panel figure to go into Supp pdf.
##'
##' Creates a single .pdf with one frame for each time step, which can be
##'  included as a pausable movie of Supp pdf for primer manuscript.
##'  See `analyse_simple_time_series.Rmd` vignette for creating the similar .gif
##'  that is shown in the vignette.
##' @param E_res List of `pbsEDM` objects as output from `pbsEDM_Evec()`
##' @param pdf.filename filename to save to
##' @param ... extra arguments to pass to `plot.pbsEDM()`
##' @return saved pdf file with one page for each time step, that can become
##'  movie in Supp of primer manuscript.
##' @export
##' @author Andrew Edwards
##' @examples
##' \donttest{
##'   E_results <- pbsEDM_Evec(NY_lags_example$N_t)
##'   plot_pbsEDM_Evec_movie_save(E_results)
##' }
plot_pbsEDM_Evec_movie_save <- function(E_res,
                                        pdf.filename = "six_panel_movie.pdf",
                                        ...){
  pdf(pdf.filename,
      height = 5.36,
      width = 9)

  num.frames <- length(E_res[[1]]$X_observed)

  for(i in 1:num.frames){
    plot_pbsEDM_Evec(E_res = E_res,
                     last.time.to.plot = i,
                     portrait = FALSE)
  }
  dev.off()
}


#' Plot Method for `pbsSim`
#'
#' @param x [pbsSim()]
#' @param ... Other arguments.
#'
#' @return Panel plot.
#' @export
#'
plot.pbsSim <- function (x, ...) {
  par(mfrow = c(3, 1), mar = c(4, 4, 1, 0), oma = c(2, 1, 1, 1))
  plot(x[, "producers"], type = "l", ylab = "Producers")
  plot(x[, "prey"], type = "l", ylab = "Prey")
  plot(
    x[, "predators"],
    type = "l",
    xlab = "Time step",
    ylab = "Predators")
  par(mfrow = c(1, 1))

}

##' Plot the size of the library for a given `T` as a function of `E` and `tstar`
##'
##' Show a coloured contour plot of how the library size changes for a
##' univariate time series with sample size `T` as chosen embedding dimension
##' `E` and target index `tstar` vary.
##'
##' @param T integer of the sample size of a univariate time series (time series
##'   itself not needed)
##' @param E_vec vector of integers to show for the embedding dimension
##' @param tstar_vec vector of integers to show for the focal point time index
##'   `tstar`, `t*` in the manuscript, from which projections are made from
##' @param annotate logical whether to add numbers along the middle of the plot
##'   (avoids the need for a colorbar, which was fiddly to do, see attempts
##'   deleted in 8f567ff, while also having a grey background), and add text to
##'   explain the grey areas.
##' @param annotate_cex text size for main annotation
##' @param annotate_tstar tstar value at which to add the annotated numbers; if
##'   NULL then is `max(tstar_vec)/2`.
##' @param annotate_extra_cex text size for extra smaller numbers in boxes at
##'   the top
##' @param words_pos vector of positioning for words in grey areas, given by `E`
##'   and `tstar` position for `Early values' text (text is placed to the right
##'   of these), then for `Late values' text (text is centred around these.
##' @param words_cex text size for above words
##' @param mgp_vals mgp values to use, as in `plot(..., mgp = mgp_vals)`.
##' @return plots the contour plot and returns the matrix (`C` in manuscript) of calculated library sizes
##' @export
##' @author Andrew Edwards
##' @examples
##' \dontrun{
##' # For manusript doing:
##' postscript("library_size.eps",
##'             height = 6,
##'             width = 6,
##'             horizontal=FALSE,
##'             paper="special")
##' plot_library_size()
##' dev.off()
##' }
plot_library_size <- function(T = 50,
                              E_vec = 2:10,
                              tstar_vec = 1:50,
                              annotate = TRUE,
                              annotate_cex = 1,
                              annotate_tstar = 23,
                              annotate_extra_cex = 0.7,
                              words_pos = c(6, 3.5, 6, 49.6),
                              words_cex = 0.9,
                              mgp_vals = c(2, 0.5, 0)){
  stopifnot(length(T) == 1,
            length(E_vec) > 1,
            min(E_vec) > 1,
            length(tstar_vec) > 1)

  stopifnot("T is too small relative to max(E_vec); change code if you want to try exceptions" =
              T - 2 * (max(E_vec) + 1) >= 0)

  if(is.null(annotate_tstar)){
    annotate_tstar <- max(tstar_vec)/2
  }

  C <- matrix(NA,
              nrow = length(E_vec),
              ncol = length(tstar_vec))

  for(i in 1:length(E_vec)){
    E <- E_vec[i]
    if(T - E - 2 >= E){
      for(tstar in E:(T - E - 2)){
        j <- which(tstar_vec == tstar)
        C[i, j] <- T - 2 * (E + 1)
      }
    }
    for(tstar in (T - E - 1):(T - 2)){    # E > 1 so always valid
      j <- which(tstar_vec == tstar)
      C[i, j] <- tstar - E
    }
    # And tstar <= E-1 and tstar >= T-1 remain as NA
  }

  image(E_vec,
        tstar_vec,
        C,
        xlim = range(E_vec) + c(-0.5, 0.5),   # colours are correctly around integers
        ylim = rev(range(tstar_vec) + c(-0.5, 0.5)),
        xaxs = "i",
        yaxs = "i",         # sets exact axis range
        xlab = expression(paste("Embedding dimension, ", italic(E))),
        ylab = expression(paste("Focal time, ", italic(t), "*")),
        mgp = mgp_vals)

  rect(0,
       0,
       max(E_vec) + 2,
       max(tstar_vec) + 2,
       col = "darkgrey")    # So NA's come out grey

  image(E_vec,
        tstar_vec,
        C,
        add = TRUE,
        col = rev(hcl.colors(max(C, na.rm=TRUE) - min(C, na.rm=TRUE) + 1,
                         "Spectral")))

  # Add extra unlaballed tickmarks
  box()
  axis(1, E_vec, tcl = -0.4, labels = rep("", length(E_vec)))
  axis(2, tstar_vec, tcl = -0.2, labels = rep("", length(tstar_vec)))
  every_five <- tstar_vec[which(tstar_vec %% 5 == 0)]
  axis(2, every_five, tcl = -0.4, labels = rep("", length(every_five)))

  if(annotate){
    annotate_main <- C[1:length(E_vec),
                       which(tstar_vec == annotate_tstar)]  # main annotation
    text(E_vec,
         annotate_tstar,
         annotate_main,
         cex = annotate_cex)

    annotate_extra_three <- which(C > max(annotate_main),
                                  arr.ind = TRUE)          # extra small annotations
                                                           # (just three for default)


    tstar_max <- max(annotate_extra_three[ , 2]) # max valid tstar, same for all
                                               #  E, namely (T-2 = 48 for default)
    # Want value at t^*=48 and values along diagonal
    annotate_extra_row <- 1:length(E_vec)

    annotate_extra_max_tstar_col <- rep(tstar_max, length(E_vec))  # values at t^* = 48
    annotate_extra_diag_col <- tstar_max - 1:length(E_vec)         # values along diagonal

    annotate_extra <- rbind(annotate_extra_three,                  # duplicates some
                                                                   #  values but keeps the
                            cbind(annotate_extra_row,              #  col headings
                                  annotate_extra_max_tstar_col),
                            cbind(annotate_extra_row,
                                  annotate_extra_diag_col))

    for(k in 1:nrow(annotate_extra)){
      text(E_vec[annotate_extra[k, "row"]],
           tstar_vec[annotate_extra[k, "col"]],
           C[annotate_extra[k, "row"],
             annotate_extra[k, "col"]],
           cex = annotate_extra_cex)
    }

    text(words_pos[1],
         words_pos[2],
         expression(paste("Early values of ",
                          italic(t),
                          "* not possible")),
         pos = 4,
         cex = words_cex)
    text(words_pos[3],
         words_pos[4],
         expression(paste("Late values of ",
                          italic(t),
                          "* not possible")),
         cex = words_cex)
  }

  return(C)
}

##' Save eps figure of the size of the library for manuscript.
##'
##' @param eps.filename filename to save to
##' @return saved .eps file to use in manuscript.
##' @export
##' @author Andrew Edwards
plot_library_size_save <- function(eps.filename = "library_size.eps"){
  postscript(eps.filename,
             height = 6,
             width = 6,
             horizontal=FALSE,
             paper="special")

  plot_library_size()

  dev.off()
}